The proposed model is evaluated by comparing the different properties of the simulated neural activity to those observed during the BMI experiments described in Carmena et al. (2003). The BMI experiments were conducted with macaque monkeys whose goal was to move the cursor to randomly appearing targets and hold the cursor on the target for 150 ms to accept a juice reward. Each experiment included three control modes: (i) pole control during which the monkey controlled the cursor using hand-held pole, (ii) brain control with hand movements (BCWHM), during which the cursor was controlled by the output of the BMI interface while the monkey continued moving the pole and (iii) brain control without hand movements (BCWOHM), during which the cursor was controlled by the BMI interface even though the monkey stopped moving the pole.

Neural activity was recorded from multiple brain area, but mostly from the primary motor area (M1) and the dorsal premotor area (PMd). The BMI interface binned the recorded spike trains in 100 ms bins, to generate the input to a linear filter. The linear BMI filter was trained with data recorded during the last 10 min of pole control, and held fixed during brain control.

The simplifying assumptions on which the proposed OFC model for BMI experiments is based are explicitly stated in Section 2.3.1. As depicted in Figure 1A, the model includes three main parts: (1) Hand and cursor model, (2) Brain model, and (3) BMI filter, as briefly explained in Sections 2.3.2–2.3.5, and detailed in Appendices A and B. Model parameters are detailed in Section 2.3.6, and the effect of process noise in pole control is investigated in Section 2.3.7.

As detailed in Section 2.3.7, we hypothesize that the observed changes in neural modulations following the transition to brain control result from the higher process noise introduced by the imperfect BMI filter. Here we investigate theoretically the effect of increasing process noise during simulations of normal reaching movements on POM defined theoretically in Equation (1), and show that under the assumptions of linear decoding and invariant internal model (see theoretical assumptions 1 and 2 below), the POM indeed increases with the process noise. The covariance matrix of the process noise, Ω^p, is assumed to increase from a nominal low level Ω^p_L to a high level Ω^p_H such that the difference Ω^p_H − Ω^p_L is positive semi-definite (psd) (Horn and Johnson, 2012) and is denoted by Ω^p_H − Ω^p_L ≽ 0 or equivalently by Ω^p_H ≽ Ω^p_L. Note that this analysis is general and does not rely on the specific structures of Ω^p_L and Ω^p_H suggested by Equations (9, 13), respectively.

We make the following simplifying theoretical assumptions (TA):

Within the framework of state estimation and OFC, (TA1) implies that Γ can be expressed as the output of an extended linear system whose state x˜=[xkx^k|k1]T includes both the actual and estimated states (of the hand and cursor). Following the details in Equations (C2–C5 in Supplementary Material), the extended system can be described as: (14)x˜k=A˜kx˜k−1+ξ˜kΓk=H˜kx˜k where ξ˜ is the process noise of the extended system, which captures both the process noise of the actual system and the estimation error (Equation C4 in Supplementary Material). Since both are assumed to be white Gaussian noise, so is ξ˜, and its covariance matrix E[ξ˜iξ˜j]=Ωip˜δij is given by (see Equation C6 in Supplementary Material): (15)Ωkp˜=[Ωkp000Kk + 1HΩkpHTKk + 1T+Kk + 1ΩmKk + 1T0000]

Given (TA2), the Kalman gains K_k are optimized for the nominal system with the low process noise and are not updated to reflect the higher process noise.

Proposition: Let Γ^L_k and Γ^H_k be the outputs of two systems described by Equations (14, 15) with the same parameters except for the process noise that is either Ω^p_L or Ω^p_H, respectively, where Ω^p_H ≽ Ω^p_L. Let POM^L and POM^H be the POM (Equation 1) of the resulting spike counts generated by DSPPs from Γ^L_k and Γ^H_k, respectively, then POM^H ≥ POM^L.

Proof: Proposition C1 in Appendix C (Supplementary Material) proves that under the conditions of the proposition, the cumulative spike-rates Γ^L_k and Γ^H_k satisfy two conditions: (i) E_m[Γ^L_k] = E_m[Γ^H_k] and (ii) E_m[(Γ^H_k)²] ≥ E_m[(Γ^L_k)²], where E_m[ · ] denotes ensemble average over different movements starting from the same distribution of initial conditions. Proposition C2 indicates that those two conditions assure that the POMs of the spike counts generated by DSPPs from Γ^L_k and Γ^H_k, i.e., POM^L and POM^H, respectively, satisfy POM^H ≥ POM^L.

