From the motor cortex to the movement and back again

The motor cortex controls motor behaviors by generating movement-specific signals and transmitting them through spinal cord circuits and motoneurons to the muscles. Precise and well-coordinated muscle activation patterns are necessary for accurate movement execution. Therefore, the activity of cortical neurons should correlate with movement parameters. To investigate the specifics of such correlations among activities of the motor cortex, spinal cord network and muscles, we developed a model for neural control of goal-directed reaching movements that simulates the entire pathway from the motor cortex through spinal cord circuits to the muscles controlling arm movements. In this model, the arm consists of two joints (shoulder and elbow), whose movements are actuated by six muscles (4 single-joint and 2 double-joint flexors and extensors). The muscles provide afferent feedback to the spinal cord circuits. Cortical neurons are defined as cortical "controllers" that solve an inverse problem based on a proposed straight-line trajectory to a target position and a predefined bell-shaped velocity profile. Thus, the controller generates a motor program that produces a task-specific activation of low-level spinal circuits that in turn induce the muscle activation realizing the intended reaching movement. Using the model, we describe the mechanisms of correlation between cortical and motoneuronal activities and movement direction and other movement parameters. We show that the directional modulation of neuronal activity in the motor cortex and the spinal cord may result from direction-specific dynamics of muscle lengths. Our model suggests that directional modulation first emerges at the level of muscle forces, augments at the motoneuron level, and further increases at the level of the motor cortex due to the dependence of frictional forces in the joints, contractility of the muscles and afferent feedback on muscle lengths and/or velocities.


