Spinal Mechanisms May Provide a Combination of Intermittent and Continuous Control of Human Posture: Predictions from a Biologically Based Neuromusculoskeletal Model

Several models have been employed to study human postural control during upright quiet stance. Most have adopted an inverted pendulum approximation to the standing human and theoretical models to account for the neural feedback necessary to keep balance. The present study adds to the previous efforts in focusing more closely on modelling the physiological mechanisms of important elements associated with the control of human posture. This paper studies neuromuscular mechanisms behind upright stance control by means of a biologically based large-scale neuromusculoskeletal (NMS) model. It encompasses: i) conductance-based spinal neuron models (motor neurons and interneurons); ii) muscle proprioceptor models (spindle and Golgi tendon organ) providing sensory afferent feedback; iii) Hill-type muscle models of the leg plantar and dorsiflexors; and iv) an inverted pendulum model for the body biomechanics during upright stance. The motor neuron pools are driven by stochastic spike trains. Simulation results showed that the neuromechanical outputs generated by the NMS model resemble experimental data from subjects standing on a stable surface. Interesting findings were that: i) an intermittent pattern of muscle activation emerged from this posture control model for two of the leg muscles (Medial and Lateral Gastrocnemius); and ii) the Soleus muscle was mostly activated in a continuous manner. These results suggest that the spinal cord anatomy and neurophysiology (e.g., motor unit types, synaptic connectivities, ordered recruitment), along with the modulation of afferent activity, may account for the mixture of intermittent and continuous control that has been a subject of debate in recent studies on postural control. Another finding was the occurrence of the so-called “paradoxical” behaviour of muscle fibre lengths as a function of postural sway. The simulations confirmed previous conjectures that reciprocal inhibition is possibly contributing to this effect, but on the other hand showed that this effect may arise without any anticipatory neural control mechanism.


Introduction
The maintenance of upright quiet stance is a challenging task for the central nervous system (CNS), and the objective is to achieve the control of an intrinsically unstable biomechanical system under the effect of gravity. Posture control is a position control problem in which the CNS, leg muscles, and different sensory systems (e.g., the muscle proprioceptors) interact to maintain the horizontal projection of the centre of mass (COM) within a region bounded by the feet (for a review of the basics of posture control, see [1]). Sensory systems, such as the vestibular, visual, and somatosensory, play a significant role in the aforementioned motor task, so that disorders in any of these systems may lead to postural instability [2,3].
A conceptual question under debate in the literature concerns the manner the CNS controls upright stance in adults (i.e., human beings who already learned to walk). Some researchers argue in favour of a negative-feedback continuous control, with the leg muscles reflexively activated in response to drifts or perturbations away from an equilibrium position [4][5][6]. Conversely, others suggest that an anticipatory (feedforward) mechanism is necessary to explain some findings, such as the low feedback gains observed in some experiments [7][8][9]. More recently, the view that human upright stance is controlled by an intermittent mechanism has grown [10][11][12][13][14][15]. The latter is based in part on results from an experiment involving the manual balancing of a virtual inverted pendulum controlled by a subject through a computer joystick [13]. Additionally, other recent studies showed that motor units (MUs) of the Medial Gastrocnemius (MG) muscle exhibited an intermittent recruitment during upright quiet stance with a pattern closely linked to COM and centre of pressure (COP) fluctuations [14,16]. It is noteworthy that the previously referred papers used a wide-sense/physiological conceptualisation of intermittent postural control, which is understood as one exhibiting discrete (bursty) actions of the neuronal controller, producing ballistic-like (phasic) muscle activations. In addition to this more qualitative/physiological view of intermittency, there is a strict, quantitative, definition of intermittent control that has also been used to analyse motor control problems, such as stick balancing and postural control [17]. In the present study we will adopt the former (physiological) view of intermittency.
From a theoretical standpoint, mathematical models have been developed to describe and to investigate human postural control. These models explicitly incorporate the above-mentioned control strategies hypothesised to be adopted by the CNS during upright standing control, but do not identify which part of the CNS the control resides, e.g. cortex, brainstem, cerebellum, spinal cord. For instance, [4,5,[18][19][20] represented the postural control system as a negative-feedback continuous control system, so that proprioceptive information of body position and velocity were fed back to the CNS. Alternatively, the models proposed by [21][22][23] adopted a continuous predictive control system, while others [11,12,24,25] represented the control of posture as an intermittent control system. Despite these fundamental differences, all these models were based on a control engineering framework, whereby the whole system was simplified and the CNS was represented by a PD/PID (proportional-derivative/proportional-integral-derivative) controller, an optimal continuous controller or an intermittent controller that activated the muscles in discrete bursts. Despite their advantages of a relative simplicity and the power of explanation, it is not easy to translate their results in terms of the underlying physiological mechanisms involved in the control of human posture.
The present study aims at providing a complementary approach to those mentioned above in that the fundamental focus is on biology (anatomy and physiology). A large amount of physiological knowledge, from the behaviour of neuronal ionic channels to the dynamics of muscle contraction, was funnelled into a large-scale mathematical model of the nervous, muscular, and biomechanical systems involved in posture control. In this vein, a biologically based large-scale neuromusculoskeletal (NMS) model was developed and used to investigate the problem of postural control from a more neurophysiological standpoint. The model encompasses: i) a spinal neuronal network, which includes conductance-based models of both a motor neurons (MNs) and interneurons (INs); ii) Hill-type muscle models to represent the viscoelastic properties of the Soleus (SO), MG, Lateral Gastrocnemius (LG), and Tibialis Anterior (TA) muscles; iii) models of both muscle spindle and Golgi tendon organ (GTO); iv) afferent fibres providing Ia, II, and Ib feedback; and v) an inverted pendulum model, which is a first approximation of upright quiet stance [26]. Here the proprioceptive feedback is provided by muscular proprioceptors (i.e., spindles and GTOs), since they seem to be largely responsible for the position sense of the limb [27,28]. The first hypothesis of the present study is that stance control might be properly achieved by a spinal-like controller (SLC, approximated here by the developed NMS model) based on a proprioceptive feedback. An associated second hypothesis is that the activation of leg muscles by this SLC is a continuous process yielding the maintenance of human body equilibrium.
Another relevant issue related to postural control that the present study analyses is related to recent experimental findings of a ''paradoxical'' behaviour of the calf's muscle fibre lengths during postural sway [29,30]. The behaviour was called ''paradoxical'' because the muscle fibre lengths were negatively correlated with COM/COP displacements [31]. These studies proposed that reciprocal inhibition from antagonistic (TA) Ia afferents might be responsible for this unexpected motor behaviour. With the availability of the detailed NMS model employed in the present study, the correlation between muscle fibre length and COM/ COP displacements could be tackled without much difficulty. Therefore, a third hypothesis of the present study is that a model structure without the reciprocal inhibition pathway (from TA Ia afferents to Triceps Surae MNs) would exhibit positive correlation coefficients between muscle fibre length and COM/COP displacements. A complement of this hypothesis is that the reciprocal inhibition neural circuit may increase the probability of occurrence of the ''paradoxical'' relation mentioned above. This combined third hypothesis is, therefore, evaluated here using the complex NMS physiological model in order to verify a hypothesis put forward by previous authors on the basis of heuristics and experimental results from humans [29][30][31].
To the best of our knowledge, this is the first study to address the control of an intrinsically unstable neuro-biomechanical system associated with the maintenance of human quiet standing by means of a complex large-scale system mostly based on known physiology. However, others have also used biologically based reductionist models of the NMS system to investigate how the CNS controls other motor tasks [32,33]. A minor part of this material was already published as conference abstracts [34,35].

