Optimal feedback control to describe multiple representations of primary motor cortex neurons

Abstract

Primary motor cortex (M1) neurons are tuned in response to several parameters related to motor control, and it was recently reported that M1 is important in feedback control. However, it remains unclear how M1 neurons encode information to control the musculoskeletal system. In this study, we examined the underlying computational mechanisms of M1 based on optimal feedback control (OFC) theory, which is a plausible hypothesis for neuromotor control. We modelled an isometric torque production task that required joint torque to be regulated and maintained at desired levels in a musculoskeletal system physically constrained by muscles, which act by pulling rather than pushing. Then, a feedback controller was computed using an optimisation approach under the constraint. In the presence of neuromotor noise, known as signal-dependent noise, the sensory feedback gain is tuned to an extrinsic motor output, such as the hand force, like a population response of M1 neurons. Moreover, a distribution of the preferred directions (PDs) of M1 neurons can be predicted via feedback gain. Therefore, we suggest that neural activity in M1 is optimised for the musculoskeletal system. Furthermore, if the feedback controller is represented in M1, OFC can describe multiple representations of M1, including not only the distribution of PDs but also the response of the neuronal population.

Appendix

1.1 Prediction of single-neuron and neuronal population responses

Assuming that a single-neuron response is proportional to the muscle control signal, with delay (Sergio et al. 2005; Ajemian et al. 2008), the motor command from the muscle dynamics model could predict neural activity in M1 in an isometric torque production task (Fig. 7a). Although this model is focused on the response of only a single neuron, it could also represent a neuronal population, as a mean activity (Trainin et al. 2007). The mean neuronal population response of muscle-related neurons in M1 can be given by

$$ Q\left( t-{t}_d\right)\propto {\left[\left( f(t)+\lambda \dot{f}(t)\right) \cos \varphi \right]}_{+} $$

(18)

Fig. 7

Prediction of neural responses. The motor command was assumed to be delayed by 100 ms, and the ordinates of the data were normalised in arbitrary units. a Single-neuron response. The top and bottom panels indicate the discharge pattern at the PD and the opposite direction, respectively. b Mean neuronal population response as a function of time and force direction, where the direction is relative to the PD of each neuron in Cartesian coordinates

where t _d , φ and λ are the response-time delay, the force direction relative to the neuron’s PD and the decay parameter, respectively. According to eq. (18) (where t _d  = 0.1 [s] and λ = t ₂), a simple dynamics model can also reproduce the mean neuronal population response observed in monkeys during an isometric force production task (Fig. 7b) (Sergio et al. 2005), similar to a previous study (see Fig. 5 in Trainin et al. 2007).

1.2 Derivation of parametric matrices of the state-space model

This section provides a derivation sequence of the parametric matrices A, B and C in eqs. (6) and (7).

1.2.1 State and input matrices

The right formula in eq. (2), which represents muscle activation dynamics, could be transformed into the following discrete-time form:

$\frac{{\mathbf{a}}_{k+1}-{\mathbf{a}}_k}{\Delta}=\frac{{\left[\left(\mathbf{I}+{\boldsymbol{\upvarepsilon}}_k\right){\mathbf{u}}_k\right]}_{+}-{\mathbf{a}}_k}{t_1}$,

${\mathbf{a}}_{k+1}=\left(1-\frac{\Delta}{t_1}\right){\mathbf{a}}_k+\frac{\Delta}{t_1}{\left[\left(\mathbf{I}+{\boldsymbol{\upvarepsilon}}_k\right){\mathbf{u}}_k\right]}_{+}$,

where ∆ is the single time step of the simulation (i.e. ∆ = 0.01 [s]). Using the state-space vector, the above equation can be represented as

$$ {\mathbf{a}}_{k+1}=\left[\kern1em \begin{array}{cc}{\mathbf{0}}_{6\times 2}\kern1em & \kern1em \left(1-\Delta /{t}_1\right)\cdot {\mathbf{I}}_{6\times 6}\end{array}\kern1em \right]\cdot {\mathbf{x}}_k+\frac{\Delta}{t_1}{\left[\left(\mathbf{I}+{\boldsymbol{\upvarepsilon}}_k\right){\mathbf{u}}_k\right]}_{+} $$

(19)

Equation (4), which describes the dynamics of joint-torque generation, could also be transformed into the following:

$\frac{{\boldsymbol{\uptau}}_{k+1}-{\boldsymbol{\uptau}}_k}{\Delta}=\frac{{\mathbf{J}}_m\cdot {\mathbf{a}}_k-{\boldsymbol{\uptau}}_k}{t_2}$,