The average cross correlations between the recorded cursor velocity and the recorded neural activity of either PMd (Figure 4A) or M1 units (Figure 4B), demonstrate different temporal relationships. While there is no dominant time shift between the recorded neural activity of PMd neurons and the cursor velocity, the neural activity of M1 neurons clearly precede the velocity (Lebedev et al., 2005). The average cross correlations were computed by first calculating the cross-correlation for each unit from each non-overlapping 2-min interval, then averaging across units to get the ensemble-mean cross correlation, and finally averaging across intervals. Error bars depict the standard deviations of the ensemble-means across the non-overlapping 2-min intervals. Lags at which the cross correlation is significantly lower than the peak at either 0 s in (A) or −0.2 s in (B) are marked with open circles and stars, for significance level of p = 0.05 (standard) or p = 0.005 (Bonferroni corrected for 10 multiple comparisons), respectively (Wilcoxon rank sum test). It is evident that for M1 units, the peak at −0.2 s in Figure 4B is significantly higher than the cross-correlation at zero lag. In contrast, for PMd units, the cross-correlation at negative or positive lags are either significantly lower than the peak at 0 s or at least not significantly different (Figure 4A).

Similar temporal relationships are also observed in the simulations, depending on whether the simulated neurons encode only the estimated state or also the control signal. Figure 4C demonstrates that the activity of neurons that encode only the estimated state do not exhibit any dominant time shift with respect to the velocity, in agreement with the temporal relationship depicted by the recorded PMd units. As detailed in Section 2.3.4, this agreement is in line with the literature on PMd neural encoding (Messier and Kalaska, 2000; Hendrix et al., 2009), so we refer to these simulated neurons as PMd-like neurons.

In contrast, the activity of simulated neurons that encode both the estimated state and the control signal, precedes the simulated cursor velocity, as demonstrated in Figure 4D. In particular, the peak cross correlation at −0.2 s is significantly higher than the cross correlation at zero-lag. This is in agreement with the cross correlation depicted by the recorded M1 units. Given this agreement, and the literature on M1 neural encoding (Evarts, 1968; Kalaska et al., 1989; Todorov, 2000), we refer to these simulated neurons as M1-like neurons. The delay in the cursor velocity can be attributed to the effect of the muscles, which are modeled as a second order over-damped filter, and thus introduce a delay in the forces and eventually the movements that are generated in response to the control signal.

Simulated neural activity generated by the model described in Section 2.3, based on the framework of OFC, successfully depicts the observed changes in neural modulations, as shown in Figure 5. Figures 5A,B depicts the results from a single session to facilitate direct comparison with the experimental results, while Figures 5C,D demonstrates the robustness of the results across a number of sessions with different sets of random tuning weights. As in Zacksenhouse et al. (2007), mean PO^M and PKM were computed by first estimating the PO^M and PKM from non-overlapping 2-min intervals of the binned spike-counts of individual neurons, then averaging across the relevant neurons to get the ensemble mean, and finally averaging across the 10 intervals of each control mode. Error bars in Figures 5A,B depict the standard deviations across the 10 intervals.

Figure 5 demonstrates that the PO^M estimated from a single simulated session increases significantly when switching to brain control, with no matching increase in PKM. Focusing first on the detailed results from a single session in Figures 5A,B, the ensemble-mean PO^M during simulated brain control is significantly higher than of ensemble-mean PO^M during simulated pole control (Wilcoxon rank sum test, p = 0.001 for both BCWHM and BCWOHM), while the ensemble mean PKM did not differ significantly (p > 0.55). This behavior is evident in the simulated neural activity of both PMd-like and M1-like simulated neurons. This reproduces well the significant increase in the ensemble mean PO^M from experimental pole control to BCWHM or BCWOHM (Wilcoxon rank sum test, p < 0.008), with no significant increase in the ensemble mean PKM (p>0.39).

