Active learning of neural population dynamics using two-photon holographic optogenetics

Fig. 2: It is unclear why the GRU is shown if it is not used later.

We include results for the GRU model to justify our use of the linear model. In particular, while similar GRU models have often appeared in the computational neuroscience literature (e.g., Pandarinath et al., Nature Methods, 2018), our results show that on this dataset, we can obtain a more effective fit using a linear model, which validates our use of this model for the active learning experiments.

Line 21: I feel like have read the first sentence many times, almost every year at NeurIPS, in papers on computational neuroscience…

We thank the reviewer for pointing this out and will revise for the final version accordingly.

The exposition in Section 4 is overly technical. It would help if the high level narrative of the paper and the essential steps where explained on a somewhat more intuitive level.

We apologize for any confusion this section may have caused and will add additional intuition and explanation to the final version using the provided extra page.

Intuitively, the chief contention of Section 3 is that the neural population data can be effectively modeled by low-rank dynamics. Section 4 then seeks to determine how measurements should be taken so that, given this low-rank structure, the estimation error of the population dynamics may be minimized as effectively as possible. The results of Section 4.1 provide a quantification of the estimation error of the nuclear norm estimator for such low-rank models. Section 4.2 then proposes an algorithm which chooses measurements that minimize the estimation error, as quantified in Section 4.1. Together, then, these results show how to choose measurements to estimate the neural population dynamics as effectively as possible, assuming they exhibit low-rank structure.

No errorbars in Figure 3 and 4.

Figures 3 and 4 do include error bars denoting 1 standard error, they are simply too small to be visible in most plots, as the variance from trial to trial we observed was quite small. See, for example, the blue curve in the “Best” row of Figure 4(b) for a plot where the error bars are visible.

Line 138: Why did the authors choose the term Causal connectivity matrix? I don’t understand why H would be a causal estimate. It would be good if the authors could elaborate or choose another term.

We use the term “causal connectivity matrix” because that matrix summarizes the causal response of each i-th neuron to photostimulation of each j-th neuron. We use the term “causal” to stress that these relationships are learned from direct causal perturbations to the neural population. By contrast, the vast majority of published work on neural population dynamics involves fitting dynamical models to passively obtained neural data. Due to the lack of causal manipulations in those studies, one cannot distinguish whether statistical relationships arise between neurons due to correlation (e.g., due to a shared upstream influence) versus causation (e.g., neuron i directly influences neuron j) . Such correlative relationships are typically referred to as “functional connectivity,” and we intended for our “causal connectivity” to convey the additional causal interpretability afforded. We will update our paper to clarify this distinction.

Wouldn’t theorem 2 be sufficient or the theorem from ref 67 for the point the authors want to make? How does Theorem 2 reconcile with the claim of the authors in line 176 that this is the first such bound.

Theorem 2 is a lower bound—that is, it shows how much estimation error any estimator must incur, i.e. it is not possible to achieve estimation error less than that given in Theorem 2. Theorem 1 is an upper bound—it bounds the estimation error achieved by the particular estimator we are considering here (the nuclear norm estimator). The lower bound is not novel, but the upper bound we provide is.

For any given estimator, it is important to understand how the measurements taken affect the estimation error for that particular estimator, so that measurements can be taken to minimize the estimation error for that estimator. The lower bound of Theorem 2 does not provide this insight, as it holds for any estimator; Theorem 1 does quantify how the estimation error scales for the nuclear norm estimator, and therefore motivates our choice of sampling.

Section 4.2: I understand that Algorithm 1 somehow falls out of the theory, but it would be helpful for me if the authors could convey the intuition between $\lambda^V$ and $\lambda^{uniform}$, and why one wants to mix them.

By combining Theorem 1 and equation (4.3), our results show that the estimation error scales with a combination of two terms: the first corresponding to the input power in the directions spanned by the low-rank subspace, and the second corresponding to the span of the inputs in all directions. This suggests that, to minimize estimation error, we should devote some amount of input energy to the directions spanned by the low-dimensional subspace, and some amount to cover every direction. This is precisely what Algorithm 1 instantiates: $\lambda^V$ chooses the inputs to target the low-dimensional subspace, while $\lambda^{uniform}$. By mixing these allocations we ensure that we play inputs targeting both relevant objectives.

Line 205: What is meant by “played isotropically”? Algorithm 1: What would be a typical range for l (log T)? Is it possible that the 2^l loop in step 5 becomes infeasibly long?

By “isotropic”, we simply mean the amount played in every possible direction. In Algorithm 1, T can be chosen by the user based on how many samples they can collect, so, for example, the real data we have contains approximately 2000 input-response segments, so this would correspond to $T = 2000$. As such, $2^{\ell}$ will only be as large as the user desires, so it will not become infeasibly large.

Line 248: Is this the same T as in Algorithm 1?

Yes.

The dataset is very simple and does not contain anything to enrich the dynamics of the neural population except for random noise…

The authors could potentially expand their repertoire of datasets to different modalities and more complex behavioral paradigms as this will help them infer whether the algorithm can be effective in more realistic datasets…

I would also add as a limitation the very specific and constrained real-dataset in which the algorithm has been tested on and whether it can have a generally good performance in other, more complex frameworks…