${\boldsymbol{\uptau}}_{k+1}=\left(1-\frac{\Delta}{t_2}\right){\boldsymbol{\uptau}}_k+\frac{\Delta}{t_2}{\mathbf{J}}_m\cdot {\mathbf{a}}_k$.

Adopting the state-space vector, the above equation could be represented as

$$ {\boldsymbol{\uptau}}_{k+1}=\left[\kern1em \begin{array}{cc}\left(1-\varDelta /{t}_2\right)\cdot {\mathbf{I}}_{2\times 2}\kern.5em & \varDelta /{t}_2\cdot {\mathbf{J}}_m\end{array}\kern1em \right]\cdot {\mathbf{x}}_k+{\mathbf{0}}_{2\times 6}\cdot {\left[\left(\mathbf{I}+{\boldsymbol{\upvarepsilon}}_k\right){\mathbf{u}}_k\right]}_{+} $$

(20)

combining eqs. (19) and (20) into the following:

$$ \left[\begin{array}{c}\hfill {\boldsymbol{\uptau}}_{k+1}\hfill \\ {}\hfill {\mathbf{a}}_{k+1}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill \left(1-\varDelta /{t}_2\right)\cdot {\mathbf{I}}_{2\times 2}\hfill & \hfill \varDelta /{t}_2\cdot {\mathbf{J}}_m\hfill \\ {}\hfill {\mathbf{0}}_{6\times 2}\hfill & \hfill \left(1-\Delta /{t}_1\right)\cdot {\mathbf{I}}_{6\times 6}\hfill \end{array}\right]\cdot {\mathbf{x}}_k+\left[\begin{array}{c}\hfill {\mathbf{0}}_{2\times 6}\hfill \\ {}\hfill \Delta /{t}_1\hfill \end{array}\right]\cdot {\left[\left(\mathbf{I}+{\boldsymbol{\upvarepsilon}}_k\right){\mathbf{u}}_k\right]}_{+} $$

(21)

Thus, comparing eq. (21) with eq. (7), the state matrix A and the input matrix B, respectively, are given by

$$ \mathbf{A}=\left[\begin{array}{cc}\hfill \left(1-\varDelta /{t}_2\right)\cdot {\mathbf{I}}_{2\times 2}\hfill & \hfill \varDelta /{t}_2\cdot {\mathbf{J}}_m\hfill \\ {}\hfill {\mathbf{0}}_{6\times 2}\hfill & \hfill \left(1-\varDelta /{t}_1\right)\cdot {\mathbf{I}}_{6\times 6}\hfill \end{array}\right],\mathbf{B}=\left[\begin{array}{c}\hfill {\mathbf{0}}_{2\times 6}\hfill \\ {}\hfill \varDelta /{t}_1\cdot {\mathbf{I}}_{6\times 6}\hfill \end{array}\right] $$

(22)

1.2.2 Output matrix

A discrete-time formation of eq. (5) is denoted by

${\mathbf{y}}_k=\left[\begin{array}{c}{\boldsymbol{\uptau}}_k\\ {}{\mathbf{F}}_k\end{array}\right]+{\mathbf{v}}_k$,

where the hand-force is given by F _k  = (J ^T)⁻¹ τ _k . J∈ R ^2×2 is the Jacobian matrix, described as

$\mathbf{J}=\left[\begin{array}{cc}\hfill -{l}_1 \sin {\theta}_2-{l}_2 \sin \left({\theta}_1+{\theta}_2\right)\hfill & \hfill -{l}_2 \sin \left({\theta}_1+{\theta}_2\right)\hfill \\ {}\hfill {l}_1 \cos {\theta}_1+{l}_2 \cos \left({\theta}_1+{\theta}_2\right)\hfill & \hfill {l}_2 \cos \left({\theta}_1+{\theta}_2\right)\hfill \end{array}\right]$.

Thus, the joint torque and hand-force are, respectively, described by the state-space vector as

$$ {\boldsymbol{\uptau}}_k=\left[\begin{array}{cc}{\mathbf{I}}_{2\times 2}& {\mathbf{0}}_{2\times 6}\end{array}\kern1em \right]\cdot {\mathbf{x}}_k\kern0.2em \mathrm{and}\kern0.2em {\mathbf{F}}_k=\left[\begin{array}{cc}{\left({\mathbf{J}}^T\right)}^{-1}& {\mathbf{0}}_{2\times 6}\end{array}\kern1em \right]\cdot {\mathbf{x}}_k $$

(23)

Combining the above equations:

$$ \left[\begin{array}{c}{\boldsymbol{\uptau}}_k\\ {}{\mathbf{F}}_k\end{array}\right]=\left[\begin{array}{cc}{\mathbf{I}}_{2\times 2}& {\mathbf{0}}_{2\times 6}\\ {}{\left({\mathbf{J}}^T\right)}^{-1}& {\mathbf{0}}_{2\times 6}\end{array}\right]\cdot {\mathbf{x}}_k $$

(24)

Therefore, comparing eq. (24) with eq. (7), the output matrix C is given by

$$ \mathbf{C}=\left[\begin{array}{cc}\hfill {\mathbf{I}}_{2\times 2}\hfill & \hfill {\mathbf{0}}_{2\times 6}\hfill \\ {}\hfill {\left({\mathbf{J}}^T\right)}^{-1}\hfill & \hfill {\mathbf{0}}_{2\times 6}\hfill \end{array}\right] $$

(25)

1.3 Prediction of specifications of muscle and neuronal PDs

This section shows that our model can predict muscle PDs and distributions of the neural PDs, which are nearly identical to experimental measurements (Kurtzer et al. 2006; Herter et al. 2007) and other models (Hirashima and Nozaki 2012; Lillicrap and Scott 2013).

1.3.1 Neural representation

The control model allowed us to investigate an algorithm controlling the musculoskeletal system. We assumed that the feedback controller imitated the role of M1 neurons in voluntary movement. However, direct comparison with the distribution of neuronal PDs was not possible. Thus, to examine the neural distribution representing the feedback controller, we adopted a simplified neural network model to transform the state vector into motor commands for muscles through a hidden layer, consisting of 2500 neurons (Fig. 8). This hidden layer connects the input and output layers via simple linear functions. The hidden layer receives the estimated state (other than the desired torque) from the input layer, because error corrections during movement are thought to be realised by associations between M1 and the posterior parietal cortex (PPC) (Desmurget et al. 1999). Moreover, PPC is responsible for the formation of motor plans and state estimations during motor tasks, derived from actual sensory feedback (Mulliken et al. 2008). Although the model was limited because of its simplified structure, it allowed us to investigate the relationships between the biomechanical features of the arm model and the distributions of neuronal PDs. Using this model, the neural activities r _k ∈ R ²⁵⁰⁰ and the motor commands u _k could be represented as neural activities of the hidden and output layers, respectively, as

$$ {\mathbf{r}}_k={\mathbf{W}}_1\left[\begin{array}{c}{\hat{\mathbf{x}}}_k\\ {}{\boldsymbol{\uptau}}^d\end{array}\right]\kern0.2em \mathrm{and}\kern0.3em {\mathbf{u}}_k={\mathbf{W}}_2\cdot {\mathbf{r}}_k $$

(26)

Fig. 8

Neural network model. The network transforms the state vector x and the desired torque τ ^d (input layer) into motor commands for muscles u (output layer) using 2500 neurons r (hidden layer). The input and hidden layers are connected by synaptic weights W ₁, such as the corticospinal neurons, and the hidden and output layers are connected by synaptic weights W ₂, such as muscle innervations

where W ₁ ∈ R ^2500×10 and W ₂ ∈ R ^6×2500 are the synaptic weight matrices and were assumed to consist of the corticospinal neurons in M1 and muscle innervations, respectively. Here, we adopted three assumptions, similar to those reported previously (Hirashima and Nozaki 2012): 1) the desired torque τ ^d is input into the corticospinal neurons via the synaptic weights W ₁ from the input layer with cosine-tuned firing rates (Georgopoulos et al. 1982), 2) each corticospinal neuron innervates multiple muscles (McKiernan et al. 1998; Griffin et al. 2008), and 3) the motor commands are constructed with a weighted linear sum of M1 activities (Morrow and Miller 2003). Then, we assumed two types of connection between the corticospinal neurons and muscles (Rathelot and Strick 2009): the cortico-motoneuronal neurons, which have direct connections with motor neurons, and indirect connections, which affect motor neurons through spinal interneurons. Thus, the weights of the muscle innervations from the neurons W ₂ are allowed to take positive and negative values. However, as no data is available for W ₂, similar to a previous study (Hirashima and Nozaki 2012), we assumed that W ₂ was constant and distributed uniformly in the range [−1/2500, 1/2500]. Here, comparing eqs. (11) and (26), W ₂ ·W ₁ = −L. Thus, the synaptic weight W ₁ for a desired torque is given using the pseudoinverse matrix of W ₂ and the optimal feedback gain L by

$$ {\mathbf{W}}_1=-{\mathbf{W}}_2{\left({\mathbf{W}}_2{\mathbf{W}}_2^T\right)}^{-1}\cdot \mathbf{L} $$

(27)

1.3.2 Muscle and neuronal PDs