The robustness of the changes in neural modulations following the transition to brain control was investigated by simulating 100 sessions, each with a different set of tuning weights. Pole control was simulated first, while brain control was conducted only in the 78 sessions in which the performance of the BMI filter was within the constraint described above (0.65 < R(v^, v) < 0.85). The scatter plot in Figure 5C demonstrates that in these cases the PO^M in BCWHM and BCWOHM are significantly higher than the PO^M in pole control (Wilcoxon rank sum test, p < 10⁻¹⁰). In contrast, the scatter plot in Figure 5D demonstrates that the PKM in BCWHM and BCWOHM are not significantly higher than the PKM in pole control (p = 0.13 and p = 0.08, respectively).

Our main hypothesis is that the higher PO^M in brain control, without matching increase in PKM, results from higher process noise introduced by the imperfect performance of the BMI filter. As detailed in Section 2.3.7, this hypothesis is investigated here by simulating noisy pole control in which the covariance of the process noise is characterized by Equation (13) instead of Equation (9). The additional term in Equation (13) describes the equivalent noise introduced by the imperfect BMI filter with magnitude that is quantified by α_rec. Figure 3 demonstrates that indeed the cursor movements during noisy pole control (with α_rec = 0.035) are more similar to those during brain control than those during pole control. In particular, the variance of the velocity in the simulated trajectories generated with α_rec = 0.0035 was 88 ± 7 [(cmsec)²], comparable with the variance of the velocity in the simulated trajectories during BCWHM (75 ± 29 [(cmsec)²]), and significantly higher than the variance of the velocity in the simulated trajectories during pole control (Wilcoxon rank sum test, p < 0.0001).

The effects of α_rec on PO^M and PKM are depicted in Figure 6 both when the internal model of the noise is not updated (in line with simplifying assumption SA3 that there is no adaptation) and when it is updated to match the actual process noise. All simulated sessions analyzed in Figure 6 were performed with the same set of randomly selected targets and neural tuning weights. At each α_rec, the average PO^M and PKM and their standard deviations were computed from a singe simulated session as detailed in Section 3.2.

Figure 6 demonstrates that as the process noise increases, due to increasing α_rec, PO^M increases while PKM does not change or even decreases. This is evident by the high positive (>0.99) and negative (< − 0.96) cross correlations between α_rec and either PO^M or PKM, respectively. Furthermore, at all the noisy pole control cases analyzed in Figure 6, the PO^M was significantly higher than the PO^M with the standard process noise (i.e., with α_rec = 0, Wilcoxon rank sum test, with Bonferroni corrected significance level of p = 0.005, corrected for the 10 multiple comparisons), as marked by stars in Figure 6. The effect is similar whether the internal model of the process noise is updated to match the actual process noise or not. The results shown in Figure 6 are averaged PO^M and PKM across all the simulated neurons, but similar trends were also observed when averaging across PMd-like and M1-like simulated neurons separately.

In this and the next section we investigate the sensitivity of PO^M to the baseline covariance matrices of the different sources of noise, i.e., the covariance matrices of the measurement and process noise during pole control. Since the internal model is assumed to be well-adapted for pole control, these are also the covariance matrices of the internal model of the noise. Each sensitivity analysis is conducted with 15 sessions simulated with the same set of neurons but different sets of targets. This is equivalent to simulating longer sessions, which improves the accuracy of the estimated POM. In particular, the mean PO^M across all the simulated neurons was estimated from each session and the average and standard deviations were computed across the 15 sessions.

The Kalman filter depends on the relative size of these covariance matrices, so the value of one of the elements can be kept constant. Furthermore, since we found that the PO^M is insensitive to the ratio between visual position and velocity measurements noise, we kept the diagonal covariance matrix of the visual measurement noise fixed and equals to the one suggested in Todorov (2005).