Model Validation in Terms of Postural Control
Typical biomechanical and neuronal outputs of the NMS model are presented in Figure 1. The model's responses resemble qualitatively those frequently reported in postural control studies (e.g., [9,29,36]). Irrespective of model structure (i.e., Model 1 or Model 2 -see Methods for details), the inverted pendulum leaned about 5 deg forward (equilibrium point), so that COM and COP displacements oscillated around a basal value of 80 mm (see Figure 1A). The basal plantar-flexion torque (negative torque) was *10% of the maximum isometric torque produced by the model. One can notice that COM and COP ( Figure 1A) oscillated in anti-

Author Summary
The control of upright stance is a challenging task since the objective is to maintain the equilibrium of an intrinsically unstable biomechanical system. Somatosensory information is used by the central nervous system to modulate muscle contraction, which prevents the body from falling. While the visual and vestibular systems also provide important additional sensory information, a human being with only somatosensory inputs is able to maintain an upright stance. In this study, we used a biologically-based large-scale neuromusculoskeletal model driven only by somatosensory feedback to investigate human postural control from a neurophysiological point of view. No neural structures above the spinal cord were included in the model. The results showed that the model based on a spinal control of posture can reproduce several neuromechanical outcomes previously reported in the literature, including an intermittent muscle activation. Since this intermittent muscular recruitment is an emergent property of this spinal-like controller, we argue that the so-called intermittent control of upright stance might be produced by an interplay between spinal cord properties and modulated sensory inflow.
phase with respect to the muscle torque ( Figure 1B), i.e. when the body leaned forward from its equilibrium position the plantarflexion torque increased (more negative). Conversely, muscle activations (EMG envelopes in Figure 1C-E) were modulated approximately in phase with postural sway. In the simulations, TA muscle was silent during postural sway (not shown).
A quantitative analysis was performed to validate the model with respect to the available data from the literature. Typical timedomain metrics were calculated from the COP time series and compared to data from normal subjects and vestibular loss patients standing on a force plate without visual information (see Table 1). Root mean square (RMS) and mean velocity (MV) of simulated COP were higher than the values observed experimentally in normal subjects, but compatible with data from vestibular loss patients. Another quantitative validation was based on a crosscorrelation analysis performed between the COM and COP time series (Figure 2A-B), as well as between COP and EMG envelopes ( Figure 2C-D). COM and COP were highly correlated (r&1) at lag zero. COP and EMG envelopes were positively correlated with the correlation peak occurring at a positive lag. Correlation coefficients (r) and cross-correlation peak lag values were compatible with experimental data from healthy subjects (see Table 1). In general, correlation coefficients were higher for Gastrocnemii in comparison to the SO, and muscles' activations (EMGs) were advanced by approximately 200-300 ms in relation to the postural sway (COP). The 50% power frequency (F 50) estimated from the COP power spectrum (see Figure 2E-F) resulted quite similar to the value from healthy subjects (see Table 1). COP power spectra of both model structures were limited to *1 Hz.
A final quantitative validation was based on the pooled histogram of COM displacements (1-mm bins) as shown in Figure 3 (data are from the simulations of Model 2). The histogram shape was bimodal, with two peaks around the equilibrium position of the inverted pendulum (value 0 in the abscissa). The Jarque-Bera goodness-of-fit test was applied to verify if this histogram could be fitted by a typical Gaussian probability density function [11]. The null-hypothesis (the histogram comes from an unimodal Gaussian function) was rejected (p~0:001). The same result was obtained for Model 1. INs, and afferent fibres were modulated during postural sway. An interesting qualitative finding was that MUs from the MG muscle were intermittently recruited/de-recruited as the inverted pendulum swayed forward/backward ( Figure 4B). This intermittent pattern of MU recruitment was similar for the LG muscle (not shown), but less evident for the SO muscle (see Figure 5A). The degree of intermittency for the MG and SO MUs was quantified  . Cross-correlation functions and centre of pressure (COP) power spectra for typical simulations carried out on Model 1 (graphs A, C, and E) and Model 2 (graphs B, D, and F). (A-B) Cross-correlation functions between centre of mass (COM) and COP. Note that for both models, cross-correlation peaks occurred at zero lag (dashed lines). (C-D) Cross-correlation functions between COP and muscle electromyograms (EMGs). Black, red, and blue curves represent cross-correlation functions for Soleus (SO), Medial Gastrocnemius (MG), and Lateral Gastrocnemius (LG), respectively. Irrespective of the model structure, there was a lag of about 300 ms between COP and EMG envelopes from the three muscles. (E-F) COP power spectra. Dashed line represents the 50% power frequency (F 50). It is noteworthy that for Model 2 there was a broader bandwidth in comparison to Model 1. doi:10.1371/journal.pcbi.1003944.g002 In order to quantify the intermittent recruitment of MG MUs, the interval between successive recruitments was computed for a subset of 30 randomly chosen MUs (10 MUs were chosen per simulation). In accordance with the method used by [14], intermittent recruitment was considered if a given MU was discharging at a rate lower than 4 Hz (i.e., interspike intervals higher than 250 ms LG MUs was recruited (less than 30) and the SO MUs were mostly tonically active during the simulation of postural control (see the activation ratios in previous paragraph), hence the intermittency of the MUs from these muscles were not quantitatively evaluated here.

Intermittent Recruitment of the Motor Units
Panels C-I in Figure 4 and panel B in Figure 5 show typical results of how proprioceptive feedback (encompassing afferent fibres and spinal INs) was modulated during sway (Model 2 was used for this simulation). The activity of the Ia afferents from the MG muscle ( Figure 4C) was highly modulated, following approx-   imately the COM/COP displacement (note the firing rate modulation for three different Ia afferents, as indicated by the thin lines). Since there was little variation in the MG muscle torque (RMS value *2.50% of the maximum MG muscle torque) and the mean MG muscle fibre length was maintained at an approximately steady value (i.e., there was little change in the static component of the muscle fibre length -see below), the activities of Ib and type II afferents were minimally modulated (see panels G and H in Figure 4). The proprioceptive pathways responsible for the reciprocal inhibition were also highly modulated during postural sway (see panels F and I in Figure 4). Inhibitory Ia INs discharged phasically when the inverted pendulum swayed backward and this contributed to a decrease in the ankle joint torque generated by the plantar flexor muscles. Conversely, Ia afferents from the SO muscle spindles were poorly modulated in the posture control task (see Figure 5B).
The intermittent recruitment of MG MUs was evaluated on the basis of two phase plots that relate angular velocity and muscle torque with ankle angle data obtained from the postural control model ( Figure 6). Figure 6A shows that most of the MG MUs (*60%) were recruited when the inverted pendulum was leaning forward from its equilibrium position irrespective of its velocity (first and fourth quadrants of the angle-velocity phase plots). Nonetheless, a large number of MUs (*28%) were recruited when the inverted pendulum was at a backward position but with a positive velocity (second quadrant), i.e., the pendulum was starting to return to a forward position. Similarly, most of the MG MUs (*50%) were recruited when the pendulum was leaning forward and producing a higher plantar flexion torque (fourth quadrant in the Figure 6B), i.e., the pendulum was at a forward position and decelerating. Almost 35% of the MUs were recruited when the pendulum was at a backward position and with a lower (more positive than the mean value) plantar flexion torque. In general, the discharges in the first and third quadrants were mainly involved in the generation of a basal torque, while the discharges in the second and fourth quadrants represented the phasic corrective torque control produced by the MG muscle.

Effect of Reciprocal Inhibition
Two model structures were adopted to investigate the effect of the reciprocal inhibition pathway on postural control (see Methods for details). According to the data presented in the previous section, there were no dramatic differences in the time-and frequency-domain metrics, as well as in the MU recruitments. The COP power spectrum calculated from the model without reciprocal inhibition (Model 1) was slightly narrower than that from Model 2. In addition, correlation coefficients between COP and EMG envelopes were lower when the reciprocal inhibition pathway was included in the model. Nonetheless, both model . It is noteworthy that MUs seem to be preferentially recruited in the first quadrant, i.e., when the pendulum is leaning forward as it comes from a backward position. (B) Phase plot showing the recruitment of MG MUs in angle-torque bi-dimensional graph. Again, the point (0,0) is a reference that was used because each simulation produced slightly different mean torques and angles. Note that MG MUs are mostly activated in the fourth quadrant, i.e., when the pendulum is leaning forward and the plantar-flexion torque produced by Triceps Surae (TS) is higher than the mean. doi:10.1371/journal.pcbi.1003944.g006 structures produced fluctuating outputs that resembled experimental data, suggesting that reciprocal inhibition from TA Ia afferents is not a strict requisite for the control of upright standing.
In the present study, we tested the hypothesis that reciprocal inhibition may be responsible for the negative correlation between muscle fibre length and COM/COP displacement [29][30][31]. This was tested by performing a correlation analysis between COM displacements and muscle fibre lengths for the SO, MG, and TA. Typical signals representing these variables are shown in Figure 7. For this typical simulation (Model 2) one can notice that MG muscle fibre length was positively correlated with COM, while TA muscle fibre length and COM displacement were negatively correlated. The latter results are a typical ''orthodox'' behaviour that has been shown for some healthy subjects during upright quiet standing [29]. For the SO muscle, a more quantitative analysis was performed (see Figure 8). Correlation analysis between SO muscle fibre length and COM displacement was performed on 3-s windows (see dashed vertical lines in Figure 7) according to the method adopted in [29]. Correlation coefficients calculated from Model 1 and Model 2 were pooled into two groups: positively correlated (rw0) and negatively correlated (rv0). For Model 1 most of the intervals (*80%) showed a positive correlation between SO muscle fibre length and COM displacement (i.e., ''orthodox'' behaviour), whereas for Model 2 the number of negatively correlated intervals increased to *50%. x 2 -test revealed a statistically significant difference (pv0:05) between the responses of the two model structures. This suggests that, at least for the SO muscle, reciprocal inhibition might be involved in the genesis of the ''paradoxical'' behaviour of muscle fibre length. Nonetheless, a small percentage of negatively correlated intervals (*20%) might also be generated in the absence of any reciprocal inhibition from antagonistic Ia afferents.

Discussion
A large-scale NMS model was applied in the present study to the problem of controlling upright standing in humans. A different feature of this approach in comparison with most of the previous studies in the literature (e.g., [4][5][6][10][11][12]18,19,[22][23][24]) is that the structure and behaviours of each element were based on known physiology, anatomy, and biomechanics encompassing important parts of the postural control system. Therefore, the control strategy employed by the modelled CNS emerged from the interaction between several neuromechanical elements involved. While the overall mechanisms that control the inverted pendulum sway are beyond an analytical understanding, the increased biological realism provides important clues regarding some putative neurophysiological mechanisms underlying the posture control task. In the following sections, the results presented earlier shall be discussed with respect to relevant experimental findings reported in the literature.

Postural Control in Quiet Standing
The results presented here showed that a SLC was effective in maintaining the equilibrium of an intrinsically unstable biomechanical system (see Figure 1 for an illustrative example). This leads to the first contribution of this study, which is to support the hypothesis that human upright quiet standing may be properly controlled by spinal mechanisms, for example, without any cortical involvement. This view is consistent with several studies, which suggest that cortical control is decreased or may be absent when a motor task is well trained (e.g., [30][31][32][33][34][35][36][37][38][39][40]). Normal quiet postural control, under no special restrictions (as standing on a narrow beam or during dual tasks), is certainly a candidate for a well-trained task that would not require cortical control. Since the scope of the present study is limited to the investigation of neurophysiological mechanisms underlying the control of quiet standing, the assumed lack of supraspinal neural structures should not limit the ensuing interpretations. A separate section below (see Model Limitations and Future Research) presents and discusses the limitations both of the modelling as well as the conclusions derived from the simulations.
Time-and frequency-domain metrics obtained from the NMS model were compatible with experimental data (see Table 1). The equilibrium values of the inverted pendulum (mean angle, torque, and COM/COP displacements) varied slightly both within and between simulations due to programmed (randomised) changes in system configuration. Between-simulation changes occurred due to changes in neuronal connectivities and intrinsic properties (e.g., action potential thresholds) that were randomly attributed at the beginning of each simulation [41]. On the other hand, withinsimulation variations were mainly related to the number of recruited MUs, which varied stochastically due to neuronal noise. Therefore, the mean equilibrium position of the body depends on the overall instantaneous configuration of the postural control system. The analysis of the COP time series showed that the model was less stable (i.e., presented a larger postural oscillation) than healthy young subjects [42]. Notwithstanding, simulated data were compatible with those from vestibular loss subjects standing on a stable surface without visual information [3]. As a consequence, the simulation results reinforce that the increased postural oscillation observed in patients may be due to the lack of other sensory inputs providing information to the CNS, such as vestibular and visual sources. Or, in other words, the proprioceptive feedback gain in such patients is not sufficient to replace the other missing sensory feedback modalities. Interestingly, the variability observed in the simulated postural sway was exclusively generated by the variability in sensory afferents and descending commands, which results in random fluctuations of motoneuronal discharges. Therefore, it is predicted that most of the biomechanical variability (sway) observed during upright standing has a neuronal origin and is less influenced by internal disturbing forces (e.g., heartbeats and respiration) as proposed elsewhere [4,6,8,11,25].
The cross-correlation analysis ( Figure 2) showed that EMGs from the Triceps Surae (TS) muscles were positively correlated with postural sway (as measured by the COP). Simulation results are in accordance with experimental data that showed higher correlation coefficients between EMGs from Gastrocnemii and COP [9,43]. Additionally, the time lags between COP and EMGs were within a range of 200-300 ms, which is also compatible with experimental data [9,43]. This is in some sense a remarkable result that emerged from the NMS model since the sum of the afferent and efferent action potential propagation delays is much smaller than this time lag between COP and EMG. Albeit qualitative, or semi-quantitative, this is a relatively strong sign that the model was able to capture at least a part of the overall system dynamics. In [9] it is argued that the existence of this lag between the mechanical and neuronal responses would be due to an anticipatory action of the neuronal controller, i.e., the postural control is mediated by a feedforward mechanism. However, theoretical and computer simulation studies [19] showed that even in the absence of any feedforward mechanism, lags between neuromechanical signals may be obtained depending on the parameters of the continuous feedback system and stochastic features of the input signals. The results presented here corroborate the latter view, since no feedforward mechanism was incorporated into the NMS model.
Another experimental finding from human postural control studies that was reproduced by the model was the bimodal distribution of COM displacements (see Figure 3). In [11] most of the data from young subjects exhibited double-peaked histograms of COM displacements with a local minimum in between (see their Figure 5). The authors of the above-mentioned study [11] showed that the bimodal distribution of COM displacements was only obtained when they represented the postural control system by a mixture of both continuous and intermittent (with a phasic controller operating at a 3-4 Hz burst rate) control mechanisms. A continuous postural control model produced unimodal Gaussian-like histograms. Therefore, they argued that the postural control in humans is not mediated exclusively by a continuous control mechanism. Our results corroborate this proposal. However, in our approach the control structure (e.g., continuous, intermittent, or a mixture of both) was not imposed a priori. The mixture of continuous (SO motor nucleus and muscle fibres) and intermittent (mostly Gastrocnemii motor nuclei and muscle fibres) control behaviours was a result of the interactions of the several interconnected neuromusculoskeletal elements of the model and their respective dynamics. Further discussion on the issue of continuous and intermittent control mechanisms operating during postural control shall be presented in a separate section below.
It is worth mentioning that a key parameter for stabilizing the inverted pendulum model was the constant level of the fusimotor activity that adjusted the sensitivity of muscle spindles for each simulation. Without a proper value for both static and dynamic fusimotor activities the pendulum fell at the beginning of the simulation. As previously mentioned, some parameters were randomly distributed in each simulation run [41], hence producing a different set of initial conditions in different runs. This explains why the mean values of fusimotor activities vary (slightly) across the different simulation trials (see Methods). In the context of human postural control, the dependence of an effective control upon the fusimotor activity suggests that the CNS must properly set the muscle spindle sensitivity for the performance of the task [27,44,45]. Another point that should be stressed is that the NMS model was not stabilised without the neuronal activity, i.e., the pendulum fell when the SLC was turned off. This is compatible with the current view that the viscoelastic properties of the muscles around the ankle joint are not sufficient to control upright standing [8,[46][47][48]. In the model, the intrinsic passive ankle joint stiffness was about 70% of the critical toppling torque (see Methods for details), which is in accordance with the experimental estimates at ankle joint rotation angles similar to those obtained in our model (about 0.60 deg on average) [46,48]. The adoption of a constant passive stiffness is a typical simplification (see for instance [5,11]) that was also adopted in the present study. In an analytical study [49], the conclusion was that the feedback provided by the muscle spindles and GTOs is not sufficient to stabilize an inverted pendulum representing the human body. However, the authors did not considered any passive mechanism at the joint level and, hence, the total torque was generated by the active neural controller. Here, as the passive properties were included, the demand of the CNS was about 30% of the necessary stabilizing torque, which is compatible with the experimental estimates [11,46].

Intermittent Recruitment of the Motor Units
The most remarkable result obtained from the proposed NMS model is shown in Figure 4. The intermittent recruitment of MG MUs is a phenomenon recently observed in human experiments [14,16] and the postural control model reproduced this behaviour with a high degree of fidelity. The experimental study in [14] reported that MG MUs were intermittently recruited with a modal frequency of *2 Hz (pooled data from 7 subjects), which is similar to the value observed in [13] for actions of a human subject manually controlling an inverted pendulum. A central hypothesis raised by several recent studies postulates the involvement of an intrinsic predictive mechanism used by the CNS in the performance of postural control [13,15,50]. This intrinsic mechanism, sometimes named as a ''refractory response planner'' [51] and involving a ''psychological refractory period'' [15,50,51], would be responsible for the intermittent actions of the neural controller during the equilibrium maintenance of an unstable load. The simulation results showed that even in the absence of any predictive mechanism (or an internal time setting neural circuit), actions of the neuronal controller occurred at a mean (modal) rate of *2 (4) Hz, i.e., MUs from the MG muscle were recruited with a mean (modal) interval of *500 (250) ms, irrespective of model structure (i.e. Model 1 and Model 2). The models adopted in this study are interpreted as representing a single subject instead of a population of subjects, and the MU intermittence rates observed in the simulations are within the experimental range (2-4 Hz) reported elsewhere (e.g., [11,13,14]). These data suggests that the interplay between a SLC and the muscles involved in the task being performed is sufficient to provide a mechanism underlying the intermittent actions of the CNS during postural control. No complex central mechanism (e.g., predictive, response planner) was needed in our model for the genesis of this control pattern.
As discussed in [14,16], the MG muscle seems to be mostly involved in balance control during standing, while the SO muscle provides a basal torque due to its mostly continuous activity (see Figure 5 and Figure 1C). Regarding the LG muscle, a recent finding showed that this muscle has a minimal or absent activation during the postural control task [16]. The simulation data are in agreement with these experimental results, and suggest that the differences in the organisation of the MG and SO motor nuclei might be responsible for their different actions during postural control. Additionally, the LG muscle was minimally activated during the simulations (the reason why we performed quantitative analysis only on MG MUs), although the recruited LG MUs (less than 30 per simulation) followed a pattern similar to the MG MUs. This different behaviour between the lateral and medial parts of the Gastrocnemius might be explained by a different number of MNs innervating each muscle. The LG muscle has approximately 60% less MUs than the MG muscle (see Methods), and hence for a similar effective synaptic current the number of recruited MUs with similar intrinsic characteristics would be lower for the LG muscle in comparison to the MG. In the model, the SO muscle is more homogeneous, having a high number of low-threshold and smaller twitch amplitude MUs, while MG and LG muscles have an equal proportion of low-and high-threshold MUs generating higher twitch amplitudes (and mostly briefer twitches). The lower proportion of low-threshold MUs (with a similar recruitment range -see [16]) might naturally produce the intermittent recruitment due to fluctuations in sensory feedback. The rationale is that for a similar lower sensory inflow (or effective synaptic current), a low number of MG/LG MUs are recruited, hence producing a low torque at the ankle joint. Conversely, a higher number of SO MUs are recruited producing a basal torque sufficient to counterbalance the static toppling torque. During a forward sway, any small modulation of sensory inflow is sufficient to recruit additional higher-threshold MG/LG MUs, counterbalancing the postural perturbation. As the inverted pendulum returns to a backward position, sensory inflow decreases and the recruited MUs are derecruited (see Figure 4). For the SO, these oscillations in sensory inflow seem to be lower during postural sway (see Figure 5B). However, the analysis is quite complex since the system is operating in closed loop, so that any argument based on causality may lead to logical difficulties. In spite of this limitation, the simulation results indicate that upright standing could be controlled by means of proprioceptive sensory information feedback and a mixture of continuous (SO muscle) and intermittent (mostly MG and LG muscles) action of the CNS. Therefore, the second working hypothesis raised in the present study (that efferent actions of the CNS are continuous during postural control -see Introduction) turned out to be half true, i.e., the leg muscles are activated by a combination of both continuous and intermittent processes.
The results in Figure 6 showed that MG MUs were preferentially recruited when the body leaned forward (panel A in Figure 6) and decelerated (panel B in Figure 6), i.e., MG MUs were mostly recruited in order to counteract the toppling torque due to gravity, pushing the body to a backward position. These simulation results are similar to those obtained from human subjects, as reported in [14]. The authors of the referred study proposed that a strategy of MU recruitment instead of MU rate modulation during upright standing would be generated by the CNS due to the postural task demands. The data from the simulations in Figure 6 reinforce this view. The preference for recruitment coding would be due to the same mechanisms discussed in the previous paragraph, i.e., structural features of the MG motor nucleus and modulation of sensory information due to perturbation from a mean equilibrium position. On the other hand, recent experimental and computer simulation studies have shown that during isometric contractions, the TS torque control relies mainly on rate coding [52] and the variability observed in both torque and EMGs is highly dependent on the MU discharge rate variability. Therefore, the same muscle group (i.e., the TS) is probably being driven according to two different laws depending on the motor task: rate coding for isometric torque control in a very stable condition, and recruitment coding (for the MG/LG muscles) in a more challenging condition, such as erect posture. Interestingly, recent experimental data relating postural sway with isometric torque variability (at similar mean torque values) in young subjects found that they have a positive correlation [53] albeit the first is much larger in magnitude than the latter. As the isometric torque control (seated subjects) involved almost certainly only continuous feedback (mostly from the SO) this experimental result gives support to the dual control mode (continuous and intermittent) that was found in the present simulations for standing posture control.

The Role of Reciprocal Inhibition
Two model structures were simulated in order to investigate whether postural control may be influenced by the reciprocal inhibition pathway (see Methods for details). Recent studies have discussed the importance of reciprocal inhibition in movement control. For instance, [29,30] hypothesised that this inhibitory pathway may be a better source of feedback control since TA proprioceptive activity is unmodulated by the homonymous muscle activation during postural sway.
An interesting result was that in comparison to the model without reciprocal inhibition (Model 1) the complete model (Model 2) showed an increased number of intervals in which SO muscle fibre length was negatively correlated with the COM displacement (see Figure 8). This ''paradoxical'' behaviour was reported in some experimental studies [29][30][31] and was used as evidence to postulate the significant role of reciprocal inhibition in the control of upright quiet standing [29]. The simulation results corroborate the hypothesis that the ''paradoxical'' behaviour of muscle fibre lengths may be generated by the reciprocal inhibition pathway. Nevertheless, no interval with negative correlation was found between MG and LG muscle fibre lengths and COM. In [29], the authors reported that two out of eight subjects showed a larger number of positively correlated intervals for the MG muscle, and they discussed that these subjects oscillated in a more forward position. For the physiologically-constrained set of parameters adopted in the present model the mean equilibrium position was *5 deg forward, which is similar to experimental findings [36]. Therefore, further studies are necessary to better understand the real significance of ''paradoxical'' muscle fibre behaviour and how it emerges during upright stance control. Yet, it is interesting that a highly complex and unexpected biological phenomenon may be partly explained/reproduced by a biologically plausible NMS model, and, therefore, providing neurophysiological clues to its genesis.
Regarding basic postural sway metrics (e.g., COP RMS, MV, and spectral contents) the simulation results did not show large differences between the two model structures (see Table 1), suggesting that reciprocal inhibition is not a fundamental mechanism for postural control.
In spite of the suggestion that TA muscle spindles must be a better (''cleaner'') source of ankle angle feedback than TS muscle spindles [29] the simulation results from Model 1 (without reciprocal inhibition) showed that even ''noisy'' sensory feedback from the TS muscle receptors is sufficient for an adequate postural control. The TS spindle feedback is ''noisy'' in the sense that the TS muscle receptors are signalling a mixture of information from ankle angle changes as well as changes in muscle length and tension due to the MN pool activation.

Model Limitations and Future Research
One conclusion that can be reached from the present simulation results is that mechanisms beyond those included in the model are not strictly necessary to reproduce experimental data from other studies. However, it is not possible to exclude that, despite theoretically not necessary, such mechanisms play a role in human postural control. Specifically, contributions from additional sensory modalities, such as foot soles, joint and skin receptors, vision, and vestibular system, certainly contribute by varying degrees to postural control depending on the particular experimental conditions [2,3,42,54]. Additionally, one cannot rule out the involvement of supraspinal centres (e.g., brainstem, basal ganglia, primary motor cortex) [51,55], specially if the maintenance of upright standing is being learned, such as in infants and adults recovering from a serious medical/neurological disease. Modulations of fusimotor [44,56] and presynaptic inhibition activities [57,58] are examples of important spinal-related mechanisms that certainly play relevant roles too.
In a general context, proprioceptive information from the legs is provided by muscle, joint and cutaneous receptors [54,59]. However, in the NMS model presented here the proprioceptive information was provided exclusively by muscle receptors (i.e. muscle spindles and GTOs), which are postulated as being primary sources of sensory information in response to limb movement [27,28]. The degree of influence of joint and cutaneous receptors on postural control is a controversial issue in the literature. Some experimental findings showed little change in postural sway after ischemia or anaesthesia [60,61], while others showed that stimulation of cutaneous afferents evoked postural changes during quiet standing [54]. The simulation results presented here are in accordance with the former view, but further theoretical/computational and/or experimental studies are required to investigate what is the relative contribution of additional neuronal structures for the maintenance of human upright standing.
Regarding the biomechanics of the human body, it is well known that the inverted pendulum is an approximation for the human body during quiet standing [4,26]. Other expansions, for instance, multi-link and/or multi-dimensional (e.g., including medial-lateral oscillations) models and analyses [24] are very interesting research avenues. However, any biomechanical expansion in large-scale models such as that used in the present study should envisage an increasing number of neuronal and musculoskeletal elements, along with the complexities of their interactions.
Finally, it is noteworthy that we simulated postural control during relatively brief periods (about 30 s). Prolonged unconstrained standing is associated with large changes in body equilibrium position along time [62]. A NMS model to provide approximate postural control during prolonged standing would probably require a reasonably higher complexity than the one employed in the present model and this is certainly an interesting challenge for future research.

Concluding Remarks
Large-scale modelling has been adopted in several fields of modern neuroscience research (e.g., [32,33,[63][64][65][66][67]). To our knowledge this is the first study aimed at modelling the NMS system involved in the control of human upright standing from a more neurophysiological perspective. The main conclusion drawn from the simulation results is that posture control might be, at least in part, mediated by spinal mechanisms, with proprioceptive information being fed back to the spinal neuronal circuitry. Additionally, the results provided evidence that complex phenomena observed in human experiments, such as intermittent actions of the CNS, might not depend on intricate control strategies of complex structures within the CNS. The structure and organisation of the spinal cord, i.e., the types of MUs, synaptic connectivities, the ordered recruitment of MUs, as well as the modulation of proprioceptive information could be sufficient to explain how the CNS controls body position during upright quiet standing in a general sense.

Mathematical Models
The model proposed in this study encompasses four main subsystems (neuronal controller, muscles, proprioceptors, and body biomechanics) that were interconnected to represent the NMS system involved in the control of human upright stance. An overview of this large-scale model is depicted in Figure 9. It is worth mentioning that the model is aimed to study body sway in the sagittal plane during unperturbed stance. In this condition, the posture control task relies mainly on afferent and efferent activities associated with the muscles around the ankle joint (ankle strategy) [9]. Figure 9A shows a schematic view of the neuronal circuitry composed of groups of spinal MNs and INs (see mathematical description below), referred to as the SLC (Spinal-Like Controller). MNs were assembled in motor nuclei associated to the TS (i.e., SO, MG, and LG) and TA muscles, which is the most relevant antagonist group of muscles involved in postural control during ankle strategy [9,29,36]. Stochastic descending commands represented part of the synaptic inputs from the brain that drive the spinal MNs during the maintenance of upright standing. Musculotendon units (MTUs) were represented by Hill-type models (see mathematical description below), which were driven by the spike trains from the spinal MNs ( Figure 9B). Muscle spindles were placed in parallel with the muscle fibres and received commands from Gamma motor neurons (c-MNs), while GTOs were placed in series with the tendon. Proprioceptive feedback was provided by Ia, II and Ib axons mediating: i) monosynaptic Ia excitation; ii) di-synaptic Ib inhibition; iii) di-synaptic II excitation; and iv) reciprocal inhibition from antagonistic Ia afferents. These are fundamental pathways frequently associated with different motor tasks, including upright standing [58]. An inverted pendulum was adopted to represent the body biomechanics ( Figure 9C), which is a first approximation for unperturbed quiet standing [4,5,10,11,20,25,26]. In the following sections, the mathematical details concerning each of these models will be provided.
Spinal neuron models. The MN pool model has been extensively described elsewhere [41,52,68,69]. Briefly, ach typespecified MN (i.e., S-, FR-, and FF-type) was modelled as a twocompartment conductance-based neuron model. Geometrical and electrotonical parameters followed those reported in the literature for lumbar MNs of anaesthetised cats. The somatic compartment included ionic conductances responsible for the genesis of action potentials (Na z and fast K z ) and afterhyperpolarization (slow K z ). Voltage-gated ionic channels were not included in the dendritic compartment (passive-dendrite model) because they are mainly involved in the generation of persistent inward currents (PICs) and the bistable behaviour of MNs [68,69]. It is out of the scope of this study to evaluate the possible role of PICs on human postural control. Motor axons were represented as simple spike conductors, i.e., if a MN generates a spike it is transmitted to the motor end-plate with a given delay depending on the axon conduction velocity (S-type: 44-51 m/s; FR-type: 51-52 m/s; FFtype: 52-53 m/s) and the distance between the spinal cord and the muscle (0.80 m). IN models were represented as point neurons (single-compartment) with the same ionic channels adopted in the MN models [41]. IN passive properties were based on the estimation reported in [70]. For simplicity, all IN models were supposed to have the same dynamical behaviour irrespective of the group they belong (i.e., Ia, II and Ib INs). However, spike thresholds varied linearly along the IN pool (from 10 mV to 20 mV) in order to represent the scattered recruitment of these cells within a given group. The time course of each conductance included in both MN and IN models was simplified [71] in order to speed up the simulation of thousands of neuronal elements. The models have been previously validated with respect to their biophysical properties, for instance, frequency-to-current (f-I) curve, afterhyperpolarization time course, and response to synaptic inputs [41,52,69]. Synaptic conductances followed the kinetic model proposed by [72] with parameters adjusted so that amplitude and duration of excitatory and inhibitory post-synaptic potentials matched the values reported in animal experiments (see [68,69] for details).
As mentioned earlier, MNs were divided in four motor nuclei depending on the muscle they command. Table 2 shows the number and types of MNs adopted for each motor nucleus. These values and proportions follow estimates from the literature [73,74] Since there is a lack of information regarding the exact number of spinal INs mediating each pathway represented in the present model, a fixed number of 350 INs was adopted for each group.
Model of the musculotendon units. MTUs produce both muscle torque and the surface electromyogram (EMG) (see Figure 9A). The phenomenological model of muscle EMG was extensively reported and validated in other studies [41,52,68]. Biomechanical properties of the muscles were represented by Hill-type models [75][76][77] containing elements as depicted in Figure 9B. A non-linear in-series elastic element (c T ) accounted for the passive property of tendon and distal aponeurosis. Two passive parallel elements (K PE and b) were adopted so as to represent the viscoelastic properties of the muscle fibres. In addition, the contractile properties of type-I (slow) and type-II (fast) muscle fibres were represented by two contractile elements (CE I and CE II ), which produce active force in response to MN activation. The pinnation angle between the muscle fibres and tendon (a) was also represented in the model. A mass (m) was introduced to provide computational stability. The Supplemental Material shows the frequency response of the SO muscle model that was estimated according to the experimental procedures reported in [78]. The results presented in this Supplemental are useful as a complementary validation for the modelling described in what follows.
Activation dynamics. Activation signals to the muscles were obtained by a filtering process (second-order critically-damped MNs from the TS motor nuclei receive constant intensity descending commands during the maintenance of upright stance. Proprioceptive feedback is provided by Ia, II and Ib afferents from muscle spindles and Golgi tendon organs (GTOs). The information carried through these afferents is fed back to the spinal cord representing fundamental pathways (or neuronal circuits) involved in the postural control task, such as Ia monosynaptic excitation, Ib di-synaptic inhibition, group II di-synaptic excitation, and reciprocal inhibition from antagonistic Ia afferents. Activity from Gamma motor neurons (c-MNs) set the sensitivities of the muscle spindles. Ankle joint torque (T A ) that drives the body biomechanics (to compensate for the gravitational toppling torque) is given by the sum of the torques produced by the muscles (T m ) and the passive ankle joint torque, represented by the passive ankle joint stiffness (K A ) and viscosity (B A ). The body angle (h A ) is the resultant output of the inverted pendulum acted on by gravity and by T A . It indirectly (by means of Equations 8 and 9) defines moment arms (m x ) and musculotendon (MTU) lengths (L MT ), which are used to define the dynamics of both MTUs and muscle receptors (see dashed lines). Additionally, h A is fed back without delay to account for the intrinsic passive joint impedance (stiffness and viscosity). (B) Hill-type model used to represent the viscoelastic and contractile properties of the MTUs. Muscle fibres are represented by parallel passive elements (muscle fibre stiffness, K PE , and viscosity, b) and two contractile elements (CE) representing the contractile properties of type-I and type-II muscle fibres. A pinnation angle (a) is adopted to represent the angle between the muscle fibres and the aponeurosis. In addition, a mass (m) is used to increase the computational stability. Passive properties of tendon and distal aponeurosis are represented by a lumped non-linear stiffness (c T ). Muscle spindle is placed parallel to the muscle fibres, while the GTO is in series with the tendon. (C) Inverted pendulum model used to represent the body biomechanics during the upright quiet stance. Arrows indicate the positive direction of each variable (see the description of each variable in the text). doi:10.1371/journal.pcbi.1003944.g009 linear system) applied on the MN spike trains. The filter output was followed by a non-linear function that provided a smooth saturation [52,79,80]. The non-linear function (Equation (1)) accounted for the mechanisms responsible for the muscle force saturation (e.g. Ca zz release saturation in the sarcoplasmic reticulum): in which, a MU SAT (t) is the saturated activation signal produced by a given MU; a MU (t) is the activation signal of a given MU before the saturation (i.e., after the second-order linear filter); and c is the shape parameter of the saturation function.
For each MU the parameter c was adjusted so that the tetanic activation was achieved at a given firing rate (see [52]), i.e., a low (high) threshold MU achieved its tetanic activation at a low (high) firing rate. There was also an amplitude scaling depending on the MU type, i.e., S-type MU produced activation signals with amplitudes relatively lower than F-type MU. For more details regarding the distribution of parameter values, see [52].
The sum of all activation signals produced by S-type MUs resulted in the activation signal to the CE I [a I (t)]. Similarly, the activation signal to the CE II [a II (t)] was generated by the sum of all activation signals produced by FR-and FF-type MUs. These global activation signals were normalised with respect to the maximum muscle activation, which was calculated as the sum of all tetanic activations produced by the MUs of a given muscle. A similar approach was adopted in [75].
Passive properties. Force-length relationships for the elastic elements were adopted from [76,77], while the viscous force was simply proportional to the muscle fibre stretch velocity. The total parallel passive force (Equation (2)) was normalised by the maximum muscle force (F 0 ) estimated from the physiological cross-sectional area of each muscle [75]. Muscle fibre length was normalised by the optimal fibre length (L 0 ). For the tendon, the adopted passive force-length relationship (Equation (3)) was the one reported by [75]. Similarly, the in-series passive force was normalised by F 0 and the tendon length was normalised by L T 0 , which is the value at which the force produced by the tendon in which, K PE and b are the elastic (in F 0 =L 0 ) and viscous (in F 0 :s=L 0 ) parameters, respectively;L L M is the normalised muscle fibre length (in L 0 ); E M 0 is the normalised muscle strain; andṼ V is the normalised muscle fibre velocity (in L 0 =s).
in which, k T is a constant that defines the curvature of the toe region in the force-length relationship; c T is the tendon stiffness (in F 0 =L T 0 ); L T is the normalised tendon length (in L T 0 ); and L T r is the length at the onset of the linear region.
Muscle fibre length and tendon length are related by Equation (4). The pinnation angle varied with muscle fibre length (Equation (5)) so as to maintain the muscle thickness constant [81].
in which, L MT is the MTU length (in m), whose value is given by inverse kinematics data (see below); and a(L L M ) is the pinnation angle as a function of the normalised muscle fibre length.
in which, a 0 is the initial pinnation angle (in rad). Values for the passive parameters of each muscle are presented in Table 3. These values were based on several data reported in the literature [75][76][77]82].
Contraction dynamics. The force generated by contractile elements is given by Equation (6). Force-length (FL) and forcevelocity (FV ) relations followed the proposal by [75] for type-I and type-II muscle fibres (see their Table 1 for details).
Muscle fibre velocity was calculated by integrating muscle fibre acceleration (Equation (7)), which was obtained by applying the Newton laws to the mechanical system described above (see Figure 9B for a schematic view). It is worth noting that positive velocity means muscle fibre stretch. Accordingly, muscle fibre length was calculated by integrating muscle fibre velocity.
Musculoskeletal geometry. MTU lengths (L MT ) and moment arms with respect to the ankle joint (m x ) were calculated according to the approach proposed by [83]. For each muscle, fourth-order polynomials were used to adjust the relations between MTU length and ankle joint (Equation (8)), as well as moment arm (m x ) and ankle joint angle (Equation (9)). The coefficients of these functions are presented in Table 4. Relations for the MG and LG (biarticular muscles) were obtained for the knee angle at zero degree (fully extended). Inverse kinematic data were reported in [82].
Models of proprioceptors and afferent fibres. A muscle spindle model previously reported by [84] was placed parallel to the muscle fibres (see Figure 9B). It represents the dynamics of bag 1, bag 2 and chain intrafusal muscle fibres. Ia mean firing rate is a non-linear function of all three instrafusal fibre outputs, hence conveying both static and dynamic information about muscle stretch. On the other hand, type II mean firing rate is produced by a contribution of bag 2 and chain intrafusal fibre outputs, which are dependent on the static component of muscle fibre stretch [84]. Normalised muscle fibre length (L L M ) is the input to the muscle spindle model. In addition, bag 1 intrafusal fibre is sensitive to dynamic fusimotor activity (c d ), while bag 2 and chain intrafusal fibres are sensitive to static fusimotor activity (c s ). In this study, biophysical properties of c-MNs were not represented. Static and dynamic fusimotor activities were both modelled as Gaussian random processes with variance equal to 3% of the mean value (empirically chosen). The mean value was adjusted in order to stabilize the biomechanical system (see Simulation Protocol below). Almost all model parameters were maintained equal to those reported by [84] (see their Table 1); however, the gain of each intrafusal fibre was adjusted (G bag1~7 000, G bag2~3 800, and G chain~3 000) in order to produce Ia and type II firing rates compatible with experimental data from humans [59].
A GTO model reported by [85] was in-series with the tendon (see Figure 9B). It represents the lumped dynamics of a population of Ib afferents in response to the force produced by the tendon (F T ). Basically, this model encompasses a static non-linear function (Equation (10)) in series with a linear dynamics (Equation (11)) resulting in the afferent firing rate [f Ib (t)]. The transfer function that represents the linear dynamics was transformed to a digital filter (bilinear transformation) so as to reduce the computational complexity of the system. Individual Ib afferent activity was obtained by dividing the GTO model output by the total number of Ib afferents in each muscle. Model parameters were chosen to produce individual Ib firing rates compatible with those reported in the literature [59].
in which, G 1 and G 2 are adjustable parameters, chosen equal to 60 Hz and 4 N, respectively.
These values were based on data from several studies [75][76][77]82]  The outputs of the modelled muscle receptors are continuous functions of time that represent the instantaneous firing rate of Ia, II, and Ib afferents. These continuous functions [Af MR (t)] (with MR being either, Ia, II, or Ib muscle receptor afferent) were converted into independent stochastic point processes to represent the spiking activity travelling along a bundle of afferent fibres with a given conduction velocity. The i-th spike train was modelled as a non-homogeneous Gamma point process with a mean intensity [Af i (t)] modulated by the output of a muscle receptor (spindle or GTO) and with an order (or shape factor) adjusted to produce a variability compatible with experimental data [86]. The total number of afferents for each muscle (see Table 5) and conduction velocities (see Table 6) were based on estimates from the literature [58,87,88]. Afferent fibres were recruited according to a recruitment law given by Equation 12 [79].
in which, Af i (t) is the mean firing rate (intensity) of the i-th afferent; Af MR (t) is the output of the muscle receptor models (spindle or GTO) for a given afferent type; RT i and IFR i are the recruitment threshold and the initial firing rate of the i-th afferent, respectively. The i-th afferent was recruited when the output produced by a specific muscle receptor afferent [Af MR (t)] crossed a recruitment threshold (RT i ). Irrespective of the afferent type (i.e., Ia, II, and Ib) the RT i was linearly varied along the afferents from 0-50 Hz, i.e., when the specific muscle receptor output reached 50 Hz all afferent fibres of a given type were recruited. At the time of recruitment the i-th afferent fibre discharged with an initial firing rate (IFR i ) that was randomly varied along the bundle of afferents following a Gaussian random number (mean = 5 Hz; standarddeviation = 2.50 Hz). The values of RT i and IFR i were empirically chosen so as to produce the recruitment behaviour of afferent fibres observed in human experiments [59,89]. Table 6 shows the connectivity ratios between sensory afferents and the spinal neurons (MNs and INs), with values based on data reported in [58]. Maximum synaptic conductances (see Table 6) were empirically adjusted in order to produce a stable and meaningful model response.
Body biomechanics. As mentioned above, body biomechanics during the upright quiet stance was represented by an inverted pendulum ( Figure 9C). Equations of motion of body COM with respect to the ankle joint (fixed pivot point) are given by Equations (13) and (14). It is noteworthy that our analysis was restricted to the sagittal plane, which represents the anteroposterior displacement of the body COM. Anthropometric parameters were empirically chosen to represent a young adult with m b = 60 kg and h b = 0.85 m.
in which, J b is the body inertia (in kg.m 2 ) given by J b~4 =3:m b :h 2 b ; m b is the body mass excluding the feet (in kg); h b is the COM height with respect to the ankle joint (in m); T A (t) is the restoring torque around the ankle joint (see Equation (13)); and g is the gravity (set to 9.81 m/s 2 ).
in which, T m (t) is torque produced by the leg muscles (in Nm); B A and K A are the viscoelastic coefficients that represent the passive property of the joint (in Nm.s/rad and Nm/rad, respectively). Muscular torque (T m (t)) is given by the net torque produced by the TS and the TA muscles multiplied by two (for simplicity, we have assumed that both legs produce the same torque and have the same control mechanism). The torque produced by the TA muscle was defined as positive, so that activation of the TS tends to restore the body from the toppling torque produced by gravity. The viscoelastic passive torque was introduced in the model because Hill-type muscle models cannot properly produce muscle stiffness [47]. A preliminary simulation was performed in order to estimate the ankle joint stiffness due to the passive properties of the modelled Hill-type MTUs. In this simulation the neuronal controller was turned off and the ankle joint angle [h A (t)] was randomly varied (Gaussian stochastic process). The resultant (passive) muscle torque [T m (t)] was measured at the output of the model. A spectral analysis was carried out to estimate the resultant passive ankle joint impedance (see [90,91]), and the passive ankle joint stiffness due to the Hill-type MTUs was estimated as the mean impedance value in the range 0-1 Hz (static component). The passive ankle joint stiffness produced by the Hill-type MTUs was estimated to be only about 5% of the critical stiffness (m b :g:h b ). An additional passive ankle stiffness (K A ) was adopted equal to 65% m b :g:h b , resulting in a total passive ankle joint stiffness of about 70% [46,48]. A small viscosity (B A equal to 5.81 Nm.s/rad [46]) was also included in the model. These viscoelastic intrinsic properties of the ankle joint accounted for part of the postural control through a mechanism frequently referred to as preflexes [92] (see Figure 9A).
COM and COP displacements. Typically in human experiments, both COM (y G ) and COP (y P ) displacements (projected on the anteroposterior axis in the present study) have been recorded to evaluate the body sway during upright stance. Here, these biomechanical variables were calculated according to Equations (15) and (16).