In total, 100 combinations of the desired torques were aligned uniformly on a circle in the shoulder-elbow joint torque plane, which had a radius of 1 N∙m: i.e. the desired directions were set to 2π(k–1) / 100 (k = 1, 2, …, 100), and the norms of the desired torque were equal to 1. Neuronal PDs were defined as the desired torque directions with the maximal activation of neurons. Then, we computed muscle PD as the weighted average direction of muscle activity by.

${MPD}_i=\sum_{k=1}^{100}{a}_i(k)\cdot {TD}_k/\sum_{k=1}^{100}{a}_i(k)\kern1em \left( i=1,\kern0.5em 2\kern0.5em ,\cdots, 6\right)$,

where MPD _i and TD _k are the muscle PDs of the i ^th muscle and the desired torque direction, respectively, and a _i (k) is the muscle activity when the desired torque is provided.

The muscle PDs were computed as a weighted average direction of muscle activity in joint-torque coordinates: SF: −17.4 ± 18.8°, SX: 148.2 ± 4.8°, BF: 92.2 ± 5.9°, BX: -164.6 ± 6.6°, EF: 109.3 ± 5.6° and EX: −79.4 ± 5.4° (mean ± SD; Fig. 9a, top). The SF and EX muscle activities were diverted from the mechanical PDs to the direction of the fourth quadrant (i.e. the shoulder-extension and elbow-flexion directions). The muscle activities of BF, EF, SX and BX were also diverted, but in the direction of the second quadrant (i.e. the shoulder-flexion and elbow-extension directions). Consequently, the muscles that preferred both directions were consistent with the experimental observations (Kurtzer et al. 2006). Moreover, the results were almost equal to those obtained with previous models (see Fig. 3b in Hirashima and Nozaki 2012, and Fig. 3d in Lillicrap and Scott 2013).

Fig. 9

Muscle PDs and distributions of neuronal PDs. The top and bottom rows show the default and non-biarticular muscle conditions. a Muscle PDs. The arrows represent the weighted average direction of the muscle activity, and the shaded areas indicate SDs of the weighted averages. b Distributions of the synaptic weightings of the input layer. c Polar histograms of neuronal PDs for the joint torque directions

We represented the feedback gain using the neural network model to compare the distributions of neuronal PDs between the model and experimental measurements. In the neural network model, synaptic weightings of the input layer were distributed with bias to the second and fourth quadrants in the joint-torque coordinates (Fig. 8b, top). The distribution of the neuronal PDs was also biased to the same direction (Fig. 9c, top), and a similar pattern has been found previously for M1 neurons in primates (Herter et al. 2007). The result is also essentially the same as in previous reports, corresponding to Fig. 3g in Hirashima and Nozaki (2012) and Fig. 3b in Lillicrap and Scott (2013).

Here, we performed additional simulations to determine the variable(s) that influenced the muscles, and what caused the bias of the muscle PDs and distributions of neuronal PDs. To evaluate the effects of the musculoskeletal properties, we removed the effects of the biarticular muscles from the moment arms (i.e. set α ₁₅ = α ₁₆ = α ₂₅ = α ₂₆ = 0), and the bias disappeared as in a previous study (see Fig. 5b in Lillicrap and Scott 2013). The muscle PDs were fitted to the mechanical PDs (Fig. 9a, bottom), so that we have SF: 0.0 ± 19.5°, SX: 180.0 ± 5.2°, EF: 90.0 ± 5.2° and EX: −90.0 ± 5.2° (mean ± SD). Then, the synaptic weightings of the input layers were distributed uniformly in a square area (Fig. 8b, bottom). Furthermore, the neuronal PDs accumulated in the second quadrant and the distribution pattern of the neuronal PDs in M1 could not be reproduced without biarticular muscles (Fig. 9c, bottom).

1.4 Codes

We provide the MATLAB (MathWorks, Natick, MA, USA) codes used in this study in ModelDB (https://senselab.med.yale.edu/ModelDB), under accession number 185338. The codes are also available as supplementary material in the online version of this article. To run the codes, it is necessary to install the Multi-Parametric Toolbox (available at http://control.ee.ethz.ch/~mpt/).

The downloaded files are composed of one data set, ‘data.mat,’ and three codes: ‘DesignMPC.m,’ ‘RunSimulation.m,’ and ‘PlotSensoryGain.m.’ DesginMPC.m designs the model predictive controller for the isometric torque production task, and replicates the dataset data.mat. RunSimulator.m loads data.mat, executes the simulation, and produces some of the figures found in Figs. 3- 4 and 9. PlotSensoryGain.m also carries out the simulation and plots the figures found in Figs. 5 and 6.