Focusing first on the sensitivity of PO^M to the baseline proprioceptive measurement noise, we found that the PO^M is insensitive to the proprioceptive position measurement noise. Hence, Figure 7A depicts only the effect of a common scaling factor between the covariance matrices of the visual and proprioceptive measurements noise, i.e., the effect of the proprioceptive measurement noise factor k_m in Table 1. Figure 7A indicates that the PO^M is sensitive to k_m only in BCWHM. When the covariance matrices of the proprioceptive and visual measurement noise are comparable (0.5 < k_m < 2), the PO^M in simulated BCWHM is significantly larger than the PO^M during the experimental BCWHM (Wilcoxon rank sum test, p < 0.004) and even larger than the PO^M in BCWOHM, in contrast with the observed relationship in the BMI experiments (Zacksenhouse et al., 2007). When the proprioceptive measurement noise is 10-fold larger than the visual measurement noise (k_m = 10), the PO^M in simulated BCWHM does not differ significantly from the PO^M during the experimental BCWHM and is lower than the PO^M in BCWOHM, as in the BMI experiments. Further increase in the the proprioceptive measurement noise has little effect on the PO^M even in BCWHM. These observations motivated the selection of k_m = 10 as the default factor in all other simulations (Table 1).

The lack of sensitivity to proprioceptive measurement noise during pole control can be attributed to the perfect internal model and low baseline process noise, which cause state-estimation to be based mainly on the internal model. The high PO^M during BCWHM when the baseline proprioceptive and visual measurements noise are comparable can be attributed to the resulting similar Kalman gains for those two modalities. Due to the imperfect BMI filter, the estimation measurement error in brain control is high, and during BCWHM the proprioceptive measurement may differ from the visual measurement. When both modalities are given similar Kalman gains, the deviation between the proprioceptive and visual measurements contributes to higher variance in the resulting estimated state. Finally, during BCWOHM the proprioceptive measurement noise is increased to infinity to model the lack of relevant proprioceptive feedback in this mode, and hence the PO^M in this mode does not depend on the baseline proprioceptive measurement noise.

During regular pole control, the covariance of the process noise is described by Equation (9) and its magnitude is quantified by α_u, or alternatively by the process noise factor k_u, as described in Section 2.3.6 and Table 1. As noted in Section 2.3.6 and Figure 2, adequate performance of the BMI filter in the simulated BMI sessions can be achieved with k_u = 0.1. Here we investigate the effect of k_u on PO^M.

Figure 7B depicts the changes in the PO^M as a function of the process noise factor k_u. It is apparent that PO^M in brain control is always higher than PO^M in pole control, in agreement with the observed changes in the PO^M during the BMI experiments. However, for k_u = 1, the PO^M in each of control mode is too large compared with the experimentally observed PO^M (Wilcoxon rank sum test, p < 0.002 for each control mode). This further supports the selection of a smaller process noise factor, and in particular k_u = 0.1, made in section 2.3.6 based on the performance of the BMI filter. For this process noise, the relation between the PO^M in BCWHM versus BCWOHM agrees well with the experimental observations.

Figure 7B also shows that PO^M increases with the magnitude of the process noise. This is in agreement with the results described in Section 3.4, though in that section the magnitude of the process noise described by Equation (13) was manipulated by changing of α_rec instead of α_u (via k_u).

The simulations presented so far were based on simplifying assumption SA3, so the internal model did not change due to transition to brain control. Since the monkey stopped moving its hand in BCWOHM, the internal model may change to account for the lack of hand movements. Figure 8 evaluates the effect of changing the parameters of the internal model, and in particular the mass (Figure 8A), or the coefficient of friction (Figure 8B), during pole control while keeping the external model the same. In both cases, the x-axis is the normalized parameter of the internal model (mass or coefficient of friction) relative to the nominal value of the parameter of the external model. The effect is evaluated both with the baseline proprioceptive measurement noise characteristic of pole control and BCWHM, and with infinite proprioceptive measurement noise characteristic of BCWOHM. All simulated sessions analyzed in Figure 8 were performed with the same set of randomly selected targets and neural tuning weights, but not the same set used in other sections. At each parameter value, and for each level of the proprioceptive measurement noise, the average PO^M and PKM in Figure 8 and their standard deviations were computed from a single simulated session as detailed in Section 3.2.