Simulation Protocol
The models described above were written in Java programming language (Oracle Corp) and simulated in an open-source webbased application developed in the Eclipse IDE (The Eclipse Foundation). Model source code is available for download at the website http://remoto.leb.usp.br/ and at a public repository (http://dx.doi.org/10.6084/m9.figshare.1029085). Differential equations were numerically solved by a fourth-order Runge-Kutta method with a 50 ms time step.
Recent experimental studies argued that the reciprocal inhibition from antagonistic Ia afferents (from TA muscle spindles to TS MNs -see Figure 9 A) has an important role in the postural control due to the ''orthodox'' modulation of TA length and the lack of TA contraction during standing [16,29]. To investigate this potential role of reciprocal inhibition pathway, two different model structures were used in this study. The first one, hereafter referred to as Model 1, excluded the reciprocal inhibition pathway so that all proprioceptive information was provided exclusively by the TS sensory afferents (autogenic pathway). Conversely, Model 2 is the complete model fully described above. All common parts in both models have the same parameter values.
Three independent simulations (30-s duration) were carried out with each model. During the first second of simulation, the movement of the inverted pendulum was restricted so that the neuromuscular system reached its steady-state condition. This prevented instabilities due to the lack of neural activity. In order to provide a basal activation, descending commands were modelled as 400 homogeneous stochastic Gamma point processes with an individual mean rate equal to 50 Hz and a shape factor equal to 25 (coefficient of variation equal to 20%). With this descending drive (without any feedback information) the MN pools produced a small basal muscle torque equal to approximately 2% of the maximum torque. As mentioned earlier, the fusimotor drive was modelled as a Gaussian stochastic process with a small variance. The mean values of static and dynamic fusimotor activities were adjusted for each simulation so as to produce outputs in which the inverted pendulum was oscillating around an equilibrium point for the whole simulation duration. For Model 1 the values of static (dynamic) fusimotor activities were 31.10 (33.30), 31.20 (33.30), and 31.50 (33.60). Similarly, for Model 2 the values of static (dynamic) fusimotor activities were set to 32 (33.80), 32 (33), and 31.70 (32.70).