We agree that it would be interesting to further characterize our dynamical models and active learning approaches across a variety of brain areas (e.g., including sensory areas) and while the animal engages in behaviors that may lead to different neural population dynamics. Our collaboration and many others around the world are particularly invested in understanding population dynamics in the motor cortex, however, and hence we have focussed this paper on motor cortical data. Nonetheless, we do hope to extend our techniques to other brain areas in future work. In that vein, we will update the limitations discussed in the paper to assert that further experiments will be required to assess how well our approach generalizes to other brain areas and behavioral paradigms. We found that low-rank autoregressive linear dynamics were sufficiently expressive to model the datasets described in the paper and did not require the added expressivity of the nonlinear dynamical models we fit (see Fig 2, “GRU”). In future work, we hope to repeat these analyses on data collected while the animal engages in a behavioral task. If these data demand nonlinear dynamical models, that could justify developing active learning approaches in those nonlinear models.

Some of the figures are not adequately explained in their legends and some of the ranks are obsolete (Figure 3 and Figure 4)

We thank the reviewer for pointing this out and will add a more detailed description to the legend of each figure for the final version, as well as working to make the figures more concise.

Part of Sector 4.1 seems more like a methods section. The authors could try to make the algorithm explanation a bit more intuitive and justify a bit more the rationale behind each equation.

We thank the reviewer for pointing this out and will add additional intuition and explanation to the final version using the extra page of space. Please see also our response to Reviewer d6y7 for additional explanation of Section 4.

Experimental analysis: Although the paper has interesting experimental results, I would like further explanations of the results. For example, what is causing the discrepancy between the best and worst cases in Figure 4?

In Figure 4, the “Best” and “Worst” plots are the best performing and worst performing train-test split, respectively. It is not entirely clear why our method performs better on some train-test splits than others, as compared with passive approaches, where the train-test split is chosen based on the unique input patterns in the dataset. For some train-test splits there may be a stronger correlation between inputs in the train and test set than for others, or the particular train-test split may make it harder or easier to learn the relevant low-rank structure, both of which could affect the performance of active learning vs passive learning. More investigation into which settings active learning provides a substantial gain is an interesting direction for future work.

We are happy to provide additional explanation on particular points, if the reviewer could bring these to our attention. We will also use the extra page in the camera-ready version to add further explanation to the paper.

Comparisons with other methods: The authors cite multiple papers on actively designing inputs for system identification. Another recent paper that could be relevant is [1]. However, the paper compares the proposed method with too simple baselines, such as random or uniform stimulation. Comparisons with other methods for active design of experiments would be crucial to demonstrate the effectiveness of the method proposed by the authors.

We make several points regarding the chosen baselines. First, for the model class we consider, i.e. AR-k linear models, the approaches we test against are the standard approaches—a fixed design where measurements are chosen to cover all directions is the standard choice for linear settings (see e.g. [Pukelsheim, 2006]), and randomly choosing inputs is also a common benchmark [Simchowitz et al., 2018]. We are not aware of commonly used approaches that differ significantly from these approaches in linear settings. Furthermore, we are not aware of any other approaches for choosing photostimulation patterns in a targeted manner that we could compare against. We hope our work motivates further research into this question, and expands the set of baselines relevant to such settings.

Second, our proposed approach is essentially Fisher Information maximization for low-rank AR-k linear models, while the approach proposed in [1] is a Bayesian instantiation of Fisher Information maximization for a particular class of latent-variable models (see e.g. [Chaloner & Verdinelli, 1995] for justification of this). Therefore, our approach relies on the same principle as that proposed in [1]—Fisher Information maximization—simply for a different model class (low-rank AR-k linear models vs latent variable models). Whether or not the exact approach of [1] is relevant in our setting is therefore primarily a question of the chosen model class. As our results in Section 3 illustrate, the data are effectively fit by a low-rank linear model, and it does not appear that latent variable models (which our GRU model can be seen as an example of) yield a significant improvement. Furthermore, [1] highlights several shortcomings of their proposed approach which would make it difficult to apply in our setting: they state that applying their approach to high-dimensional outputs is difficult (we consider output dimension in the range 500-600), the input design only allows selecting inputs from a discrete set of candidates (while we require continuous inputs), and their approach is not applicable when the state-transition depends on the inputs (which is the case in our setting). Given this, we do not believe that the exact method of [1] is particularly relevant in our setting.

If the reviewer has additional benchmarks that would be appropriate to consider with that we have missed, we would be happy to consider.

Pukelsheim, Friedrich. Optimal design of experiments. Society for Industrial and Applied Mathematics, 2006.

Simchowitz, Max, et al. "Learning without mixing: Towards a sharp analysis of linear system identification." Conference On Learning Theory. PMLR, 2018.

Chaloner, Kathryn, and Isabella Verdinelli. "Bayesian experimental design: A review." Statistical science (1995): 273-304.

Minor comment: The authors could cite more recent uses of data-driven dynamical models for neural data…

We thank the reviewer for noting these relevant references, and we will certainly incorporate them into the final version of the paper. We appreciate the suggestion of using the model from [2], and hope to investigate this model in future work.

How is the lag hyperparameter chosen for the AR models?

We chose this via standard cross-validation procedures: we fit models to the data with several different lag parameters, and then evaluated their performance on held-out data. In practice we found that increasing the lag parameter past 3 or 4 did not substantially improve the model’s fitting ability, so all results in the paper are with a lag parameter of $k = 4$. We will add details on this to the final version of the paper.

In Figures 3 and 4, it's hard to know whether certain estimation error is good enough. Some visualizations of the learned matrices or drawing horizontal lines in the graph to represent the "best" case scenario (e.g., parameters fitted to an extremely large set of input/observation pairs) would be helpful.

Please see the rebuttal pdf for visualizations of the performance for different error measures. In particular, we illustrate the causal connectivity matrix for different levels of the estimation error given in Figure 3, and the neural time-series behavior for different levels of the MSE given in Figure 4.