Figure 8 indicates that when the internal model deviates from the external model, the PO^M increases significantly, while the PKM does not change significantly or even decease over most of the range. Parameter values at which the PO^M (generated in simulations with infinite process noise) is significantly higher than the PO^M at the nominal parameter are marked by stars (Wilcoxon rank sum test with Bonferroni corrected significance level p = 0.007, corrected for 7 multiple comparisons). The effect of increasing the covariance matrix of the proprioceptive measurement noise to infinity, as in BCWOHM, is negligible. Hence, changes in the parameters of the internal model may further contribute the observed change in PO^M upon switching to brain control, and especially to BCWOHM.

The implemented OFC model successfully moves the hand and cursor in simulated pole control to reach randomly selected targets. Neural activity is generated to simulate neural recording and perform brain control. The simulated neurons are assumed to encode the signals that are relevant for OFC, i.e., the estimated state and the control signal, and spike-counts are generated as doubly stochastic Poisson processes. The model parameters were taken from the literature or selected to match the experimental observations. In particular, the magnitude of the process noise and the number of neurons were selected so the performance of the BMI filter trained on the simulated neural activity matches the performance of the BMI filter trained on the recorded neural activity. The resulting model successfully generates the main phenomena that it was developed to explain: the POM of the neural activity in brain control is higher than the POM in pole control with no matching increase in PKM.

OFC involves different computational tasks that can be mapped to different brain regions. In particular, it has been suggested (Shadmehr and Krakauer, 2008) that expected costs and rewards are evaluated by the basal ganglia, the forward model is implemented by the cerebellum, and the predicted sensory feedback is combined with the actual sensory feedback at the parietal cortex to estimate the state. Given the estimated state, the pre-motor and primary motor cortices are hypothesized to generate the control signals associated with visual and proprioceptive feedback, respectively. As mentioned above (Introduction and Section 2.3.3) the OFC is simulated in the state space, so all these computational tasks are included regardless of where they are implemented in the brain. In contrast, neural activity is generated to simulate neural recording from M1 and PMd only, since those are the main brain regions recorded in the BMI experiments considered here.

Neural encoding is an open and important research question. Common strategies to investigate this issue involve quantifying the correlation between the neural activity and different relevant signals, or assessing how well those signals can be predicted by the neural activity. However, those strategies are usually restricted to measurable signals. In particular, many studies investigated the correlation between the activity of cortical motor units and different kinematics and kinetics signals (Georgopoulos et al., 1984, 1986, 1992; Kalaska et al., 1989; Ashe and Georgopoulos, 1994; Ashe, 1997; Moran and Schwartz, 1999; Paninski et al., 2004; Kalaska, 2009). Since correlations do not imply causality, the observed correlations with any of the measurable signals do not imply that the neurons indeed encode these signals. Instead, the neurons may encode internal, hidden signals, which are also correlated with the measurable signals (Zacksenhouse et al., 2014) including, for example, muscle activation patterns (Todorov, 2000).

Here we suggest that cortical motor neurons encode signals that are relevant for performing state estimation and control, i.e., the estimated state and the control signal. Since these signals are internal signals generated by the brain, direct correlation cannot be evaluated. Nevertheless, during normal reaching movements, the actual movement should agree well with the estimated state, and thus the neural activity would appear to be correlated with the movement. However, in other conditions, and in particular during brain control, the estimated state may deviate from the measured state. This work suggests that the higher POM during brain control with no matching increase in PKM may arise from such deviations.

We noted that the neural activity recorded from M1 and PMd units exhibited different cross-correlations with the velocity. In particular, the average cross-correlations between the velocity and the recorded neural activity during pole control indicate that there is no time-shift between the neural activity of PMd units and the cursor velocity, while the neural activity of M1 neurons precedes the velocity. These temporal relationships are reproduced well by simulated neurons that encode only the estimated state or also the control signal, respectively. The time shift in the latter case can be ascribed to the delay in the response of the muscles. Hence, the different types of observed cross-correlations, with and without time shift, can be attributed to encoding different signals within the framework OFC, including and excluding the control signal, respectively. Encoding the control signal in M1 but not in PMd neurons agrees well with evidence that the activity of PMd neurons is modulated mainly by the direction and magnitude of movements while the activity of M1 neurons is correlated also with the forces (Evarts, 1968; Kalaska et al., 1989). Further support is provided by previous investigation of the BMI experiments considered here, which indicated that PMd neurons can predict well the position and velocity but not the force, while M1 neurons can predict well the force too (Carmena et al., 2003). While this specific distinction between M1 and PMd neurons may result from the specific sample of units recorded during the BMI experiments considered here, it is well reproduced by the simulated model based on OFC.