Introduction
Even simple arm movements such as reaching require complex interactions among the central and peripheral nervous systems, and skeletal muscles to generate the intended arm movements. Over three decades, a lot of effort has been put into understanding neural mechanisms a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 controlling reaching movements [1][2][3][4]. Reaching is broadly defined as the arm's movement starting at some initial position in space and ending at a target position. In experiments, unperturbed reaching movement usually occurs along a straight-line trajectory with a bellshaped velocity profile [5]. Dynamically, reaching movements result from complex concurrent or sequential activation patterns of multiple muscles used to accelerate and then, slow down and stop the arm along the intended trajectory. To generate the required muscle activation patterns, the motor cortex needs to solve a corresponding "inverse problem" and, based on this solution, provide the appropriate dynamical inputs to the spinal circuits [6][7][8].
In reaching tasks, the relationships between neuronal activity in the motor cortex and movement parameters have been widely debated, and remain controversial [9][10][11][12]. It has been suggested that the neuronal activity in the primary motor cortex (M1) encodes such movement parameters as direction [2,[12][13][14], hand position [15][16][17][18], velocity [19], acceleration [20], and reaching distance [13,18]. However, other studies have argued that neural activity in the motor cortex correlates with kinetic variables, such as forces and torques [21,22]. In 1982, Georgopoulos et al. demonstrated for the first time a correlation between neuronal activity in the motor cortex and the direction of reaching movement [2]. They showed that the average firing rate of M1 neurons during reaching movements varied with the direction of movement, and that each M1 neuron had a preferred direction (PD) for which its average firing rate was maximal. Since then, directional tuning has been ubiquitously considered as a key property of neural activity in the motor cortex. However, it remains controversial whether directional preference is the fundamental property of cortical neurons or a side effect of muscle activity or other movement features [11,[23][24][25][26][27].
Although many movement parameters correlate with cortical activity, the cause of these correlations is still not clear. For example, it has been suggested that the directional sensitivity of cortical neurons is the result of a specific organization of inhibitory interactions within and between neuronal columns in the motor cortex [28,29]. A competing viewpoint is that cortical activity is related to the activity of corresponding muscles that have anisotropic properties and thus, form cortical directional tuning [27,[30][31][32]. Moreover, the contribution of the spinal cord circuitry to directional modulation is not well understood.
Mathematical models have been used to better understand the relationship between the activity of neurons in the motor cortex and movement parameters during reaching [32][33][34]. However, previous models did not consider the spinal cord network and/or length/velocitydependence of contractility of the muscle controlling the movement. In the present study, we have developed an integrative mathematical model of a motor control system that incorporates cortical neuronal populations, complex spinal neural circuits controlling arm muscles and receiving afferent feedback from them, and a two-joint arm actuated by these muscles that performs reaching movements in a 2D space. The arm model includes the shoulder and elbow joints whose movements are generated by pairs of flexor and extensor muscles controlling a single joint and two bi-articular flexor and extensor muscles controlling both joints at the same time. The arm model is very similar to the model developed by Lillicrap and Scott [32]. Unlike Lillicrap and Scott who used a complex cortical neuronal network including a learning system to control the arm, we focused on a spinal cord network which receives aggregate signals from six supraspinal neuronal pools. The cortical activity (hereinafter referred to as a motor program) was calculated by solving an inverse problem based on a straight-line trajectory to a target position and a predefined bell-shaped velocity profile. This approach allowed us to generate and analyze a variety of activity profiles of cortical neurons that may be used to perform movements to the same or different target positions. We evaluated our findings in the context of different theoretical hypotheses concerned with the neural control for reaching [30][31][32]. Particularly, Lillicrap and Scott [32] showed that limb geometry, intersegmental dynamics, and the force length and force velocity properties of muscle are the main causes of directional preferences of cortical neurons during reaching movements. In agreement with Lillicrap and Scott, our results show that the muscle length and velocity dependent contractile components of the muscle forces are the root causes of the directional preference of spinal motoneurons. Furthermore, we showed that afferent feedback has significant impact on the directional behavior of cortical neurons and proposed that directional dependence of the mean firing rates of M1 neurons primarily results from afferent feedback signals carrying information muscle lengths and velocities. Moreover, our results show that the spinal cord circuit may play a comparable or even more significant role in directional tuning of M1 neurons than the contractile components of the muscle forces. The model reveals the mechanisms by which directionally indifferent torques in the arm joints imply directionally tuned cortical activity.

Results
The model presented here comprises three main modules: Arm, Spinal Cord, and Motor Cortex ( Fig 1A). The Arm module is modeled as a mechanical system of two rigid segments and two joints (shoulder and elbow) controlled by six Hill-type muscles, namely: the shoulder flexor (SF) and extensor (SE), the elbow flexor (EF) and extensor (EE), and the two-joint extensor (BE) and flexor (BF) (see Fig 1C). Arm movements are restricted to a horizontal plane and are produced by coordinated activation of the muscles. The Spinal Cord module has neural circuits that receive descending signals from the Motor Cortex, and relay them through spinal motoneurons to the corresponding muscles. The Spinal Cord neurons receive afferent feedback from the muscles and form local reflex circuits. These include (1) monosynaptic excitation of homonymous motoneurons by Ia muscle afferents; (2) reciprocal inhibition between the antagonistic flexor and extensor motoneurons via Ia interneurons receiving Ia afferents; (3) non-reciprocal inhibition of homonymous and synergistic motoneurons by Ib muscle afferents via the Ib interneurons; and (4) recurrent inhibition of motoneurons via corresponding Renshaw cells (see Fig 1B). The Motor Cortex module projects the activity of six cortical neuronal populations (M1-N) down to six motoneuron pools in the Spinal Cord to control movements of the arm (see Fig 1A and 1B). This descending command quantitatively represents the aggregate cortical input rather than specific mono-or poly-synaptic projections. A detailed description of the model including mathematical formulations is provided in the Methods section.

Directional modulation of cortical neurons
Using our mathematical model, we simulated center-out reaching tasks in 8 different directions and analyzed the activity of each M1 neuron (Fig 2). We found that the firing rate of these neurons depended on the movement direction, such that it was maximal for one direction (preferred direction, PD, denoted by Θ PD ), and minimal for the opposite direction (anti-PD, 180 o apart from Θ PD ). Fig 3 shows an example of the neural activity of a population of M1 cortical neurons that projects to the elbow flexor motoneuron and is therefore responsible for the contraction of the elbow flexor muscle during reaching movements in all 8 directions. The activity profiles (outer traces in Fig 3) show dynamical and quantitative changes when the movement direction changes. This particular population of M1 neurons showed the highest amplitude of activity when reaching movements were made in the 270 o -315 o direction, and the lowest amplitude of activity when reaching movements were made in the 90 o -135 o direction. The firing rate of this particular neuron averaged over the entire reaching movement exhibited the highest value for the 270 o direction and the lowest value for the 90 o direction (polar curve in the center of Fig 3). This type of directionality has been previously observed in The motor cortex sends motor signals to the spinal cord neuronal network which sends its outputs to the muscles. The spinal cord combines motor signals with afferent feedback to generate the motorneuron outputs. (B) Organization of interconnections between Renshaw cells (RC), motoneurons (MN) and other interneurons in the spinal cord network. Motoneurons send their outputs to their corresponding arm muscles. Ia and Ib inputs are the feedback signals from the muscles. (C) A 2-joint arm with six muscles: four major flexor and extensor muscles about the shoulder and elbow joints, and two biarticular muscles controlling both shoulder and elbow joints. SF, EF and BF represent shoulder, elbow and biarticular flexors, and SE, EE and BE represent shoulder, elbow and biarticular extensors, respectively. distribution for the 8 movement directions. All six M1 neurons-corresponding to SF, SE, EF, EE, BF, and BE-showed similar changes in their activity profiles and averages in conjunction with the change in movement direction, regardless of the torque distribution. We computed single-trial averages of the M1 neuron's firing rate over the duration of each movement. For each reaching direction, we computed the mean of the single-trial averages to produce a multitrial average. We found that multi-trial averages varied in an orderly fashion with movement direction (Fig 4, black curves). Multi-trial averages were maximal for one unique movement direction, which was the PD for that particular M1 neuron. Each M1 neuron had a PD distinct from the other five neurons. The multi-trial activity of four M1 neurons (corresponding to SF, EF, EE, and BF) demonstrated coefficients of determination (R 2 , a measure of the precision of the regression fit) greater than 0.7 when fitted to cosine tuning curves (Fig 4, green (Fig 4). We have found that there is no change in preferred directions when the model uses only the shoulder and elbow muscles (by removing the biarticular muscle), which is in an agreement with Lillicrap and Scott [32]. However, the simulations show that biarticular muscles have significant impact on the coefficients of determination of the tuning curves.

Directional modulation of spinal motoneurons and afferent feedback
The analysis of activity of both spinal motoneurons and Ia afferents showed directional tuning properties similar to M1 neurons (see Fig 5). However, the correlation coefficients between  [13]). Moreover, the PDs of Ia afferents were opposite to the PDs of their corresponding M1 neurons ( Fig  6A and 6B).
In the model, directional modulation of the Ia afferent feedback signals emerges from their explicit dependence on muscle lengths and velocities (see Eq (5)). This is consistent with the results of Jones et al. [35] who showed that directional tuning of individual muscle afferents during voluntary wrist movements in humans can be predicted from the length changes of the corresponding muscle.
To quantify the effects of feedback on cortical activity, we increased the amplitude of feedback signals to spinal motoneurons 2-, 5-and 10-fold. The outcome showed that as the amplitude of feedback signals increased, cortical activity became more directionally tuned to a PD opposite to the corresponding feedback signal's PD. For example, when feedback inputs were amplified by 2, 5 and 10-fold, the coefficient of determination for the directional tuning curve of SE-M1 activity increased from 0.47 to 0.69, 0.85, and 0.9, respectively. Moreover, removing the feedback signals decreases the coefficient of determination to 0.1. These results support the idea that the cortical activity may be modulated by afferent feedback signals. Obviously, amplification of the feedback signals does not affect the directional tuning curves of the spinal motoneurons.
To examine the effect of reaching distance on directional modulation of cortical and feedback activity, we tested 6 different reaching distances and observed no significant change in the directional tuning curves of M1 neurons and Ia afferents. The PDs of M1 neurons and Ia afferent feedback remained relatively constant over the 6 reaching distances used (Fig 7), suggesting that PD is independent of reaching distance. Similarly, Fu, Suarez et al. [13] showed that PDs of neurons in the premotor and primary motor cortices of monkeys were  [13] showing that the PDs of 17 recorded cells over varying reaching distances are preserved). (C) The difference between the PD of each population and that of its corresponding Ia feedback is about 180º over six different reaching distances.
https://doi.org/10.1371/journal.pone.0179288.g007 independent of reaching distance (see Fig 3A in [13]). The preservation of PD across different reaching distances and speeds was also observed by Churchland, Santhanam et al. [36]. Moreover, the relationship between PDs of M1 neurons and Ia afferents was also independent of reaching distance; the difference between PDs of M1 neurons and their corresponding Ia afferents remained approximately 180 o over all 6 reaching distances ( Fig 7C).
Directional modulation of the contractile components (force-length and force-velocity) in the muscle force While reaching movements are performed in different directions, six arm muscles in our model contracted (shortened) or stretched (lengthened) with different magnitudes depending on the direction of the movement. Fig 8 shows the directional modulation of the muscle force. The averaged force-length ( Fig 8A) component of each muscle was highest for a distinct (preferred) direction, and lowest for the opposite (anti-preferred) direction, with a smooth transition in between. Cosine functions provided excellent fits (R 2 > 0.95, see Fig 8A, green curves) for the force-length components of all six muscle. The contractile components of the antagonist muscles demonstrated opposite PDs and tuning curves (Fig 8A), consistent with the findings of Cherian A. et al. [37], who showed that EMG activity of biceps and triceps had opposite directional tuning.
Similar to force-length components, the averaged force-velocity component of each muscle strongly correlated with the movement direction in our simulations. Specifically, it exhibited the same PD as that of the averaged force-length component, and had a similar tuning curve ( Fig 8B).

Directional modulation of other movement variables
An identical velocity profile, with fixed peak velocity, was used to model the arm's movements in all directions. Accordingly, the magnitude of the acceleration of the arm endpoint was not dependent on the direction; neither were the magnitudes of the joint torques (which primarily depend on joint angular accelerations). To determine if muscle forces had any correlation with movement direction, we performed 50 simulations where the torque distribution parameter d was varied from trial to trial in the range (0.5~1) using a uniform probability distribution. The results demonstrated that muscle forces were very weakly directionally tuned compared to cortical input tuning, muscle contractile component tuning, and feedback tuning (Fig 9); for all six arm muscles, the correlation coefficient between force and movement direction was relatively small (0.05 to 0.4). Muscle forces had large standard deviations (2.5 to 10 N), and although average muscle forces slightly changed with direction, the change was not significant. We further checked that weak correlations shown in Fig 9 emerge from the effects of joint viscosities (frictional forces, see Methods), as eliminating viscose friction forces in the joints decreased the six muscle forces' correlations to almost zero (not shown). In summary, muscle forces are weakly modulated by the direction of the reaching movement compared to the modulation of spinal motoneurons and cortical neurons.

Shoulder-centered reference frame
Although Cartesian coordinates have been frequently used as a reference frame for center-out reaching tasks, coordinate systems connected to the joints (shoulder or elbow) can also be used to represent the direction of the arm movement. It was previously suggested that cortical neurons, both in the motor and premotor areas, encode movement directions with an intrinsic coordinate system representing the shoulder-centered reference frame [38,39]. Moving the workspace without changing the shoulder location can potentially affect the PDs and tuning curves of cortical neurons because of the changes in arm posture. To investigate this, we shifted the workspace in our model to the left and right, using the shoulder joint as a reference frame ( Fig 10A). The whole workspace (including initial position and 8 targets) was translated with respect to the shoulder joint. The vector from the shoulder joint to the initial position (Fig 10A, black vector) was used as a reference for this transformation. The default workspace, used in most of our simulations, had coordinates (0.0, 0.4) for its initial position, and the vector from the shoulder joint to this initial position was parallel to the positive y-axis (Fig 10A, Center). Fig 10B (Center) shows the PD distribution of M1 neurons for the default workspace. When the workspace was rotated counter-clockwise by 45 o with respect to the shoulder joint (Fig 10A, Left), all preferred directions were shifted in the same direction by approximately the same angle (44.83 o ± 2.04 o ) (Fig 10B, Left). Similarly, when the workspace was rotated clockwise by 45 o with respect to the shoulder joint ( Fig 10B, Right), all PDs were shifted in the same direction by approximately the same angle (45.67 o ± 3.8 o ) (Fig 10B, Right). The PD shift of the EF and EE neurons was exactly 45 o for both clockwise and counter-clockwise workspace transformations. Transforming the workspace by 20 o clockwise or counterclockwise also shifted the PDs by 20 o clockwise or counter-clockwise (not shown). In other words, a shift in the reference frame resulted in a similar shift in the PDs of M1 neurons, which means that the relationship between movement direction and activity of M1 neurons is invariant once the direction is defined relative to the line connecting the shoulder and the initial point of the movement. In summary, rotating the workspace about the shoulder joint shifts all PDs of M1 neurons by the same angle in the same direction. Similarly, Lillicrap and Scott [32] also showed that preferred directions shifted to the left or to the right when the work space was shifted to the left or right, respectively (see Fig 6 in their paper). The changes in muscle lengths during the movement largely depend on the relative direction to the shoulder joint because of the rotational symmetry of the arm geometry. Hence, in our model the directional modulation manifests as dependence on the angle between the direction of movement and the direction from the initial position to the shoulder joint. This is consistent with an idea that the motor cortex encodes direction based on the shoulder reference frame, suggested in other experimental [38,39] and modeling [33] studies.

Directional modulation and the movement distance
To evaluate the relationship between cortical activity and reaching distance, we simulated center-out reaching tasks in 8 directions with 7 reaching distances per direction. Since the reaching time was fixed based on previous experimental studies [13,40], both the peak velocity and peak acceleration increased as the reaching distance increased. In our model, the average activity of each M1 neuron increased monotonically with the reaching distance (Fig 11). However, the changes in average firing rate with respect to the distance were different across M1 neurons, and depended on the movement direction as well. Our simulation results are consistent with experimental findings of Fu et al. [13] who also showed a direction dependent increase in the average firing rate of cortical neuron with increasing reach distance (see Fig 4 in [13]). Due to constant reach time, greater reaching distances result in higher peak acceleration of the arm endpoint. In our model, the dependence of cortical activity on the reaching distance emerges from strong correlation of the latter with peak acceleration. Greater acceleration requires stronger torques at the joints, and hence, stronger activation of M1 neurons. To evaluate this hypothesis experimentally, it would be sufficient to vary the reaching time in proportion to reaching distance, such that longer reaching distances require a greater reaching time. In this case, the peak velocity would be fixed, and contrary to the task presented here, the peak acceleration would decrease as the distance increases. In our model, such a change results in the decline of average cortical activity with increased reaching distance (not shown). We suggest that correlation of cortical neuronal activity with reaching distance may depend on experimental constraints, e.g. reaching time limitations. Hence, the correlation of M1 neurons with distance emerges from the correlation between M1neurons and peak velocity or peak acceleration.

Discussion
Most motor cortical neurons were found to have specific preferred directions in which their mean firing rate during reaching movement was maximal [2,13]. Some studies showed that the activity of cortical neurons correlated with torques at joints and muscle tensions [41,42], or with other kinematic variables [13,19,43]. It has also been suggested that the activity of a particular cortical neuron may correlate with more than one movement parameter [13,14,18,44,45]. Although many previous studies addressed these phenomena (see for review [9,[46][47][48]), the origin of the directional tuning properties continues to be subject to debate. Neuronal activity in the motor cortex may reflect the complexity of movement dynamics [49] as well as higher-order features such as cognitive information related to a motor task or function [50]. This complexity has led to different viewpoints and interpretations [51].
During center-out reaching movements, average values of muscle forces and joint torques do not have pronounced directional dependence, so that anisotropy of motor cortical activity is truly non-trivial. The model presented here suggests that directional tuning in the primary motor cortex is the result of movement geometry and muscle dynamics. Arm muscles contract and stretch during reaching. Geometrically, the magnitude of muscle contraction or stretching depends on movement directions. As a result, muscle state variables such as muscle lengths and their rate of change (muscle velocities) have significantly different average values when the movement is performed in different directions. Our simulations have demonstrated that averages of muscle state variables fit perfectly the cosine tuning curves. This leads to similar directional tuning of afferent feedback and muscle contractile components, which are functions of muscle state variables.

Mechanisms of directional modulation
Most studies agree that M1 neurons are sensitive to movement directions, but fail to explain how directionally tuned motor commands produce non-directionally tuned endpoint kinematics (arm endpoint acceleration and velocity) and kinetics (net joint torques), or why nondirectional endpoint variables require directionally tuned motor commands.
Different factors contribute to directional modulation at different levels of the system's hierarchy. At the lowest level, frictional forces are proportional to joint angular velocities. The latter are controlled by muscle stretch/contraction and, hence, have directional modulation similar to muscle velocities. Therefore, since muscle forces have components working against frictional forces, they should exhibit some directional modulation. However, the resulting correlation between muscle forces and movement direction is pretty weak (Fig 12, blue bars) due to relatively weak friction in the joints. As an obvious implication, one can expect more pronounced directional modulation of muscle forces at higher viscous friction, e.g. if reaching movements are performed under water.
At the next level, we found that motoneurons showed stronger directional modulation than muscle forces (Fig 12, green bars). Muscle forces induced by inputs from motoneurons depend on muscle lengths and velocities (see Eq (4) in Methods). This implies that for muscle forces to be directionally indifferent, motoneuronal activities should have directional modulation opposite to those of contractile components of the corresponding muscles (see Eq (15) in Methods).
Further, motoneurons are influenced by afferent feedback from the muscles arriving at the spinal cord neural network. This feedback carries information about muscle lengths and velocities, which are strongly directionally dependent. That is, supra-spinal cortical signals have to be directionally tuned in the directions opposite to those of their corresponding feedback to provide appropriate activation patterns of spinal motoneurons. Interestingly, all factors described above are synergistic in a sense that all of them-frictional forces, muscle contractility, and afferent feedback-have the same directional modulations as the corresponding muscle lengths and velocities. This explains why directional tuning becomes more and more pronounced when ascending on the schematic in Fig 1A from net joint torques through muscle forces and motoneurons to the cortex (Fig 12).
Previous studies have also suggested that directional tuning in the motor cortex can be attributed to muscle-state variables [30], muscle-state-dependent variables [31,52], or limb biomechanics [27,32]. The hypothesis that directional tuning is the result of directional dependence of muscle length dynamics is not new or unique. However, our model provides additional supporting evidence and reveals several remarkable differences from previously published studies. Mussa-Ivaldi [30] described cortical activity as a linear function of muscle velocity, and suggested that directional tuning of cortical activity emerges from muscle state variables. Our results support this idea in general, but disagree with the notion that activity in the primary motor cortex is a simple linear function of the muscle velocity. Moreover, the linear relationship suggested by Mussa-Ivaldi [30], does not clearly determine whether cortical activity and muscle velocity have the same directional tuning properties or opposite ones. Our findings specifically show that motor cortical neurons and their corresponding muscle state variables, velocity and length, have opposite directional modulation. Moreover, this modulation nonlinearly propagates across the spinal cord network and affects all spinal neurons (see Fig 12).
Todorov, using a simplified mathematical model [31,52], demonstrated that activity in the primary motor cortex encodes arm muscle activity. He argued that correlation between cortical activity and high-level parameters, such as movement direction, speed, and target position, was a secondary effect. He suggested that cortical neurons are directionally tuned as a result of directional tuning of muscle state-dependent variables (forces generated due to muscle length and velocity). Similarly, contractile components in our model (force-length and force-velocity) contribute to the directional tuning of motoneurons and, thus, cortical neurons, which is in agreement with Todorov's findings. In addition, our model illustrates the interplay between cortical signals and afferent feedback, which leads to augmented directional tuning of neuronal activity in the motor cortex compared to the spinal cord. In contrast to Todorov's assumption that the spinal cord circuitry was not involved in directional tuning [31,52], our model provides evidence suggesting that spinal circuitry may play a significant or even major role.
The arguments above also explain why directional modulation of cortical neuron activity is independent of the biomechanical redundancy of the arm, which we tested by simulating reaching tasks with randomly distributing the torques between mono-and bi-articulate muscles. Indeed, additional directional modulation emerging at each level is solely defined by variable dependent on muscle lengths and/or velocities, i.e. by the kinematics of movement in a particular direction. The same reasoning explains the results of Lillicrap and Scott [32], who used an optimized artificial neuron network to show that limb geometry, intersegmental dynamics, and force-length/velocity properties of muscles are dominant factors in the emergence of directional modulation.
In light of the above, the directional modulation in the cortical activity can be thought as a product of motor learning which has interesting implications for the emerging neuro-prosthetic technology, a.k.a. brain machine interface (BMI). BMI utilizes neural signals of motor area in the brain to control external devices such as a robotic arm. Even though BMI decoders attempt to use natural neuronal tuning to actuate the movement, the initial motor performance is usually poor, but improves with training engaging the same neuroplastic mechanisms that are involved in motor learning. Since proprioceptive feedback responsible for the initial directional tuning is no longer relevant for the artificial actuator, one could predict that the neuronal tuning should change during the learning process to eventually reflect the BMI decoder properties. This prediction finds substantial experimental support and may explain difficulties in interpretation of neuronal tuning during BMI control (see [53] for review).

Directional modulation is invariant in the "right" coordinate system
We showed that transforming the workspace shifts the tuning curves (PDs) of M1 neurons and other directionally tuned variables. Similar changes in tuning properties have been observed in previous experimental studies. When the workspace was transformed counterclockwise from the left side to the middle and then to the right side of the animal's body, the PDs of cortical neurons shifted clockwise [38]. Tanaka and Sejnowski also reported shifts of PDs when the workspace was transformed [33]. Lillicrap and Scott [32] also showed that optimal preference movement distributions change as a function of limb posture and musculoskeletal organization. In addition, it was observed that directional tuning of motor cortical neurons was affected by the posture and orientation of the monkey's arm [23] even when the workspace remained fixed. Our analysis suggests that directional modulation remains virtually invariant with respect to the transformations above if defined in terms of the angle between reaching direction and the line connecting the shoulder joint and the initial arm position.

Directional tuning may vary depending on task conditions
It was proposed that directional tuning in the motor cortex, as defined by Georgopulos, is not robust; it varies depending on task conditions [54,55]. In particular, some studies showed that directional tuning changes with the limb's velocity or reaching distance [36,49]. Similarly, our simulations suggest that the fit of cortical activity by a cosine tuning curve gets less accurate when the movement speed is increased (not shown) which may result in significantly different (more variable) PD estimates. Neuronal activity patterns both in the cortex and the spinal cord have two components-the dynamic component aimed at generating a specific movement, and the compensatory component that counteracts the effects of directional modulation of muscle length and velocity dependent variables (frictional forces, contractile components, and afferent feedback). The latter is responsible of the directional modulation of neuronal activity. With increasing speed and/or distance of the movement, the dynamic component begins prevailing the compensatory component, and, thus, partially destroys cosine-like fits.

Higher-level and supra-spinal structures
The spinal cord circuitry we described in this study is the low-level circuit close to output motoneurons. In our model, we assumed that, M1 neurons directly control spinal motoneurons to activate muscles in the arm. Some pyramidal neurons in M1 of monkeys do directly project to motoneurons in the spinal cord [56][57][58]. However, a majority of M1 neurons indirectly control them through interneurons [59]. Various types of these interneurons with complex interconnections may be involved in relaying motor commands from the motor cortex to the spinal motoneurons. The Modularity Theory, for example, suggests that all motor behaviors may be constructed by combining motor primitives generated by composite motor modules in the spinal cord [60][61][62][63][64][65]. According to this theory, the brain does not control individual muscles directly, but rather activates sets of muscles via motor primitives to produce complex movements. Therefore, directional properties of the cortical neurons controlling motor primitives would depend on geometry of multiple muscles, and our finding that cortical neurons have directional preferences opposite to the corresponding muscle lengths and velocities may not be entirely valid. However, our main conclusion about directional dependence of the neuronal activity in the higher-level and supra-spinal structures due to directional dependence of the muscle geometry still holds.

Comparison with other models
Several previous models of arm reaching movements were developed to understand activity patterns in the motor control system and their relationships with movement parameters. In some studies, simple mathematical equations have been suggested to describe relationships between motor cortex activity and movement characteristics [30,66]. These models, however, were synthesized to illustrate particular directional tuning hypotheses, and did not have specific biological justification.
Activity patterns of motor cortical neurons have also been modeled based on endpoint velocity, acceleration and limb position to understand tuning properties derived from musclestate dependent variables [31,52]. However, spinal cord circuits and biomechanics of the arm were not considered in these studies. Tanaka and Sejnowski [33], for example, developed a three-joint arm model performing reaching movements in 2D space, where muscle forces were described as a rectified sum of motor cortical activities taking into account neither muscle contractile components nor spinal reflexes. Lillicrap and Scott [32] developed a two-joint arm model controlled by six muscles for reaching movements in 2D space, which is very similar to our arm's model. However, there are major differences between the model structures.
For example, their model did not include any spinal circuits, but included complex cortical neural computations to adjust parameters in the cortical network via a learning procedure using the feedback from the arm dynamics. We can argue that this learning procedure is aimed at finding a particular solution of the same inverse problem we are addressing in our study. However, the cost function used and a specific implementation of the learning process define which particular solution is selected. That may potentially induce the directional preference of its own. In our study, we analyzed a family of possible motor programs for the movements in each direction, and confirmed that all members of each family possess similar directional properties.
In summary, Lillicrap's and other models used optimal control theory based on specific assumptions that cortical activity is optimized to achieve certain behavioral objectives [32,34,67]. We do not use any assumptions about optimality. Instead, we focus on a particular type of movement, i.e. goal-directed reaching, explicitly assimilating information that reaching movements follow straight-line trajectories with a bell-shaped velocity profile as it was shown in previous experimental studies [36,68,69].

Conclusion
In the present study, we proposed mechanistic explanation of the fact that during planar reaching arm movements directional dependence of time-averaged activity of cortical neurons originates from the anisotropy of average muscle lengths and velocities. Specifically, cortical neurons have opposite preferred directions compared to the "preferred" directions of lengths and velocities of corresponding arm muscles. This anti-correlation builds up as one ascends along the hierarchy of the motor control system. First, the directional preference emerges at the level of muscle forces as a consequence of viscous friction in the joints. Second, it augments at the level of spinal motoneurons to compensate for the muscle contractile component (which are functions of the muscle lengths and velocities) directional dependence. Third, neuronal directional modulation is further amplified at the cortical level to counteract the directional dependence of afferent feedback inputs carrying information about muscle lengths and velocities to the spinal cord network. To conclude, we showed that the motoneurons and afferent feedback have significant impact to modulate the directional tuning behavior of cortical neurons.

Biomechanical model of the arm
The Arm (see Fig 1C) is modeled as a 2-joints mechanical system of rigid segments connected by revolute joints, which operates in the horizontal plane. The dynamics of the arm's motion are derived from the Lagrange equations and include angular velocities and accelerations, Coriolis and centrifugal forces, muscle forces, and viscoelastic forces at the joints. The parameters of mechanical system such as masses and lengths of correspondent arm segments are based on human biomechanics [70]. Kinematics of the segments is described by the following differential equations: where: q 1 = q 1v − q 2v + q 1M and q 2 = q 2v + q 2M are generalized forces (torques), which include joint friction forces (q 1v , q 2v ) and torques created by muscles (q 1M , q 2M ); θ 1 , θ 2 are the generalized coordinates (the angles between positive y-axis and the upper arm or forearm, respectively); I 1 ¼ m 1 L 2 1 =3 is the moment of inertia of the upper segment around an axis through the point of suspension; I 2 ¼ m 2 L 2 2 =12 is the moment of inertia of the lower segment about an axis through its center of mass; m 1 = 1.79 kg, m 2 = 1.55 kg are the masses of the upper and lower segments, respectively; L 1 = 0.34 m is the length of the upper segment, L 2 = 0.31 m is the length of the lower segment; d 1 = L 1 /2 is the distance from the shoulder joint to the center of mass of the upper arm, d 2 = L 2 /2 is the distance from the elbow joint to the center of mass of the forearm.
After simplifying Eq (1), the movement of the 2-joint arm can be described by angular accelerations at the shoulder and elbow joints given by: where: The joint viscous friction forces are given by: where the viscosity parameter η v = 0.05.

Muscles
We considered six muscles, which control arm movements in horizontal plane. Specifically, four 1-joint muscles (shoulder flexor (SF), shoulder extensor (SE), elbow flexor (EF) and elbow extensor (EE)) and two 2-joints muscles (shoulder and elbow extensor (BE) and shoulder and elbow flexor (BF)) control the arm's movements. Muscles, which perform the similar action and have similar anatomy, were substituted by one muscle with averaged parameters. Total muscle torques, q 1M and q 2M , are produced by all muscles at the shoulder and elbow and calculated as: where: F and R represent muscle forces and averaged moment arms [71,72] for all muscles (see Table 1). Muscle indices are defined as follow: SF and SE represent the shoulder flexor and extensor; EF and EE represent the elbow flexor and extensor. BFS, BES represent 2-joint flexor; BES, BEE represent 2-joint extensor, respectively.
Description of muscle forces is based on the Hill-type model proposed by Harischandra and Ekeberg [73] was also adapted to human biomechanics. The total force F for each muscle is calculated as: where: F max is the maximum force (see Table 1) based on experimental data [71,72]; MN is the   Table 2.

Afferent feedback
Muscle afferent feedback activities (Ia and Ib feedback signals) are derived and modified from the formulas suggested by Prochazka [75]: where: v norm = v/V max is the normalized muscle velocity, d norm = (l − 0.2 Á L opt )/L opt is the normalized muscle lengthening if l > L opt or 0 otherwise, y is the output activity of the corresponding motoneuron, F norm = (F − 0.1 Á F max )/F max is the normalized muscle force. Parameters V max , L opt and F max are based on averaged parameters of correspondent muscles of upper of human male dataset [71,72] and specified in Table 1. The coefficients k v are optimized in order to normalize Ia signal from each muscle operating within locomotor range (see Table 1). The remaining coefficients are defined as follows: k dI = 0.8, k nI = 0.05, k F = 1.0, const I = 0.01.

Spinal neuronal network
The Spinal Cord circuit provides direct control of the arm biomechanical model via corresponding motoneurons (Fig 1A). Although the activity of the motor cortex is correlated to muscle activity during voluntary movements such as reaching task [9][10][11], spinal reflexes still play an important role in arm kinematics [76]. In our model of the spinal cord we implemented several basic reflexes such as: stretch reflex; autogenic inhibition reflex; and recurrent inhibition of motoneurons via Renshaw cells (see Fig 1B).

Muscle
Length, l (m) Velocity, v (m/sec) A simplified scheme of the stretch reflex including monosynaptic excitation of synergist motoneurons and disynaptic Ia reciprocal inhibition circuitry is based on previously published work [77]. The additional Ia-interneurons, which receive Ia afferent input and mediate inhibition to antagonist motoneurons, are introduced for all pairs of antagonist muscles. These interneurons also receive mutual inhibition from antagonist Ia-interneurons.
Because of autogenic inhibition reflexes, activation of Golgi tendon organ Ib afferents evokes inhibition of synergist motoneurons [78,79]. The inhibitory pathway includes an additional Ib-interneuron that receives an excitatory Ib feedback from the corresponding muscle and relays the inhibitory signal to the same motoneuron and its synergists. Schematic of this reflex is based on previously described circuitry [77,80].
The activity of motoneurons is also regulated by Renshaw cells [81]. These interneurons innervate and inhibit the very same motoneurons which activate them as well as synergetic motoneurons [82]. Renshaw cells also inhibit ipsilateral interneurons which provide Ia reciprocal inhibition [78]. In addition, the activity of Renshaw cells is regulated by the activity of contralateral Renshaw cells excited by antagonist motoneurons (see Fig 1B).
Each neuron in our model represents the activity of a population of corresponding spiking neurons and is based on a non-spiking description of the neuron model. The output activity (firing rate) of each neuron (y) is represented by a sigmoidal function of the aggregate input: x i is the i-th input and w i is the synaptic weight (or strength of connection) between the i-th input and the neuron, b is a bias term that controls the excitability of the neuron, and s(v) is an activation function. The parameters v 1 = 2 ¼ 0:5 and k = 0.1 in Eq (6) specify the threshold and the slope of the activation function, respectively. As mentioned above, the spinal cord neuronal network consists of six local circuits which are responsible for controlling corresponding muscles. Each circuit includes a motoneuron (MN), Renshaw cell (RC), Ia-interneuron, and Ib-interneuron. The bias term was the same (b = −0.28) for all neurons in the spinal cord.
The gains of the feedback inputs were chosen in such a way that all reflexes modulate the motoneuron activity within 15% of maximum. This roughly corresponds to existing experimental data on contribution of reflexes in the EMG activity during voluntary limb's movements [83,84]. The connections in the network are defined as follows [80] (see also Fig 1B). The synaptic weights are indicated in brackets.
• MNs and Ia interneurons are excited by corresponding Ia feedback (+0.15).
The spinal cord neuronal network contains recurrent connections, and, hence, the firing rates cannot be calculated in a feed-forward way. We assumed that the firing rate variables of the network are in instantaneous equilibrium, which can be found by an iterative procedure explained below. The "current" state of the network, Y ! k , that includes spinal motoneuron activities as well as firing rates of all spinal interneurons, is described by: where W is the matrix of connection weights between elements (including ascending feedback, neuron outputs, and descending signals from upper levels, see above) in the network, vector B ! consists of bias terms for all neurons in the network, k is the iteration number, MC ! is the input from the motor cortex, FB ! is a vector of afferent feedback from the muscles, and s is the sigmoid activation function (6). Eq 7 is iterated until Y ! k reaches a steady state with a preset tolerance. Its equilibrium value, Y ! , is accepted as the network response to the cortical input MC ! given the feedback FB ! . We verified that with the synaptic weights used (W) the map (7) always had a unique stable equilibrium. Hence, the network state can be considered as a function of the cortical inputs and afferent feedback signals: calculated by iterating (7).

The cortical controller
We model the execution of a reaching task by moving the arm endpoint (the wrist) from a given initial position to a desired target position. Based on the specified target position, the model calculates a set of activity profiles for the motor cortical neurons (a motor program), which provide the wrist movement along a defined trajectory, ending at the target position. Given the initial and target positions of the arm's endpoint in the (x, y) orthogonal coordinate system, and preset reaching time, we first calculate muscle forces required to generate the desired motion along a straight-line trajectory with a defined velocity profile. Using the muscle forces, we then calculate required motoneuron activity (motoneuron signals), and finally supra-spinal input generated by motor cortical neurons.

Arm trajectory and joint angles
Let the initial and target positions of the arm's endpoint be (x 1 , y 1 ) and (x 2 , y 2 ), respectively. The velocity profile along the trajectory, v(t), of the arm's endpoint is defined as follows: where T is the reaching time (time to reach the target), L is the distance between the initial and target positions in meters (reaching distance), v(t) is in m/s and time t is in seconds. The velocity has a bell-shaped profile and is zero at the initial and target positions. The peak velocity is controlled by T and L; the inverse relationship between reaching time and peak velocity has been experimentally reported [85]. Using (9), we can calculate the dependence of the endpoint coordinates on time as follows: Here U x and U y are the components of the unit vector along the trajectory. We then calculate the coordinates of the elbow joint (x e , y e ) by finding the intersection of two circles with centers at the origin (shoulder) and at the endpoint of the arm with radii L 1 and L 2 , respectively (see Fig 1C): Given the time dependencies of the wrist and elbow coordinates, we calculate the joint angles and their first and second derivatives using these simple geometrical relationships: x e ¼ L 1 siny 1 ; y e ¼ À L 1 cosy 1 ; x À x e ¼ L 2 siny 2 ; y À y e ¼ À L 2 cosy 2 ð12Þ by differentiating Eqs (10), (11) and (12) twice with respect to time and solving the resulting equations for y 1 ; y 2 ; _ y 1 ; _ y 2 ; € y 1 ; € y 2 . We omit the explicit formulas due to their complexity. As a result of these procedures, we obtain the values for the joint angles and their first and second derivatives (angular velocities and angular accelerations) at every time t during the reaching movement.

Muscle forces
Using Eq (2) and the angular accelerations, we calculate total torques q 1 and q 2 at the shoulder and elbow joints: which allows us to calculate the torques generated by the muscles after subtracting frictional torques: where q 1,M and q 2,M are the total muscle torques at the shoulder and elbow joints, respectively; q 1,v and q 2,v are the frictional torques in the shoulder and elbow joints, respectively, that are defined by angular velocities (see Eq (2) and the text below it).
Since torque values (q 1,M , q 2,M ) are created by 6 muscles, there may be significant redundancy in possible muscle activation patterns. To limit possible solutions, we made the following assumptions: 1) positive torque values correspond to flexor activity, and negative torque values correspond to extensor activity; 2) since bi-articular muscles are concurrently active with shoulder and elbow muscles, there are infinitely many ways to distribute the load over the muscles to provide the same aggregate torque. Therefore, we introduced a dimensionless control parameter d which dictates how the torque at the shoulder joint is distributed between the bi-articular muscles and single flexor and extensor muscles. Thus, muscle forces were calculated using the following equations: where d represents torque distribution about the shoulder, and other notation are similar to those in Eq 6. We chose to include the torque distribution parameter at the shoulder joint, rather than at the elbow, because the shoulder joint experiences significantly larger torque loads compared to the elbow joint; 3) the calculated forces cannot exceed the maximal possible force values associated with the muscles (see Table 2). The above assumptions allow us to find one parameter family of solutions for the system (3) (a family of motor programs) for all six muscle forces. Parameter d defines the bi-articular muscles' level of participation, such that d = 1 indicates that the torque at the shoulder is fully generated by the single-joint shoulder muscles, while d = 0 indicates that the torque at the shoulder is fully generated by the bi-articular shoulder muscles. In all the performed simulations, the torque distribution parameter d was randomly varied in the interval (0.5, 1) (unless otherwise stated), using a uniform probability distribution, to construct an ensemble of possible motor programs for every particular reaching movement.