Data Analysis
Several output variables were analysed in the present study: i) Spike trains from MNs, INs and afferent fibres; ii) EMGs from the TS muscles; iii) COM and COP anteroposterior displacements; iv) Muscular torque; v) Muscle fibre lengths; vi) Joint angle and angular velocity. These raw variables are available for download at a public repository (Model 1 -http://dx.doi.org/10.6084/m9. figshare.1027609; Model 2 -http://dx.doi.org/10.6084/m9. figshare.1029084). Analyses were performed on 22.50-s duration time series, so that the initial 5 s and final 2.50 s were removed due to transient responses and data filtering (see below).
COM, COP, joint angle, angular velocity, muscle fibre length, and torque time series were detrended and downsampled to 2 kHz. In order to compare model behaviour with experimental data recorded with force plates, time-and frequency-domain metrics were estimated from the COP, for instance, RMS, MV and 50% power frequency [42]. Power spectra were estimated by Welch's method. EMGs were downsampled to 2 kHz, rectified, and lowpass filtered at 2 Hz (zero-phase fourth-order Butterworth filter). Table 6. Sensory feedback parameters adopted in the model.

Parameter
Value Unit Similar to the studies reported in [9,19,43], cross-correlation analysis was performed between COM and COP time series, as well as COP and EMGs. In addition, similar to the analysis reported in [29], correlation coefficients were computed between COM and muscle fibre lengths. For the latter analysis COM and muscle fibre length time series were binned to 3-s duration intervals (without overlap) so that the correlation coefficients were calculated for each bin (see Figure 7).
In order to quantify the degree of intermittency of a given muscle, the activation ratio of individual MUs was measured. This ratio was calculated as the percentage of time a given MU had interspike intervals lower than 250 ms [16]. An activation ratio equal to 0 indicates that the MU was inactive during the whole simulation time, whereas a ratio equal to 1 indicates a continuous activity. Therefore, lower activation ratios represent an intermittently discharging MU. Additionally, phase plots were constructed to evaluate the MU recruitment during quiet standing. Recruitment phase plots were constructed by binning joint angle, angular velocity, muscular torque and spike trains in 100-ms nonoverlapping windows. Within each window, the number of recruited MUs was counted and a circle with diameter propor-tional to the counting number was plotted in 2-dimension anglevelocity (see Figure 6A) or angle-torque graphs (see Figure 6B). Data from the three simulations were pooled together in the same phase plot. A similar analysis was performed in [14].
When applicable, statistical analysis was performed with a significance level set at 0.05. All analyses were performed in Matlab (The MathWorks Inc) and SPSS Statistics (IBM Corp).

Supporting Information
Supporting Information S1 Frequency response of muscle model. This supporting information file provides additional validation data for the muscle model. (PDF)

Author Contributions
Conceived and designed the experiments: LAE AFK. Performed the experiments: LAE RNW. Analyzed the data: LAE RNW AFK. Contributed reagents/materials/analysis tools: LAE RNW. Wrote the paper: LAE RNW AFK.