Our main hypothesis is that the observed increase in POM following the transition to brain control reflects increasing variance of the encoded signals due to higher process noise. The latter is assumed to occur in brain control due to the limited accuracy of the BMI filter. This hypothesis was proved theoretically and by simulations. Theoretical proof was limited to the case when the estimated process noise in the internal model is not updated, and the neurons encode linear combinations of the estimated state. The simulations demonstrate that increasing process noise during pole control results in higher POM even when neurons encode non-linear functions of the estimated state (e.g., speed and magnitude of the control signal) and regardless of whether the internal model is updated or not.

The proposed model attributes the increase in POM following the transition to brain control to the higher process noise introduced by the imperfect performance of the BMI filter. This suggests that the POM is a proxy for the performance of the BMI filter: High POM indicates that the BMI filter performs poorly while low POM indicates that the BMI filter performs well. As a proxy signal for the performance of the BMI filer, the POM could be used to adapt the BMI filter using reinforcement learning (DiGiovanna et al., 2009; Mahmoudi et al., 2013), especially in cases when an error signal cannot be determined (Bensmaia and Miller, 2014). Furthermore, the proposed model developed here, provides a framework for investigating this and other suggestion for improving BMIs. Finally, the proposed model can be used to investigate other changes that occur following the transition to brain control and in particular the changes in neural tuning (Lebedev et al., 2005).

Our main hypothesis is that the observed changes in neural modulations following the transition to brain control are related to changes in the process and measurement noise induced by the transition, and in particular, the increase in process noise due to the imperfect BMI filter. OFC provides a coherent and concise framework for relating the changes in process and measurement noise to changes in the estimated state and control signals, and finally to changes in neural modulations. However, the key elements are state estimation and the encoding of the estimated state, so other models that include those components may produce similar results. In particular, state estimation and control may not be optimal, but only satisficing (Simon, 1956; Brown et al., 2004). The OFC framework was adopted here since optimization constrains the parameters of the estimation filter and the gains of the control law, resulting in fewer free parameters to be tuned, thus granting the resulting model higher credibility.

An alternative hypothesis that may contribute to the changes in neural modulations involve changes in the internal model. As was shown in Section 3.6, deviations between the internal model and the external system, and in particular deviations in the mass or the coefficient of friction, result in higher PO^M with no matching increase in PKM. Thus, changes in the internal model of the mass or coefficient of friction, may contribute to the increase in PO^M beyond what is attributed to the increase in process noise. Since the transition to BCWOHM involves stopping hand movements, it is conceivable that it is associated with a switch to a different internal model, with lower mass and coefficient of friction. However, the change in the internal model is not justified in the transition to BCWHM, when the monkey continues to make hand movements.

Alternative frameworks for explaining the changes in neural modulations may be provided by intermittent predictive control (Doeringer and Hogan, 1998; Gawthrop et al., 2011) and active inference (Friston et al., 2010, 2011). Intermittent predictive control alleviates the computational load of prediction that is necessary to address time delays (Doeringer and Hogan, 1998; Gawthrop et al., 2011). Due to the imperfect BMI filter, the transition to brain control is expected to result in more frequent updates of the predictor, and hence in higher variability in the predicted state. Since intermittent predictive control is based on OFC, this is an extension of the current model, and will be considered in future work. Indeed it could be argued that due to the inaccurate movements of the cursor during brain control, the monkey is making explorative movements to adapt to the new environment. While the explorative movements may contribute to higher POM, they are expected to also increase the PKM, in contrast with the observation.

In summary, we demonstrated that the observed changes in neural modulations following the transition to brain control can be successfully explained under the assumption that the neurons encode the estimated state and control signal. To be concrete we used the framework of OFC, but similar results are expected even if the controller is not optimal, as long as it relies on state estimation that deteriorates upon switching to brain control.

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.