Muscles Recruited During an Isometric Knee Extension is Defined by Proprioceptive Feedback

Proprioceptive feedback and its role in control of isometric tasks is often overlooked. In this study recordings were made from upper leg muscles during an isometric knee extension task. Internal knee angle was fixed and subjects were asked to voluntarily activate their rectus femoris muscle. Muscle synergy analysis of these recordings identified canonical temporal patterns in the data. These synergies were found to encode two separate features: one concerning the coordinated contraction of the recorded muscles and the other indicating agonistic/antagonistic interactions between these muscles. The second synergy changed with internal knee angle reflecting the influence of afferent activity. This is in contrast to previous studies of dynamic task experiments which have indicated that proprioception has a negligible effect on synergy expression. Using the MIIND neural simulation platform, we developed a spinal population model with an adjustable input representing proprioceptive feedback. The model is based on existing spinal population circuits used for dynamic tasks. When the same synergy analysis was performed on the output from the model, qualitatively similar muscle synergy patterns were observed. These results suggest proprioceptive feedback is integrated in the spinal cord to control isometric tasks via muscle synergies. Significance statement Sensory feedback from muscles is a significant factor in normal motor control. It is often assumed that instantaneous muscle stretch does not influence experiments where limbs are held in a fixed position. Here, we identified patterns of muscle activity during such tasks showing that this assumption should be revisited. We also developed a computational model to propose a possible mechanism, based on a network of populations of neurons, that could explain this phenomenon. The model is based on well established neural circuits in the spinal cord and fits closely other models used to simulate more dynamic tasks like locomotion in vertebrates. Conflict of interest statement The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Introduction
Execution of a motor task is considered to be modular in nature and is modified by sensory inputs from the periphery and descending input from the brain. The role of proprioceptive feedback in the effective recruitment of muscle fibres to counter load experienced during a given task is well studied. However, its role in control of interactions between muscles is still poorly understood during an isometric tasks. In general it is assumed that proprioceptive feedback is only marginally involved in isometric tasks, especially if static and at a single joint. In this paper, we will show that this assumption is unwarranted, along with a mechanism of action for its role in shaping the recruitment of the muscles.
There is great variation in the way a motor task can be performed, even at a single joint. Each variation is produced from a combination of muscle recruitment patterns, often described in terms of time. These patterns are commonly referred to as muscle synergies. Their use by the central nervous system (CNS) to alleviate the degrees of freedom problem is accepted but little is known about the mechanism of their recruitment (Grillner 1985;Bizzi, Mussa-Ivaldi, and S. Giszter 1991;Tresch, Saltiel, et al. 2002). Similar synergies are reported across species, especially for routine repetitive tasks like locomotion in vertebrates (Yang, Logan, and S. F. Giszter 2019), where antagonistic pairs are recruited at and across joints to generate a coordinated alternating pattern of activity (Dominici et al. 2011). Although not all neurons associated with this recruitment pattern have been identified, their functional grouping into populations within the spinal cord is undisputed, as is their ability to produce all locomotor output patterns observed under fictive conditions (Martin et al. 2007). The current accepted model based on experimental results, is a three layered central pattern generator (CPG) (McCrea and Rybak 2008;Rybak, Dougherty, and Shevtsova 2015). Although this simulation can reproduce the rhythm, pattern and output from motoneurones during locomotion (Markin et al. 2016;Shevtsova et al. 2016), the component parts of this network model have rarely been applied to human activity. In this study, we report experimental data that clearly shows an influence of proprioceptive feedback on the observed synergies. Based on our experimental data we propose a model for isometric tasks by humans that uses key elements of the model proposed by Rybak's lab (Rybak, Dougherty, and Shevtsova 2015). This model is described at the level of populations of neurons and is simulated using population density techniques (PDTs) which have been shown to accurately model population aggregates (like firing rates) while retaining a close correspondence to spiking neurons -more so than neural mass models -without producing the overhead of simulating thousands of neurons (De Kamps, Lepperød, and Lai 2019).
To investigate the effect of afferent feedback on muscle synergy recruitment, we examined the change in interactions among a group of selected muscles, at fixed knee joint angles during an isometric knee extension. The interactions between homonymous and heteronymous muscles are well known (Pierrot-Deseilligny and Burke 2005) but rarely examined in the same isometric task. We recorded from seven muscles in healthy young subjects while they performed an isometric knee extension focusing on voluntarily activating the Rectus femoris muscle of the quadriceps. We hypothesised that this would result in the same synergies being recruited at all angles of the task but changes in synergy recruitment would be due to static proprioceptive feedback from muscle stretch. Observed changes to the recruitment of the quadriceps and hamstrings, in combination with our model, suggest that afferent feedback does affect the recruitment of muscle synergies during an isometric task. This clearly demonstrates that assumptions regarding proprioceptive feedback in motor control have not been sufficiently investigated.

