Supplementary Material: Uniform and Non-uniform Perturbations in Brain-Machine Interface Task Elicit Similar Neural Strategies

Supplementary Fig. 1 shows examples of changes in preferred direction (PD) distribution, measured as the angle differences between the neurons, before (black bars) and after (purple bars) the decorrelation (DeCorr) perturbation. We produced significant shifts in the angles of the rotated pairs, even though only one neuron of each pair was rotated (Fig. 1B). However, we did not observe any significant changes in the overall angles of all the cell pairs (Fig. 1A). The average shift between the baseline and the perturbed PDs for all the neurons were between 4-15°.


Changes in tuning preferences after perturbations
Supplementary Fig. 3 displays weighted PDs contribution (gray arrows) to the population vector, in order to produce movement (red arrow) towards a specific target (dark circle). With each cursor movement, we can update the direction to the desired target (black arrow) to re-estimate possible changes in the neurons PDs. Figure 3A shows arbitrary PDs with the corresponding movement vector, and target direction. It is possible for neurons to take a number of strategies after a perturbation is introduced to the system. Some of the proposed strategies are re-weighting of the neurons contribution to the movement, re-aiming to a latent/virtual target, or re-mapping (Jarosiewicz et al., 2008). We expect to measure a mix of these strategies in our neurons, however in the visuomotor rotation (VMR) trials we consider that the compensation would be lead by re-aiming strategies. Figure 3B shows what we would observe if the PDs did not show any adaptation after the rotation: the cursor movement would shift in the direction of the rotation, and the  updated direction to target would shift in the direction opposite to the rotation. Figure 3C shows changes in the PDs under the re-aiming strategies, where the neurons aim towards a virtual target (dashed light circle) which produced a cursor movement in the direction of actual the target, and a matching direction to target.
Supplementary Table 1 summarizes the total number of sessions ran for each monkey, the number of sessions where a new calibration map was used ("New map"), and the weeks from the beginning of the experiment when these updates occurred, and the number of sessions where the calibration map from the previous session was utilized ("Fixed map"). For each subject, the beginning of the experiments marks "week 1", parenthesis numbers indicate two or more decoder updates in that week (e.g. week 15(2), means two decoder calibrations in week 15).

Preferred directions in baseline trials
We expected the neural units to change their tuning throughout the sessions and across several days, but considered that, similar to the behavior, these neural tunings would stabilize during the baseline unperturbed trials across days. Taking the action preferred directions (aPDs), estimated during active brain control, we Most sessions had larger shifts between the first and second days, in both subjects, and the angle shifts followed a decreasing power curve trend (R 2 ≥ 0.8) for most of these sessions (monkey O: 38/62; monkey M: 13/20).  Figures 4B and 5D show similar trajectories and pooled data during the DeCorr trials, respectively. Rotated neurons are highlighted (magenta traces/bar) from non-rotated pairs (blue/black).

Latent target estimation
In order to estimate possible re-aiming strategies used in either task, we implemented the algorithm proposed by Chase et al. (2010), where a new virtual or latent target direction is estimated, such that it would better explain the changes in the firing activity of the neural units ( Supplementary Fig. 3 C). The method assumes that the neural units firing activity follows a cosine tuning function, it calculates a new set of preferred directions for each neural unit based on the average trajectories for each target direction during online active brain control. After these initial movement-based tunings were calculated, an iterative process was used to estimate the latent target positions which most likely described the firing rates. The estimated latent target direction d * j was computed such that the error between the firing rates measured for the jth target and the firing activity predicted by the new tuning curves was minimized, as stated in 1 (Chase et al., 2010). where y i,j refers to the measured firing activity of the ith neuron in the direction of the jth target, λ i ( d) is the re-estimated movement based tuning curve of the ith neuron evaluated at d, and σ 2 i is the variance of the ith neuron.
After estimating initial latent target directions for each target, we re-estimated the tuning curves using these new directions. Then repeated the latent target and tuning curve estimation steps until the error in the average tuning curves estimation was less than 1%, similar to that used in Chase et al. (2010).
Similar to the preferred direction calculation, we used sets of trials across different days, only using successful trials, we computed these for the baseline, VMR and DeCorr task of both subjects. Supplementary  Fig. 6 shows exemplary locations of the new estimated targets during baseline trials and during both perturbation tasks for both subjects.

Changes in neural population dynamics
Supplementary Fig. 7 shows examples of variations in the average peak cross-correlations values during DeCorr trials for both subjects. We see significant shifts after first introducing the perturbation. However, as we track the peak correlations across sessions, we observe eventual convergence to baseline values in both subjects. Table 2 shows results from one-way analysis of variance test of peak crosscorrelations for both subjects. We tested whether task type (baseline vs DeCorr) had an effect in the peak cross-correlation values across all cell pairs. We found the majority of the sessions had significant shifts in these peak values, and higher cross-correlations were measured during baseline trials, when compared to DeCorr trials.

Factor analysis methods
We used factor analysis (FA) algorithm to reduce the dimensions of the neural space, and estimate latent dimensions that explain the variability of the neural firing activity. We based this analysis using the expectation-maximization (EM) algorithm developed by Rubin and Thayer (1982), and applied to neural events by Yu et al. (2009). This algorithm allows to estimate factors or unobservable dimensions from observable measurements. The EM algorithm applied to FA allows to assume that the observable measurements, the firing activity of neurons, will depend in the activity of unobservable dimensions. Their relationship is described in (2).
where u represents the firing activity of the recorded neurons, z refers to the latent variable, which follow a Normal distribution with zero mean and unitary variance (see (3)). The firing of the measured neural units is conditioned by the latent variables, following a Normal distribution with mean Λz + µ, where Λ  Figure 7. Changes in peak cross-correlations. Average shift from baseline for monkey O (purple circle) and monkey M (black cross) across multiple sessions. Initial perturbation shows the first 4 blocks (128 trials) after DeCorr perturbation was introduced. Then 300-400 trials after perturbation in the same set of neurons, and after 800 trials.

Supplementary
contains the weights that relate the latent dimensions with the neural units in its column space, and Ψ is the covariance, as shown in (2): The EM algorithm allows to estimate the Λ, which contains the manifold that relates the latent variables with the neural activity in its column space. The subspace Λ ∈ R p×q , where p is the number of recorded neurons, and q the number of latent dimensions, can help describe neural activity and task parameters.
This EM algorithm has been implemented for FA by Ghahramani (1996) in MATLAB (Mathworks, Inc.), and more recently applied to neural process by Yu et al. and Cowley et al. (2013). We formatted our neural data as specified in the tutorials for the DataHigh visualization tool, and extracted the estimated intrinsic manifolds from sets of 32 trials. We further use the cross-validation option in the algorithm to estimate the number of dimensions which would maximize the log-likelihood (LL) function of the EM algorithm, and selected the dimensions where, across different days, the LL value had plateaued. Supplementary Fig. 8 shows these LL values across different latent dimensions for the initial set of days. The blue dashed lines highlights the selected dimensions (n = 12).