Motoneuron pool signals
Activities of six motoneuron populations (MN) corresponding to and controling the six muscles actuating the biomechanical arm is given by: where F is the muscle force calculated in the previous step, and F max , F p , F l , and F v are muscle specific parameters explained above (see Eq (4) and comments below the equation).

Activities of cortical neurons
The six cortical neurons form the supra-spinal input drives (MC ! ) to the motoneurons. Calculation of their desired activity patterns is achieved by numerically inverting the dependence of motoneuron activity on cortical inputs (8) using the adapted secant method implemented as the following iterative procedure: where j ¼ 1; 6 is a component number, i is the iteration number, and (8)). The map (16)

Simulation setup and regression analysis
Following many previous experimental studies, the movement direction was defined as the direction from the initial center position to a peripheral target position [2,13,14]. The 0 o direction is defined as moving the arm's endpoint (wrist) to the right along a horizontal line, and moving the wrist in the opposite direction is the 180 o direction. Moving the wrist away from the body along a vertical line is the 90 o movement direction, and moving the wrist in the opposite direction is the 270 o movement direction. To understand the relationship between cortical activity and movement direction, our arm control model simulated several center-out reaching movements with an identical initial center position and equal reaching distances (radii) to 8 peripheral target positions in 8 different directions, at 45 o intervals from each other (Fig 2). The same type of center-out task has been used extensively in experimental studies [2,13]. Average cortical activity, average muscle length and velocity, and average feedback signals were fitted to a cosine tuning curve using a regression function. The cosine tuning curve is given by: which is equivalent to: where b 0 , b 1 , b 2 , and c 1 are regression coefficients. θ PD is the preferred direction in which the function has a maximal value. The coefficient of determination R 2 was used to measure the degree of the regression fit. Note that the average of f across all directions (8 directions in this case) is equal to b 0 . The same regression analysis of cortical activity using a cosine tuning curve was previously used in several experimental studies [2,13] to investigate the directional tuning of cortical neurons. The index of directional modulation, which describes the proportional increase or decrease over the mean activity level, is given by I = c 1 /b 0 . For most of the center-out reaching tasks, except where noted, the following simulation parameters were used: integration time step of 0.001s; reaching time of 1s; reaching distance of 0.2m; and initial center position with Cartesian coordinates (0.0, 0.4). Note that we further investigated the center-out tasks with 16 different directions at 23 o intervals and found no significant differences in the obtained results. Therefore, we elected to use 8 directions for our study.

Model validation
In order to evaluate the robustness of the model we independently varied basic parameters of the mechanical system (masses and lengths of arm segments) within ±15% range. In addition, we investigated the model's behavior by varying maximal forces (F max ) for each muscle in the same range (±15% around mean). For all simulations, there were no qualitative differences in the model's performance. Specifically, the tuning curves of the cortical activities were similar to the corresponding tuning curves for the default values of the biomechanical variables. The statistical analysis based on the 50 simulations also did not show any significant differences in PDs for the chosen ranges of the basic biomechanical parameters (see Table 3).

Software
The model of the motor control system was implemented under MATLAB 8.4. The code used for the simulations presented here will be made available. Supporting information S1 File. Directional modulations and tuning curve fittings. The compressed zip file contains the data in '.mat' format and the Matlab codes which can be accessed by MATLAB. The data and codes in this file produce reaching movement in 8 directions, directional modulation curves, tuning curves, preferred directions and coefficient of determinations. (ZIP)