Ethics
The study was conducted according to the Declaration of Helsinki and all experimental protocols were approved by the University of Leeds Research Ethics Committee (reference number BIOSCI 16-004). 17 healthy subjects of mixed gender (male = 9, female = 8) with an age range of 18-30 (24.4± 2.57 years) were recruited to participate in this study. Exclusion criteria included previous knee or leg injuries, if participants had done exercise within 48 hours prior to testing, knee stiffness or self-reported pain, use of recreational or performance enhancing drugs, ingested alcohol in the previous 24 hours or were unable to provide informed consent. Subjects provided informed written consent to the study, noting possible risks associated with the activity.

Data Collection
Surface EMG was recorded from seven muscles of the subjects dominant leg; rectus femoris (RF), vastus lateralis (VL), vastus medialis (VM), semitendinosus (ST), biceps femoris (BF), medial gastrocnemius (MG) and tibialis anterior (TA) -of which the MG and TA were discarded due to low signal to noise ratio. Data analysis was therefore performed on the five remaining muscle recordings. The skin was prepared for electrodes with shaving, cleaning with alcohol wipes and then application of conductive electrode gel. Data was sampled at 2 KHz using wireless Delsys Trigno IM electrodes. Electrodes were placed on the muscle belly, defined by landmarks as described in Rainoldi, Melchiorri, and Caruso 2004 based on anatomical observations: VL-between the greater trochanter and the lateral epicondyle; VM -on the distal fifth of the medial knee joint; RF -between the greater trochanter and the lateral epicondyle; VM -on the distal fifth of the medial knee joint; RF -between the anterior superior iliac spine and the superior pole of the patellar, and MG belly located in distal third of the medial knee joint.

Movement Protocol
Subjects were asked to lay on a standard medical examination bed. They were then shown how to perform an isometric knee extension with the leg brace attached to their dominant leg. Subjects were shown the resulting EMG output recorded using a Delsys Trigno system. Subjects were asked to perform an isometric knee extension at maximal voluntary effort for five seconds, attempting to maximise RF activity. This was repeated six times with a three minute rest between contractions. The dominant knee was fixed at one of four angles using a Donjoy TROM locking knee brace at 0°, 20°, 60°and 90°. The angle of the knee was always measured against the hip joint and the bony prominence on the outside of the ankle. Data was collected in two different positions and sessions for each subject. In position one the participant was supine with both legs flat against the bed. In position two the contralateral leg was kept bent such that the foot is flat against the bed so that both the knee and hip are fully flexed. The position selected for each subject was randomized for their first session. In the second session the subject performed the task in the other position.

Data Preprocessing
Data was prepared for NMF via a three stage filtering process. Signals were initially band-pass filtered (high pass = 20 Hz, low pass = 450 Hz, second order Butterworth filter), rectified and then finally zero-lag high-pass filtered (5 Hz, second order Butterworth filter) to remove frequency changes induced by rectification. Each EMG channel was normalized to the maximum value for that channel across all six contractions. Data was visually inspected and segmented into equal sections containing one burst. NMF was performed on all contractions across all EMG channels.

Synergy Extraction
We used NMF to identify muscle synergies during the task. Information theory shows that this dimensionality reduction reflects latent structure in the data, in this case muscle synergies (Lee and Seung 2001). NMF's chief advantage compared to other approaches is the constraint of non-negativity aligning with muscle activity i.e. muscle activation is never negative. NMF is also more effective at identifying latent structure in the data when compared to other techniques such as principal component analysis (Ebied et al. 2018). The five raw (but smoothed) EMG time series were combined into a matrix D of size 5 × n where n is the length of the time series. We used iterative NMF decomposition algorithms (Tresch and Bizzi 1999;Lee and Seung 2001) to reduce D to a combination of two matrices, W and C such that, C is an N × n matrix where N is the chosen NMF rank factor, in this case 2. Each row of C represents some structure in the time series similar to a PCA component. W is a 5 × N matrix which, when multiplied by C, approximates matrix D. Each column of W quantifies the amount that the corresponding row in C contributes to the original data in D (Lee and Seung 2001;Donoho and Stodden 2004;Berry et al. 2007;Torres-Oviedo and Ting 2007). Each synergy, s, is represented by the corresponding column W * s and row C s * . We describe C s * as the activation pattern of the synergy as it represents some underlying structure of the original EMG time series. We refer to W * s as the muscle contribution vector of the synergy as each component value indicates the contribution of the synergy's activation pattern to the associated muscle activity.
Selection of rank factor is critical to produce dimensionality reduction such that C has fewer rows than D. Rank factor was chosen consistent with previously literature such that rank factor was increased to the minimum required for the variance accounted for (VAF) by W C compared to D was greater than 90%. VAF was calculated for each synergy profile for both the individual muscle and for all muscles collectively. If VAF was below 90%, the resulting synergies were discarded (Tresch, Cheung, and d'Avella 2006). The iterative optimization algorithm used was initialized using singular value decomposition to reduce calculation time and to ensure a unique result (Boutsidis and Gallopoulos 2008).
Each row in C and column in W was normalized to its maximum value. Cosine similarity analysis was used in a pairwise fashion to determine the similarity between subjects' synergy vectors and activation coefficients W * s and C s * (Rimini, Agostini, and Knaflitz 2017).

Neural Population Modelling
Our aim was to create a neural population model such that applying NMF to the firing rate activity of the motor neuron populations would yield the same synergy patterns as those identified in the EMG data. We did not attempt to reproduce simulated EMG signals. Instead, we assumed that the cumulative activity of multiple motor units described by the average activity of distinct motor neuron populations would serve as a proxy for EMG. We first considered rate-based models which represent a population metric, for example the average firing rate (Wilson and Cowan 1972) or oscillation frequency (Kuramoto 1991), abstracted from the underlying individual neurons. Rate-based models are suitable for reproducing firing rates in neural circuits, but there is no clear relationship with the state of the underlying neural substrate. Although not essential for this study, in light of detailed modeling with Rybak's group, as well as future development of the modeling work, we are interested in a technique that retains a closer relationship with the state of the spiking neurons that comprise the neural circuit. PDTs are such a technique: they retain information about the state of neurons in the circuits but calculate population level aggregates directly.

Population Density Techniques
Population Density Techniques (PDTs) model neural circuits in terms of homogeneous populations of neurons. The individual neurons are described by a model, such as leaky-integrate-and-fire. The model of an individual neuron is characterised by a so-called state space: the values that determine the state of individual spiking neurons. For a simple neuron model this can be its membrane potential. For more complex models variables representing the state of e.g. a synapse can be included. PDTs represent a population by a single density function that represents how neurons are distributed across the neuron's state space.

MIIND
MIIND is a neural simulator De Kamps, Baier, et al. 2008 which implements a version of a PDT to simulate multiple interacting populations of neurons. It can provide a visual representation of the probability density function by displaying the density during simulation. Figure 5 shows an example of this visual representation. A network of populations can be built in MIIND using a simple XML style code format to list the individual populations and the connections between them. Populations in the network interact via their average firing rates, which are assumed to be Poisson distributed spike trains. For each connection, the firing rate of the source population becomes the average rate of the Poisson distributed input spikes to the destination population. The connections defined in the XML code, have three parameters: the post synaptic potential or instantaneous synaptic efficacy, the number of individual connections between source neurons and target neurons, and a delay which can be used to approximate time taken for spike propagation and synapse transmission.

The Spinal Circuit Model
We used MIIND to build a network of populations of exponential integrate and fire (EIF) neurons according to the connectivity diagram in Figure 1. Table 1 shows the connection parameters for all populations in the model. All populations use the same underlying neuron model as described in Equation 2.
Where v is the membrane potential, v rest = −70 mV, ∆ T = 1.48, v thres = −56 mV, and τ = 3.3 ms. The parameters were chosen so that populations could produce a wide range of average firing rates between 0 and 200 Hz to exhibit typical neuronal frequencies. We chose to use an EIF model in contrast to more commonly used Hodgkin Huxley style neurons. This is because the objective was not to reproduce the EMG signals exactly, but to provide a concise explanation for the overall synergy patterns. We expected that any particular description of activation of ion channels (as  in a Hodgkin Huxley style model) would have no significant impact on the population level activity or synergy patterns in this task and would therefore dilute the power of the model. The main structure of the network consists of two neural populations, named "Extensor Interneurons" and "Flexor Interneurons", connected together in a network with five motor neuron populations, one for each muscle. The Extensor and Flexor Interneuron populations represent combinations of excitatory and inhibitory neurons and therefore can project both kinds of connections to other populations in the network. Other studies have previously described the connection motif of agonist inhibition with antagonist excitation and this is utilised here to connect the interneuron and motor neuron populations to elicit the agonist/antagonist relationship between the five muscles (Sherrington 1909;Doss and Karpovich 1965;Bigland-Ritchie 1981;Pierrot-Deseilligny and Burke 2005). These features, including the mutual inhibition between the two interneuron populations, also appear in the central pattern generator (CPG) model of McCrea and Rybak 2008. Although the rhythm and pattern formation parts are not included, the implications for applying a CPG model to an isometric task are discussed later. All supraspinal activity comes from the Cortical Drive input and is responsible for the "contraction". There is a direct connection to the MN-RF motor neuron population indicative of the muscle which is being maximally contracted in the task. Cortical Drive also projects to the Extensor and Flexor Interneuron populations. As there are more excitatory than inhibitory connections from the Extensor and Flexor Interneuron populations to the motor neuron populations, the Cortical Drive indirectly causes excitation of all motor neuron populations as well as MN-RF. During the simulation, the input to the two interneuron populations begins at 0Hz before increasing to 260Hz over 1 second, then five seconds later, dropping back to 0Hz over 1 second. We hypothesised that the most important factor in shaping the observed synergies would be the connectivity of the network model and that they could be modulated with a proprioceptive input. To simulate changes in proprioceptive feedback due to the knee angle, the Extensor Afferent Feedback input to the Extensor Interneuron population was introduced and for different trials was altered between 110Hz and 180Hz. The Flexor Afferent Feedback input was held constant. The average firing rate of each of the five motor neuron populations was recorded at a rate of 10 KHz (corresponding to the 0.1ms time step of the simulation) then sampled at 2ms intervals. NMF was performed on the resultant time series as described for the experimental recordings. Finally, we augmented the network (The greyed area of Figure 1) to generate a knee flexor bias for RF and an extensor bias for ST which was observed in the experimental results. An additional excitatory connection was added to the model from the Extensor Interneuron population to MN-ST, and from the Flexor Interneuron population to MN-RF. This is equivalent to increasing the number of excitatory connections overall between those populations. In order to modulate the effect of the afferent feedback input on these connections, two additional populations of inhibitory neurons were added to the model: InhibST and InhibRF. This network motif of an additional excitatory drive coupled with a controllable inhibitory input has previously been used to reproduce observed activity in Semitendionsus and Rectus Femoris of a cat (Shevtsova et al. 2016) and further supports the use of CPG models for human studies.
In the task experiment, two positions were used to identify the effect of passive insufficiency on the synergies. We expected, however, that the contralateral hip position would serve only to influence the degree of muscle stretch which in the model would already be accounted for in the activity of the afferent feedback input. There is, therefore, no analogue to position defined in the model.

Statistical Analysis
All statistical analysis was performed in Python 3.6.2. Cosine similarity analysis was used to compare EMG profiles, synergy activation patterns and muscle contribution vectors. Cosine similarity analysis is sensitive to differences in vectors that may have equal variation.

Code Accessibility
NMF analysis and cosine similarity analysis was performed using a custom designed program in Python 3.6.2. MIIND is available at http://miind.sourceforge.net/ and the model files and simulation results are accessible at https://github.com/hugh-osborne/isotask. * During the task, Cortical Input transitions from 0Hz to these values then back to 0Hz ** Afferent input remains constant throughout the activity but is set between 110Hz and 180Hz to produce different synergy patterns.

Results
We recorded from 7 different muscles of the leg, but only 5 of these were used for further analysis of activity and patterns as on examination the muscles TA and MG were always inactive, as expected due to the nature of the task. We extracted muscle synergies from the EMGs recorded from the 5 different muscles across the two different positions to examine how proprioceptive feedback alters muscle synergy recruitment. The contralateral hip was flexed or relaxed to induce passive insufficiency in the recorded leg to highlight differences in proprioceptive feedback. A spinal population circuit model was created with connections between interneurons, motor neurons and afferent feedback, based on current CPG models and accepted neural circuits (Pierrot-Deseilligny and Burke 2005). Similarities between the synergies observed experimentally and those extracted from the model were compared. Figure 5 shows the EMG signals for all five muscles after rectification and smoothing, and the average firing rates of the motor neuron populations in the simulation. It is difficult to discern the meaning of differences between the time series, highlighting the need for analysis techniques such as NMF.

NMF identifies two muscle synergies from the EMG activity
To identify synergies appropriate for experiment-model comparison, NMF was performed with a range of rank values, one rank per synergy. The appropriate rank to use was chosen as the number required to raise the VAF above 90% (Figure 2). In this case rank two raised VAF above this threshold. Although 90% is an arbitrary threshold, and there are other methods for choosing appropriate rank, patterns identified by three or more synergies were less consistent across participants. As described in section 2.4, each synergy consists of a column of matrix W with length five (one value per muscle) and a row of matrix C representing a time series describing some underlying structure of the original data. For each muscle, the corresponding component of W * s multiplied by C s * , gives the contribution of synergy, s, to that muscle's EMG. Cosine similarity analysis was performed on the synergy rows and columns across participants for each position, synergy and angle. There is high correlation between synergy 1 results among the participants, regardless of position and internal knee angle (Table 2). Though not as high as synergy 1, there is also high correlation between participants for synergy 2. Despite some variation, there is a common pattern of muscle synergy recruitment across all participants.
Figure 2: Average VAF scree plot for rank one to five NMF dimensionality reduction across all angles and both positions of the isometric knee extension task. The 90% VAF threshold indicates that two is the appropriate rank to use and therefore the number of synergies to extract. Error bars show Standard Error of the Mean (SEM).  Table 2: Synergy rows and columns (as defined in section 2.4) were compared across all pairs of participants using cosine similarity analysis giving a value between 0 (uncorrelated) and 1 (Highly correlated). For both positions (activating or inactivating the contralateral hip flexors) and for all internal knee angles, there is high correlation between subjects indicating that, during the task, the same synergy patterns are being recruited by the majority of subjects.

Proprioceptive feedback alters muscle synergy recruitment in a manner reflecting muscle stretch
The NMF process generates, for each of the two synergies, a time series activation pattern and a vector of five values, one for each muscle. Figure 3 shows the vector and time series of synergy 1 (A) and 2 (B) for both positions across different internal knee angles. The activation patterns (line plots) should be considered in conjunction with the five value muscle contribution vectors shown in the bar charts. Synergy 1 represents general coordinated muscle recruitment as would be expected in an isometric extension and contributes to the majority of the observed EMG activity. Because of this, the activation pattern closely matches the overall profile observed in the raw EMG data (the transition from low to high to low activity during the contraction). The high muscle contribution values for all five muscles indicates that this activation pattern is present in all five EMG recordings. Both the activation pattern and muscle contribution weights are well conserved across all angles, positions, and muscle groups. Among the muscle contribution values, there is a slight shift between 0°and 90°from a quadriceps bias to a hamstrings bias, however this change is less apparent in position 2. The shape of the activation pattern of synergy 2 indicates that this synergy captures a difference between the activity of each muscle during the resting and transition periods at the beginning and end of the contraction. During the contraction, the contribution of synergy 2 reduces to near zero because all muscles reach their maximal activity which is the feature captured solely by synergy 1. In position 1, where both hips are extended, synergy 2 muscle contributions show a strong hamstring bias at 0°which reduces as knee angle is increased to 90°. This shows that the synergy 2 activation pattern contributes strongly to the hamstring muscle EMGs at low internal angle but contributes more evenly across all muscles as angle is increased. In position 2, where the contralateral knee and hip are fully flexed, the same strong bias is indicated at 0°but is eliminated once the knee is set at 20°. The RF and ST muscles have higher contributions in synergy 2, when compared to the other quadriceps and hamstrings respectively. At 0°, 20°and 60°in position 1, and at 0°in position 2, the ST component is consistently higher than BF. Across both positions and all angles, RF has a higher component than VM and VL.

The model reproduces synergies observed experimentally
During each MIIND simulation, the two afferent feedback inputs were kept constant and the Cortical Drive input was changed from low to high activity to reproduce the contraction behaviour. All populations produced average firing rates which were either passed to connected populations in the network or recorded for analysis. The activity of the five motor neuron populations, MN-RF, MN-VL, MN-VM, MN-ST and MN-BF, was analysed. The raw output from these populations is shown in figure 5 (top right). The output is a great deal smoother than the overlaid EMG recording data due to MIIND's simulation technique and the lack of many of the experimental sources of noise. Though there is undoubtedly a great deal more information available in the EMG traces, the model is designed only to explain how the two synergies are produced and, as we will see, the smooth rise and fall of activity is enough to do that.
The heat plots in Figure 5 (top left) show the probability density functions produced by MIIND for each population in the network. As shown in section 2.7, the density function describes the likelihood of finding a neuron from the population with a given membrane potential. The top density plot shows the state of the MN-RF population during the period before the action begins. The lower density plot shows the state when the input is maximal. In the lower density plot, there is a higher probability of finding neurons at the threshold (-51mV) indicating that the average firing rate of that population is higher. The population transitions to the top density once again after the Cortical Drive returns to zero. These transitions are also visible in the probability density functions of the other motor neuron populations due to the indirect excitation from Cortical Drive via the Extensor and Flexor Interneuron populations. Therefore, for all motor neuron populations, as with the EMG signals, the average firing rate output shows an increase to a high level of activity followed by a decrease to rest. In the same manner as the EMG recordings, rank 2 NMF was performed on the time series of average firing rates of the motor neuron populations in the model producing a five value muscle contribution vector and time series activation pattern for both synergies. The afferent feedback inputs to the Extensor Interneuron, InhibRF and InhibST populations were altered and the results were compared to those of the isometric task. Figure 4 shows the results from the NMF process.
For synergy 1 ( Figure 4A), the activation pattern matches the shape of the descending input pattern from the Cortical Drive input (5 seconds of maximal activity with a 1 second ramp up and down). The five muscle contribution values are all well above zero indicating that the activation pattern is a component in the activity of all the motor neuron populations. The synergy 2 ( Figure 4B) activation pattern drops to zero during the contraction and returns to its starting value at the end. The degree to which this pattern contributes to each motor neuron population's activity changes with the level of afferent input to the Extensor Interneuron population. With a large amount of afferent activity, the activation pattern contributes to the antagonist motor neuron populations significantly more than the agonists as was observed in the experimental results. The additional excitation from the Extensor Afferent input causes an imbalance in activity between the Extensor Interneuron population and Flexor Interneuron population. The resultant higher firing rate of the Extensor Interneuron population causes additional excitation of MN-ST and MN-BF and inhibition of MN-RF, MN-VL and MN-VM. Therefore, during the entire contraction task including the resting period, the antagonist motor neuron populations have a consistently higher firing rate than the agonist populations when there is increased afferent input to the Extensor Interneuron population. This is the source of the synergy 2 pattern. As the afferent activity is reduced, the contribution of this pattern lowers until it is almost eliminated across all five populations. When the synergy 1 and 2 time series are weighted according to their respective contribution vectors and summed, MN-ST and MN-BF both show higher activity at rest compared to the agonist populations. Synergy 2 contribution values can therefore be considered an inverse measure of the range in activity of the original output. The similarity of the activation patterns to those from the experiment indicates that the same features are being captured and that the model provides a mechanism for producing them.
In the experimental results, synergy 2 vector values show a small bias towards RF compared to VM and VL, and towards ST compared to BF. As the internal knee angle is increased, the RF bias shows a minor upward trend and the ST bias, moves from a clear difference with BF to almost no difference at 90°. To capture these features, we introduced the populations and connections in the greyed area of Figure 1 to the model to produce the synergy results in Figure 4C. Although it is possible to discern the difference in biases at 0°and 90°, the linearity across the range is questionable due to the lack of data points and variability. At 20°, for example, the contribution to ST/BF appears to be greater than at 0°a nd it is difficult to see if there is a change in RF bias between 60°and 90°. In order for the model to produce these variations, and because it seems unlikely that there would exist additional excitation of RF and ST without some form of a modulation mechanism in the network, we introduced the two inhibitory populations InhibRF and InhibST. The populations have excitatory inputs from the Extensor and Flexor Interneuron populations and also from the external afferent inputs. To reproduce the matching ST bias at 20°, the input to InhibST can be reduced. Likewise, for 60°, the input to InhibRF can be reduced to produce a higher RF contribution. Finally, the biases introduced in the model without InhibRF and InhibST are much larger than observed experimentally therefore a further use of these Inhib populations is to reduce the RF and ST contribution values to match the magnitudes in the experimental synergy across all levels of proprioception. Figure 5 shows the side by side comparison of the synergies extracted from the simulation (including the RF/ST bias) and average synergies from the EMG recordings for position 1. In synergy 2, the activation pattern from the model drops to zero during the contraction whereas the experimental pattern holds at a non-zero level. This is likely due to variability in the EMG time series being represented in synergy 2 by the NMF process whereas the simulation produces no such variability. Comparing the synergy 2 muscle contribution vectors shows the matching trend as internal knee angle is increased. There is no association in the model between angle and afferent input so a selection of four decreasing levels of afferent activity were chosen to illustrate the comparison with the experimental results at 0°, 20°, 60°, and 90°. Across all levels of afferent activity, the model shows no muscle contribution for VL and VM in contrast to the experimental data which show consistent low values. Again, we attribute this to variability in the experimental results where there is none in the simulation.

Discussion
The major findings of this study show that synergy recruitment during an isometric knee extension is affected by proprioceptive feedback, and that these synergies can be reproduced by a neural population model integrating afferent feedback. These synergies are well conserved across conditions and individuals, consisting of simultaneous muscle activation and the balancing of agonist-antagonist recruitment relative to internal angle of the knee. These results are further supported by findings from a novel simulation of the local spinal circuits showing proprioception from muscles contribute to synergy level organization of motor control in humans during isometric tasks. Whilst our results are largely in line with previously established literature regarding synergy recruitment, we differ in our finding that joint angle and therefore afferent feedback regarding joint position alters synergy recruitment (Torres-Oviedo, Macpherson, and Ting 2006;Roh, Rymer, and Beer 2012;Sohn and Ting 2016). We propose that in this case, because the available synergies are constrained due to the nature of an isometric task, adaptation is achieved through changes in recruitment. Inspection of the synergies extracted during NMF analysis matches well with biomechanical interactions expected of the muscles, with synergy 1 reflecting coordinated contraction of all muscle groups and synergy 2 reflecting the agonist-antagonist pairing of the hamstrings and quadriceps. The similarity of each synergy's activation pattern supports that these are the same synergies being recruited at each angle. The changes observed in synergy 2 at different internal angles of the knee demonstrates a clear effect of proprioception on synergy recruitment. Quadriceps and hamstring activity in synergy 2 is balanced at a smaller angle when the contralateral hip is flexed. This further reinforces that these The probability density function of the MN-RF population in the model before input from cortical drive (upper) and during the contraction (lower). Colour brightness indicates the probability of a neuron in the population having that membrane potential. The y axis of the plots represents an arbitrary value for simple exponential integrate and fire neurons. A higher probability at the threshold of -51mV indicates a higher average firing rate for the population. Top Right: The average firing rate of each motor neuron population in the model and the rectified and smoothed EMG output for each muscle. Bottom: The muscle contribution vector and activation pattern time series for synergy 1 (left) and 2 (right) for the model (hatched bars) and for experimental data from position 1 (solid bars with std errors). Indicative levels of afferent activity in the model have been used for comparison with each internal knee angle of the experiment.
changes reflect the influence of afferent feedback due to change in relative muscle stretch requiring re-balancing of the agonist-antagonist activity. This finding challenges long held assumptions that proprioception does not play a role in isometric tasks and synergy recruitment.
Our results contrast with findings that muscle synergies present during isometric force generation in the hand are insensitive to changes in position or load (Roh, Rymer, and Beer 2012). The difference in findings may be reflective of either a different mode of control for muscles of the upper limb or a difference in requirement and outcome of the task resulting in recruitment of more stable synergies. Further examination is required to investigate the cause of these differences. More striking still is the comparison to synergy recruitment during postural perturbations (Torres-Oviedo, Macpherson, and Ting 2006); these dynamic movements required five synergies to adequately explain variance and yet none of these synergies were significantly altered with changes in posture or initial limb configuration. It would seem reasonable to assume that these conditions would be more sensitive to changes in afferent drive due to the greater complexity of the movements involved. We hypothesised that a fixed isometric task should isolate changes in muscle synergies to those due to proprioceptive feedback. Greater restriction of movement may necessitate alteration of synergy recruitment to produce changes in end-point forces, whereas more dynamic movements might allow subtle postural changes to produce the same outcome. To support our experimental findings, we identified a candidate network which, when modulated by afferent feedback reproducing the experimentally observed synergies. The mechanism of action is hard to determine in a human spinal cord due to the lack of direct recordings from these circuits, but using MIIND, we built a population network model based on biological evidence to propose the likely mechanism of action. This simulation reproduced the same synergies as derived from the EMG data, changing the afferent input to the neural populations in this model, the synergy recruitment trends matched those observed with increasing internal knee angle. This model then demonstrates how muscle synergies can be encoded in neural popula-tion circuits and furthermore that synergy analysis of experimental data can be used to directly drive model development.
This model's similarity to previously studied spinal circuits (Hultborn et al. 1987;Shevtsova et al. 2016;Pierrot-Deseilligny and Burke 2005), supports conclusions, (Saltiel et al. 2001;Dominici et al. 2011) that synergy encoding takes place in the spinal cord. Supraspinal input was provided to the MN-RF motor neuron population and Extensor and Flexor Interneuron populations only and so specific bias e.g. in the antagonist muscles, was introduced through the circuitry of the model itself leading to the patterns observed in synergy 2. During this study, it became clear that altering the circuitry of the model (and the level of afferent input) mainly affected the activity of each motor neuron population at rest and its activity during maximal contraction. The contribution value of synergy 2 in the simulation was found to be inversely proportional to the difference between the maximal and minimal activity and this can also be observed in the experimental results.
We have thus demonstrated that a model based on well understood afferent inputs to spinal circuits (Pratt and Jordan 1987;Jankowska et al. 1967) that fits into existing CPG models (McCrea and Rybak 2008) of locomotion in cats, can be effectively used to predict synergies in simple isometric tasks in humans. This speaks to the robustness of the mechanisms in the model. In particular, it is clear that synergy 2 relies on the reciprocal inhibition between extensor and flexor populations which is also an important feature required for locomotion. What remains unanswered is firstly whether the cortical drive presented in this network bypasses the rhythm generation and pattern formation layers of the CPG, and if the extensor and flexor populations can be co-active when mediating a common drive to all muscles for voluntary tasks. Secondly, whether the full CPG model can be applied to understanding other cyclical limb motion like locomoction in humans.
Instead of using the traditional technique of direct simulation of individual neurons, we have instead demonstrated the use of the MIIND simulation package, an environment allowing easy simulation of populations of neurons. It requires only the definition of connectivity at the population level, making it easy to setup and adjust a population network during development. Parameter tweaking is an inevitable part of the modelling process requiring cycles of adjustment followed by simulation. Reducing the need for adjustments to the neuron model itself was one reason why we used the simple exponential integrate and fire instead of a more complex Hodgkin-Huxley style neuron. We were able to reproduce the desired synergy patterns without the need for such complexity. While building the network model we experimented with different connection configurations between populations. MIIND's XML style code, used to describe the network, made it simple to add, remove or adjust connections, as well as to add further populations for the RF and ST bias. For one dimensional neuron models, MIIND can simulate a population network with much greater speed than direct methods and this allowed simulations to be run on a local machine without the need for high performance computing, significantly improving the turnaround time between changing and testing the model. From our experience here, we advocate the use of simple neuron models where appropriate, i.e. reduce the dimensionality of the neural model as far as possible. First, this increases simulation speed and second, this forces thinking about which are the essential neuronal mechanisms before simulation starts.
In conclusion, we have identified a clear effect of proprioceptive drive on the pattern of muscle synergy recruitment during a voluntary task. We have also proposed the likely mechanism of action using a population model which reproduces the same synergy patterns as those observed experimentally, thus pointing to the spinal cord as the site for synergy encoding. Finally, this population network derived from earlier central pattern generator models points towards the fact that spinal circuit components can act in both rhythmic tasks such as walking and in voluntary static tasks.