Active learning of neural population dynamics using two-photon holographic optogenetics
Authors: Andrew Wagenmaker
Lu Mi
Marton Rozsa
Matthew Storm Bull
Karel Svoboda
Kayvon Daie
Matthew D. Golub
Kevin Jamieson
Venue: NeurIPS 2024 poster (2024)
Keywords: active learning
experiment design
neural system identification
neural behavior
Primary area: neuroscience_and_cognitive_science
TL;DR: We develop active learning methods to guide two-photon photostimulation for the purpose of reducing the amount of data needed to estimate an accurate mode of the neural population dynamics.

Abstract: Recent advances in techniques for monitoring and perturbing neural populations have greatly enhanced our ability to study circuits in the brain.  In particular, two-photon holographic optogenetics now enables precise photostimulation of experimenter-specified groups of individual neurons, while simultaneous two-photon calcium imaging enables the measurement of ongoing and induced activity across the neural population. Despite the enormous space of potential photostimulation patterns and the time-consuming nature of photostimulation experiments, very little algorithmic work has been done to determine the most effective photostimulation patterns for identifying the neural population dynamics. Here, we develop methods to efficiently select which neurons to stimulate such that the resulting neural responses will best inform a dynamical model of the neural population activity. Using neural population responses to photostimulation in mouse motor cortex, we demonstrate the efficacy of a low-rank linear dynamical systems model, and develop an active learning procedure which takes advantage of low-rank structure to determine informative photostimulation patterns. We demonstrate our approach on both real and synthetic data, obtaining in some cases as much as a two-fold reduction in the amount of data required to reach a given predictive power. Our active stimulation design method is based on a novel active learning procedure for low-rank regression, which may be of independent interest.

Active learning of neural population dynamics using two-photon holographic optogenetics
Andrew Wagenmaker" University of California, Berkeley
Lu Mi' Georgia Tech
Marton Rozsa Allen Institute for Neural Dynamics
Matthew S. Bull Allen Institute for Brain Science
Karel Svoboda
Allen Institute for Neural Dynamics
Kayvon Daie Allen Institute for Neural Dynamics
Matthew D. Golub University of Washington
Kevin Jamieson University of Washington
Abstract
Recent advances in techniques for monitoring and perturbing neural populations have greatly enhanced our ability to study circuits in the brain. In particular, two-photon holographic optogenetics now enables precise photostimulation encerimenter-specfffee groups of individual neurons, while simultancous two- photon calcium imaging enables the measurement of ongoing and induced activity across the neural population. Despite the enormous space of potential photostimula- tion patterns and the time-consuming nature of photostimulation experiments, very little algorithmic work has been done determine the most effective photostimu- lation patterns for identifying the neural population dynamics. Here, we develop methods efficiently select which neurons to stimulate such that the resulting neural responses will best inform a dynamical model of the neural population activ- ity. Using neural population responses to photostimulation in mouse motor cortex, we demonstrate the efficacy of a low-rank linear dynamical systems model, and develop an active learning procedure which takes advantage of low-rank structure to determine informative photostimulation patterns. We demonstrate our approach on both real and synthetic data, obtaining in some cases as much as two-fold reduction in the amount of data required to reach given predictive power. Our active stimulation design method is based on novel active learning procedure for low-rank regression, which may be of independent interest.
Introduction
Neural population dynamics describe how the activities across a population of neurons evolve over time due local recurrent connectivity and inputs the population from other neurons brain areas. Identifying these population dynamics can provide critical insight into the computations performed by a neural population [1]. Dynamical systems models have enabled neuroscientists to generate and test multitude of hypotheses about how specific neural populations support the neural computations that underlie, for example, motor control [2-4], motor timing [5, 6], decision making [7-10], working memory [1 111, social behavior [12], and learning [13-16].
Eqqually contributing first authors.
Eqqulll contributing senior authors.
Correspondence to: ajjagen@berkeley. edu or 1mi70gatech..n
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

The traditional approach to data-driven modeling of a neural population typically involves two sep- arate stages. First, neural population activity recorded while an animal performs a task of interest. Then, dynamical systems model is fit to the recorded neural responses [17-33]. This approach suffers from two key limitations. First, any inferred structure is purely correlational, and cannot be interpreted with any notion of causality. Second, the experimenter has limited control over how the neural population dynamics are sampled, which can lead to inefficient data cetiectioon-veeramppiinn in some parts of neural activity space while altogether missing others. Given constraints on time and resources in neurophysiologica experiments, there is a strong need for techniques that minimize the amount of experimental data required to identify the neural population dynamics.
We seek to overcome these limitations by actively designing the causal circuit perturbations that will be most informative to learning a dynamical model of the neural population response. For circuit perturbations, we employ two-photon holographic photostimulation (Figure 1), which provides temporally precise, cellular-resolution optogenetic control over the activity of ensembles of neurons [3441]] When paired with two-photon calcium imaging, photostimulation protocols can provide insight into network connectivity by enabling the measurement of the causal influence that each perturbed neuron exerts on all other recorded neurons [36, 39, 42-46]. This platform enables targeted excitation of the neural population dynamics, thus providing the experimenter with unprecedented control over the data collected for informing a model of the neural population dynamics.
Here, we develop active learning techniques for designing photostimulation patterns that allow for efficient estimation of low-rank neural population dynamics and the underlying network connectiv- ity. First, we introduce a low-rank autoregressive model that captures www-dimesiiona structure in neural population dynamies and allows inference of the causal interactions between recorded neurons. We then propose an active learning procedure which chooses photostimulations to target this low-dimensional structure, and demonstrate it in two settings: estimating the underlying causal interactions when using the learned autoregressive model adaptively selecting which samples to observe from our dataset of neural population activity recorded via two-photon calcium imaging of mouse motor cortex in response to two-photon holographic photostimulation. In both cases, we show that our active approach obtains substantially more accurate estimates with fewer measurements compared to passive baselines. Our methodology is based novel analysis of nuclear-nore noreaa--ormmaa regression with non-isotropic inputs. To the best of our knowledge, this is the first approach to demonstrate significant gains applying active learning to low-rank matrix estimation problems, and thus we believe this may be of independent interest.
simulator of the true dynamics, and
2 Related Work
Modeling Neural Responses to Stimulation. Many studies have applied direct electrical or optical stimulation to neural populations to probe the dynamical properties of neural circuits and their relation to circuit function [4, 10, 26, 47-49]. However, these stimulation techniques lack the spatial specificity needed to precisely probe the causal influence of individuals neurons on the population dynamics, and these experimental designs were passive in that the stimulation protocols were specified prior to data collection (with [44, 48] as notable exceptions). Other work has explored a related but separate problem of minimizing off-target effects when photostimulating individual neurons [41].
Low-Rank Matrix Recovery. Low-rank matrix recovery has been intensively researched over the last decade and half [50-52]. However, existing analyses rely critically on the assumption that the set of measurements taken are highly symmetric and satisfy some notion of the restricted isometry property (RIP) or incoherence. The matrix recovery problem of our setting departs from the classical literature in several ways. First, the set of feasible mcasurements we can take is constrained by the physical limits of the photostimulation system. Second, matrix coefficients, we should expect that our resulting set of measurements should be highly skewed by design. Motivated by this, we develop, to the best of our knowledge, the first bounds on low-rank estimation using the nuclear norm heuristic that gives quantification of the estimation error in terms of the precise individual measurements taken (i.e., in contrast more global property like RIP).
aim to adapt and actively learn these
Active Learning and Low-Rank Estimation. The active learning literature is vast, and full survey is beyond the scope of this work. We focus in particular on active learning for dynamical systems, and problems with low-rank structure. The estimation of dynamical systems-the the system identification problem- is central to many areas of engineering and science [53]. The problem of actively designing

---

W. two-photon holographic stimulation
in-sivo neunal population dynamics recording stimulation inputs
10 singular values neural activity data
lnser
imaging laser
nouron Indor
35
singular values of shuffled neural activity
singular values of neural artivity
singuler valure
neural activity recordings
geuron
singular value Figure 1: (a) Two-photon imaging and holographic photostimulation platform (left) and representative image frame (right). Purple circles indicate neurons photostimulated immediately before frame acquisition. Red and blue indicate increases and decreases of firing activity, respectively, relative to before photostimulation. (b) Example time series photostimulation inputs (top) and neural responses (bottom) from 100 randomly selected ncurons (out of d 663 recorded neurons identified in the FoV). (c) Neural responses y/ occupy low-dimensional subspace. Singular values from representative dataset's demeaned neural activity data matrix (blue) indicate substantially more data variance residing
time (3) 10
10g
few dozen dimensions (out of full d = 663 dimensional neural activity space) than is expected by chance (orange, singular values when removing low- low-dimensional structure by shuffling time indices independently for each neuron; note clipped horizontal axis).
inputs to effectively estimate the parameters of a dynamical system has been studied extensively for decades [54-61]. More recently, a variety of provably efficient approaches have been developed for both linear [62, 63] and nonlinear [64, 65] systems. Other related work has considered active learning for latent variable models [66], which are often effective models of neural dynamics. As compared to these works, a key feature of our setting is the low-rank structure present in the data, which to our knowledge has not been previously studied within the active system identification literature.
Beyond dynamical systems, some attention has been devoted to active learning with low-rank struc- ture, in particular works on low-rank bandits [67-70]. While the setting considered in these works is somewhat different- -heen aim to solve a bandit problem, while we are interested in regression they similarly seek to develop active learning approaches which make efficient use of low-rank structure. Also related is the work of [71], which shows that in the related sparse estimation setting, there does not exist more than a logarithmic gain to being adaptive. The results of this work are minimax, however-onle applying certain "hard" problems-and and do not address the matrix recovery problem.
3 Preliminaries
Dataset Details. Neural population activity was recorded in mouse motor cortex using two-photon calcium imaging 20Hz of Immx Imm field of view (FoV) containing 500-700 neurons. Each recording spanned approximately 25 minutes and 2000 photostimulation trials. In each trial, 150ms photostimulus was delivered and was followed by a 600ms response period before the next trial began. Each photostimulus targeted group of 10-20 randomly selected neurons, and total of 100 unique photostimulation groups were defined for each experiment (= 20 trials per group). We evaluate our techniques on four such datasets.
3.1 Fitting Low-Rank Dynamical Models
We first seek to develop effective dynamical models of the neural activity in our photostimulation datasets. Obtaining such models will provide insight into which photostimuli are most informative, and gives us a means to evaluate the effectiveness of our active learning methods. We consider three classes of models: autoregressive (AR) models, low-rank AR models, and nonlinear RNN models. Results from fitting these models are shown in Figure We describe the model details next.
At discrete time t E N, we denote the true neural activity across the d imaged neurons as E Rd the noisy, measured activity as Yt E Rd. and the photostimulus intensity applied across those same d neurons as u, E Rdd Applying stimulus u, influences the measured neural activity at the next timestep y++1- However, just the snapshot Yt may not capture the full true state of the neural population. which may include not just the current neural activity, but potentially also multiple orders of temporal derivatives. To capture these effects, we consider an AR-k model defined as:
Ft-1 -.- (A,Tr-a + B,4+-s) + v, Yt I4 + w. with w. ~ ~NO, ~N))
(3.1)
3

---

direct stimulation to neuron
stimulation inputs
indirect stimulation
AR-K model (full rank] 663
measured activity
gated recureent network (GRU)
imaged activity nouron
prediction
(AF//))
pul aod
AR-K model (low rank)
AR-K KR-K model flow rank)
No
(8)
time (8)
100
0.85
ROC indiract responses
ROC responses
responsess
AUROC
1.0
1.0
true pasiiiv rate (TPR)
positive rate (TPR)
ndirect responede
0.65
663
d
7.5e-4
AR- rank
AR-K rank
MSE
AR-K rank
AR- rank
0.0
0.0
GRU
1.0
0.0
false posilive (FPR) 1.0
5.5e-4
663
100
50
10 GRU
false posilive rate (FPR)
AR-K model rank dim
Figure 2: Example data and cross-validated model predictions. (a) Roll-out predictions of the activity of an example neuron using low-rank AR-k models (A: 4) and GRU networks for 22 example data segments (3.3s per segment; segments separated by brief horizontal spaces). Each model's predictions are seeded with the first k = 4 timesteps (200ms) activity from d = 663 neurons and then unrolled to predict the activity across all d neurons the next 66 timesteps, given the full 70-timestep sequence of photostimulation all neurons. Most responses of neuron i are tied to "direct" photostimulation of neuron i (pink, first row of panels). Several "indirect responses" to stimulation of other neurons j * i that influence neuron through the population dynamics. To avoid showing all indirect stimuli (to d neurons), only select indirect stimuli shown (green, second row of panels). (b) Receiver operator characteristic (ROC) curve of true-positive rate and for response detection caleulated indirect responses only (left) and all direct and indirect responses (right). (c) Area under ROC curve (AUROC) and (d) mean square error (MSE) for all predictions.
tied
balse-positive
where A, E Rddd and B, E dRd describe the coupling between neurons and stimulus at the time lag of timesteps, O...., 1, and offset RRd accounts for baseline neural activity. Given input-observation pairs {uu,
((1,,,))6) ( J. the coefficients {(A,, B=)== / o} of (3.1) can be fit using least squares. Despite its simplicity, this linear model reproduces the recorded neural activity remarkably well (see "full rank" model of Figure 2).
Neural population dynamics are frequently reported as residing in subspace lower dimension than the total number of recorded neurons [2, 7, 17, 72-77J. The population dynamics in our datasets are consistent with such low-dimensional structure, as indicated by the singular value spectrum Fig- ure 1(c). Inspired by this observation, we introduce set of low-rank dynamical models, where each matrix of [(A., B. )=)==)) } is re-defined as diagonal plus low-rank. Explicitly, we parameterize A, DA, + UA, VÀ, and B, = DB. + Un. Vi where 1) E dxd with Dij = 0 for all * j. U e Rdxr, and V d for predefined rank r. The diagonal matrices account for substantial autocorrelation in each neuron's activity (Da.) and for the reliable response of each neuron to direct photostimu- lation (DB,). ). whereas the low-rank matrices (UVT) confer coupling between neurons. To fit these parameters, we optimize the following objective function with gradient descent over all parameters:
minimize A. -=0...
Σt (Ve+1 - k-1 A.
k-1 ===0
B,ut, u)?.
(3.2)
Figure shows that these low-rank models perform comparably to the full rank versions terms of predictive performance; indeed the rank 35 model appears almost indistinguishable from the full

---

rank model. From a statistical perspective, low-rank models have far fewer degrees of freedom, and hence require less data to fit.
To assess whether more expressive nonlinear models could be advantageous, we also fit a gated recurrent unit (GRU) network model, adapted from [22], as shown in Figure 2. Interestingly, the GRU model did not perform as well as the AR-k models, potentially due to the complexities of hyperparameter tuning. Therefore, we focus on linear models in the analysis that follows. Additional details on model fitting are provided in Appendix B.2.
3.2 The Causal Connectivity Matrix
While we require dynamical models to predict the temporal evolution of the neural population activity, we are also interested in inferring how the activity of one recorded neuron causally influences the activity of the other recorded neurons. To address this need, we define causal connectivity matrix, H E Rgddd to be the mapping such that Hu E IR quantifies the total response (across time) of each neuron to a single-timester photostimulus u. That is, > L Hu, where 21, the neural activities generated by the population dynamics if un = u, u,21 = O, and Te5o = 2ooo is the steady state or resting state of the system subject to no photostimulation. If the dynamics are linear, or more specifically follow I), such a matrix H is guaranteed to exist, and can be formed by simply rolling out (3.1) with the appropriate initializations. While H is not explicitly constrained to be low-rank, it is obtained from a low-rank AR- model, it too will exhibit low-rank structure.
-3,... F3, are
In our experimental paradigm, photostimulation acts as a causal perturbation to the population dynamies, and as such, our statistical framework is able to capture causal interactions, as opposed to merely correlative interactions. This is in contrast to the majority of work on neural population dynamies, which involves fitting dynamical models to passively obtained data. Due to the lack of causal manipulations in these 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": we instead use the term "causal connectivity" to convey the additional causal interpretability afforded in setting.
To fit H. we could first fit (A., B. =- and then use these as plug-in estimates for their true values to compute H. Alternatively, we take a more direct approach inspired by the definition of H itself. By inspecting the raw data of Figure 2(a) and observing the rate at which each stimulated neuron returns to baseline activity, it is clear that the system mixes (i.e., forgets the past) quickly. This suggests that the total response due to input asymptotes after some finite number of timesteps т. Thus, we can apply some photostimulus u E Rd at time = 0 and then measure the total response z = Ei= yt+ where y E IRd is the noisy measurement of the true neural response xt. If we repeat this for many pairs {un, 2n)} In then we can approximate H as
H arg ming 1n 1 H'',H, 112.
In this work we adopt this latter approach. Since we believe to be low rank, this amounts to a low-rank matrix recovery problem with matrix-vector observations. In the next section, we will describe how to adaptively choose (un Jn to estimate H using as few (stimulus, response) pairs as possible. Subsequently in Section 5, we will demonstrate that actively designing inputs to accelerate the learning of H effectively accelerates the learning of the full dynamies as well.
4 Active Learning of Low-Rank Matrices
In the previous section, we saw that estimating the causal connectivity matrix H induced by the neural population dynamics amounts to low-rank matrix recovery, where we apply some photostimulus u E dR and observe the neural population response z Re Hu plus noise. In this section we seek to understand how we should choose the photostimuli to estimate the causal connectivity as quickly as possible. To this end, in Section 4.1 nuclear norm regression estimator, and in Section 4.2 present algorithm motivated by these results which seeks to actively estimate low-rank matrices. These results will directly motivate a procedure for designing photostimulation inputs.
we present novel results characterizing the estimation error the
To demonstrate the generality of our results, Section 4.1 we consider a general matrix regression setting. In particular, let O, dd xd2 xd2 be a rank r (potentially non-square) matrix, Pn. B2 xd2
5

---

some input matrix, and assume scalar observations:
Zn (0+, + Tas Th ~ (N,
(4.1)
where (p.) = tr(@, Pa)) for tr(-) the trace a matrix. Note that the setting considered in Section 3.2 is special case of this observation model with θ, <- H and, for each imput stimulation u, measuring the response of (4.1) to inputs by of the form Pj = eyu for = i....,d. d.
Matrix Notation. We let || ||, || llop, || II- denote the Frobenius, operator, and nuclear norm of a matrix, respectively. i denotes the pseudo-inverse of a matrix. vec(-) denotes the vectorization of a matrix, and mat(-) the inverse of the vectorization. We also let Au denote the simplex the set of distributions over a set U.
4.1 Constrained Nuclear Norm Estimator under Non-Isotropic Measurements
We are interested in understanding how we can effectively take into account the low-rank structure of e., if our goal is estimate + from the observations of (4.1). To this end, we consider the following nuclear-norm constrained least-squares estimator for O,
θ arg min (e) =112 Em= (4fn: ) zn)² for K. : (0 ||0||_ ||0|| [|++11-), (4.2) θεK
where here we let ((+.) E T denote the vector where the nth element is (y O.), and ((O)) + 7 the vector of observations, for 7y the vector with elements The Define θ UEV as the skinny SVD such that U Rd xr, V E Rd and consider the linear projection operators Pi, P1
Rd xd2 + Rd xd2 defined as:
P PM = (1-UU MM
and Pi(M) = -
Pi (M),
for any M Rddx We call PI the projection onto the tangent space of O. Note that the dimension of the range of Pi is equal to just r(d; + d2) - 7² d1d2. We are now ready state our main result on the estimation error of è, for Ô as defined in (4.2).
Theorem 1. Define := () (12 o $ P1) 1/ ll Then with probability at least 26:
| --1lr S 1 V (P ())P 21 2 2 Hlp log didz ( 4)
8
P1
+ d + 2 #);
where here 40 (' = En Pn (pn M) and tr(-) describes the sum of the eigenvalues of the linear
operator (P) T P1) Rdixd22 Rd xds
Theorem 1 provides a precise bound on the estimation error of the nuclear norm estimator under arbitrary inputs (vp.)n To the best of our knowledge, this is the first such characterization of this estimator. This characterization is particularly essential in active learning problems, such as the prob- lem considered here, where it is critical that we understand precisely how the estimation error scales with different inputs, in order to determine which inputs will most effectively reduce the estimation error. As the observation model of Section 3.2 is special case of the setting considered in (4.1) with O, + H. Theorem provides quantification of how quickly we can estimate the causal connectivity matrix given some set of inputs; we expand on the implications of this connection in Section 4.2.
Theorem 1 states that the estimation error of the estimator (4.2) scales (predominantly) with the strength of our inputs f., in the tangent space of θ Indeed, if [w.,-.., Wa,] and [e,,..., Udg] are the left and right singular vectors of the full SVD of O.. and IRdidaxr(ds +dz]-r² is matrix with orthonormal columns vec(w.e. for (i,j) (i $ r) U < r), then
tr ((
tr( T vec()),,)e TT))),
so we see that the estimation error depends only on the scaling of Σ vec((px)vec(p)) T in the space spanned by vec(uuu)) for < r or j < r-the tangent space to O., As an example of how this scales, assume that for n 1,..., N the entries of each pn, are IID ( 1) and N 2 r(d, -d2) -

---

Then j 0,
R rid + I(PO'&P))' N. and [| (t (10))1/111np R
VN. This translates to bound of || III - r(ds +d2) 110((1/6) Critically, we see that this does not scale with the total number of parameters, d1 4,d2, but instead with r(dj d2), which could be much smaller. The following result, due to [78], provides a lower bound on the estimation error of any unbiased estimator, and shows that the rate obtained by Theorem 1 is essentially unimprovable.
Theorem 2 (Corollary 1 of [78)). For unbiased estimator θ I8|l O. 1. 2 tr (r( o- P)))).
4.2 Active Learning for Low-Rank Matrix Estimation
Given the above characterization, we turn now to the active learning problem: how can we best choose our imputs 9n, to speed up estimation error of 9,? For simplicity, rather than the general matrix regression setting of (4.1), we consider here the vector regression case, as this is the setting of interest in learning the causal connectivity. In particular, assume that we play some Un E IRda and observe 2n 2, + Aus for 7/,, (0, Id. ). single vector observation corresponds to observing di observations from (4.1). 1), the responses the matrix inputs yoj eyu, for 1, .,..., d1. Assume that O is rank r and let Vo= == [Fi,.... U/] denote the first r right singular vectors of the full SVD of O. Then we have that:
.
(r((P)''PP))
(d1 r)- tr (16 ΣN +r tr((())()))
Dx==Em un (4.3)
This calculation, combined with Theorem I, shows that the estimation error of O scales with a weight- ing two terms: one quantifying the amount of input energy we put into directions spanned by the top- right singular vectors, and one that quantifies the amount of input energy played isotropically (that is, in all directions). Note, however, that the input energy played in directions Vú is weighted by a factor of d1 - r z d1, much larger weight than the weight of given to the term quantifying the isotropic in- put energy. This suggests that, to minimize the estimation error of θ,, we should focus a large portion of our sampling budget target the directions spanned by the top-r right singular vectors of θ,
This strategy admits a transparent intuition. If O, rank-r and some vector u is orthogonal to the top- right singular vectors of O,, then O, u = O. Thus, if we know what subspace the top-r right singular vectors of O, span, playing u orthogonal to this subspace gives us no additional information about O.; in this case we should instead play aligned with this subspace. This is precisely what the first term in (4.3) quantifies, while the second term reflects the fact that we must also estimate the subspace spanned by the top-r right singular vectors of O., for which playing inputs isotropically is optimal.
In general, as we do not know θ., we do not know Vo, and so cannot directly compute inputs minimiz- ing (4.3). To circumvent this, we consider the following iterative procedure, which alternates between obtaining estimate of O., θ, and then playing the inputs that would minimize the estimation mmniniinimiie (4.3) if θ were the true parameter. We present this procedure in Algorithm
Algorithm Active Estimation of Low-Rank Matrices
1: input: horizon N, feasible inputs U, rank feasible set IC
2: θ₁ < I,D 4- 0
3: for = =,,, 1,,, [log2 N do
4: Let Vo denote the top-r right singular vectors of θ and A(A) = Encuduu solve:
X + mig tr((Po (AA) Aumif + arg min tr tr(A())t
5: For 2 steps, play input Un A + A add observations D
6: Update estimate of θ,: θ441 + arg min E(u,)eD ||z Oul
7: return θ
At every iteration f, Algorithm computes two distributions over inputs: XL, which targets the top-r right singular vectors of our current estimate of θ., and Xumir, which plays inputs isotropically, covering all directions. Rather than playing these distributions according to the precise weighting given (4.3), we instead found it most effective to mix them at an equal rate. As we do not initially know which directions are spanned by the top-r right singular vectors of O.. XL is not guaranteed to target the correct directions, especially in early iterations. Xnni plays inputs in every direction,
7

---

however, and thus, even if Ax is not aligned to the top-r right singular vectors of O., will ensure sufficient energy is still being played in the correct directions to allow for learning. Given this, we increase the weight of playing Xpmit relative to that prescribed by (4.3).
Note that the computation of the optimal inputs is form of A-optimal experiment design [79], which in general can be efficiently solved by, for example, the Frank-Wolf algorithm [80]. Furthermore, efficient procedures for solving nuclear-none norm regression problems exist, allowing us to estimate θc+1 on line 6 efficiently [81]. We remark that Algorithm takes as input r, the rank of O. and K, which requires knowledge of |/+ll.- In general, when these quantities are unknown, they can be chosen via standard cross-vall volos--vliittionn procedures.
We emphasize again that the setting considered here corresponds precisely to the setting considered in Section 3.2 with O. f HI, Un the input stimulation patterns, and 2n the observed neural response to input un As such, if the causal connectivity H is low rank, Algorithm and the preceding results provide methodology select input stimuli most efficiently estimate H. In the following section, we will apply this to our photostimulation datasets.
5 Active Learning for Estimating Neural Population Dynamies
We return now to the problem of photostimulus design for learning neural population dynamics, and seek to apply the insights of Section 4 this setting. We present two sets of experiments. In Section 5.1 we use real data to fit a model of the population dynamics, treat this fitted model as simulator for the true dynamics, and then demonstrate that we can learn the causal connectivity matrix И of this simulator faster using active inputs versus passive inputs. Then, in Section 5.2 we split our real data into 750ms long trials of (stimulus, response) pairs (see Section and demonstrate that our active learning algorithm is able to improve the performance of learning dynamical models on real data by adaptively selecting which trials observe, training model on the observed trials, and evaluating on a hold-out set of unseen trials. Here we find that our approach is able to learn an accurate model of the dynamics more quickly than adaptive approaches.
5.1 Active Learning on Data-Driven Neural Population Dynamics Simulator
In Section 3.1, we demonstrated that photostimulation data can be effectively reconstructed using an AR-k dynamics model. Given the effectiveness of these models at fitting our data, in this section treat them as a simulated representation of our true dynamics, allowing us to query them arbitrarily as stand-in for the ground truth dynamics, and seek to determine whether carefully choosing the photostimulation pattern allows for efficient estimation of the causal connectivity matrix H.
Experiment Details. To obtain models of the population dynamics to use for simulation, we fit an AR-A model to each dataset as described in Scction 3.1. In all cases we use an AR-k model with order A = 4. We do one run of the experiments using low-rank model parameters UV with rank r r= 15, and then repeat the experiments using r = 35. In each case, we simulate N = 10000 trials, where each trial corresponds to applying a photostimulus and observing the response for = 15 timesteps, simulating our true data generation process. To simulate measurement noise and other trial-to-trial variability in neural responses, we corrupt the observations with Gaussian random noise. Motivated by the empirically observed fast decay of population dynamics in our datasets, we reset the initial state of the simulator each new trial.
In practice, both the magnitude of the stimuli and number of neurons stimulated at each timestep are constrained by the photostimulatio platform. To reflect this limitation in our simulator, we constrain our inputs to lie in [0, 1]. and also impose a sparsity penalty. Precisely, we choose the input set 2/ in Algorithm to be li= := {u E [o, 1]d []u][a 57), for some value y >> 0 (which we set to y = 30). While this does not explicitly constrain inputs to be sparse, it can be efficiently optimized over, and we found in practice that the optimal inputs within this constraint set are in general at least 2y-sparse. As baseline methods, we consider the following:
Random Stimulation: At each trial n, choose y neurons at random, and set corresponding elements of Un to I.
Uniform Stimulation: Compute unif as in Algorithm I and play inputs un Auilf for all n.
8

---

1.0
Rank = 15
5tmmdan Emna
tatimation Cemmetionn Eanoe
82.4
2005 4003 6205 8000 10100
2000
4005
0905
8000
10000
2303
4500
6200
2000
10500
2305
4000
6000
5000 15000
Number
1.0
0.8
Rank = 33
Baimatoon trea
Aelive (Ouns Random
Eatimaaton Brew
Edimatonn Breee
tatmaaton tree
Unifor
0.5
10300
10300
15000
Number Piais
Number Thoi
Number Traii
Number Thoss
(a) Mouse
(b) Mouse
(c) Mouse 3 (FoV A)
(d) Mouse (FoV B)
Figure 3: Performance of active stimulation design on estimating learned dynamics model. For each mouse dataset, we fit low-rank AR-k model as described in Section 3.1 (for ranks of and 35, and k 4). Treating this simulator of the true dynamics, we compare our active stimulation design procedure (Active, Algorithm I) to randomly choosing groups of neurons to excite (Random), and uniformly allocating stimulation across all neurons (Uniform), and plot how effectively each is able to estimate the connectivity of the simulator dynamics. For each figure and method we average over 20 trials, and plot the mean performance with error bars denoting standard error (note that the error bars are barely visible as the standard deviation is very small).
Our goal is to estimate the causal connectivity matrix H induced by our learned dynamics (see Section 3.2). In practice, we are most interested in estimating the off-diagonal elements of H, as these correspond to causal interactions between different neurons. To this end, we consider the error metric [ - for Ā our estimate of H, M a matrix with all entries except its diagonal, |M
which is 0, and @ element-wise multiplication.
Experiment Results. We present results in Figure 3. As can be seen, across all learned simulators and rank levels, our active learning approach yields a non-trivial gain over both baseline approaches. In particular,
Mouse and both datasets for Mouse 3, observe gain of between 1.5-2x over baselines- is, to achieve given estimation error, our approach requires between 1.5-2x fewer samples than baseline methods. This demonstrates the effectiveness of active learning procedure for estimating low-rank matrices-ov method is able to exploit the low-rank structure present in the underlying dynamics speed up estimation, as compared to methods which do not take into account this structure. Furthermore, it shows that on a realistic simulation of neural population dynamics, we can effectively design stimuli speed up the estimation of the dynamics.
5.2 Active Ranking of Real Data Observations
As described in Section each of our datasets consist of roughly 2000 (stimulus, response) trials. In an online photostimulation experiment, we would choose the photostimulus actively for each trial. Here we seek simulate this process using real experimental data, but offline, by choosing the ordering of the trials available in our pre-collected datasets. This serves testbed for active learning procedures: we can more efficiently learn models in this offline setting, that is strong indication that we should also see gains in online experiments. Indeed, those gains may be even greater online because in our offline setting we are severely restricted choosing from only 100 candidate stimulation patterns. Thus, we interpret the results in this section as a lower bound on the performance we might expect online.
To validate this approach, we randomly choose 20 (out of the 100 total) unique photostimulation patterns and aside a test containing all 20 repeated trials of those photostimuli. This creates an 80%/20% train-test split of non-overlapping stimulus patterns. For D Dtrain and D Dest our train and test datasets, respectively, we consider the following query model:
D
2: for trials n = 1, 2. D trainl do
3: Choose input trajectory D Drssin set D D [T]. D train + Dirain' V[T}
4: Estimate model using data in D, and compute prediction MSE of model on D test.

---

0.259
0.138
NSE on 404t0o T18S
Average
0.236
00137
5.15
22580
0.25
0113
2.575
$0.252
0.10
0.135
1500
1030
1000
1000
1100
500
1302
Kumber Titls
Number 77ais
Namber
Number
M85 wett 200
1022500
USE on Ho Meo Thal
M33 on Ho area Thea
0.550
MAg on the sase Theas
1022225
6.347
0.545
2.515
Best
0.2500
0:146
0.540
02277
2:14
0.535
0.575
1003 Number of WWal
15b0
0.144
509
1000
1500
1600
1550
Number Thall
Number of Tfaass
Numberof Tiass
HSE on weguet Thels
0.248
Active (Dun)
MSE on Hetlow Mtae
MSL on Heddout Ttal
Mat en reeeaat Thals
6246
Rando m
0.470
0.555
Worst
8246
0.124
8465
0.50
0.24
B.12
0.555
0.122
0.400
502
1003
2303
302
1003
1500
1200
1305
1503
1500
(a) Mouse
(b) Mouse
(c) Mouse 3 (FoV A)
(d) Mouse (FoV B)
Figure 4: Performance of active learning estimating photostimulation response on held-out trials.
Fach mouse dataset is split into trials corresponding to a stimulus-response pair, and we consider how these trials might be ordered to obtain more effective estimates with fewer training data trials, simulating the active learning process. Our approach (Active) is motivated by the low-rank excitation criteria Algorithm (see Appendix B.4 for more details) and we compare with randomly choosing which trial to observe next (Random). We plot the accuracy of the learned model in predicting neural responses on held-out test trials. We consider 20 different train-test splits (with 20 trials per split), and include plots of average performance across these splits, well splits where Active has the largest and smallest improvement over Random. We plot error bars denoting standard error (note again that the error bars are barely visible as the standard deviation is very small).
We fit a dynamics model to the current set of observed trials, as described in Section 3.1, and use this model to predict the response of the true system on the held-out test inputs, computing the mean-squared error these predictions as our metric. We apply a variant of Algorithm 1, described in more detail in Appendix B.4, and adapted to the query model above. In particular, to apply Algorithm to learning a full dynamical system, we choose our inputs to target the right singular vectors of B, in (3.1). As baseline method, consider the procedure which randomly chooses unobserved segment from Duran at each iteration.
We run the above experiment for 20 different randomly generated train-test splits each dataset, and present our results in Figure 4, providing the results for the average performance over the train-test splits, well as the best- and worst-case splits for active learning performance. As these results illustrate, though active learning does not give a substantial gain in all cases, in many cases is able to give gain of up factor of 2 in the number of samples required over the random bascline, and in the worst case, matches the baseline performance. This further confirms that taking into account low-rank structure when choosing which measurements to take improve estimation rates, and, believe, is a strong indicator that our active learning procedure would speed up estimation of neural population dynamics in online settings.
Discussion
In this work, we have developed a principled approach to active learning of photostimulation imputs for the identification of neural population dynamics and connectivity. We discuss three limitations of our approach, which each suggest potential future directions. First, we have considered active learning of the causal connectivity matrix and minimization of prediction error, both uniformly across all recorded neurons. Future work may focus on more specific scenarios, such as targeting particular dimensions of the neural activity space changes in connectivity due to learning. Second, while found that linear dynamics fit our data remarkably well, this may not always be the case. Does our methodology effec- tively scale nonlinear dynamies? Finally, our real-data experiments were performed offline. Future work may explore running our a algorithm online during closed-loop photostimulation experiments.
10

---

Acknowledgments
This work was supported by NSF DMR award 2308979 to the University of Washington Materials Science Research Center (AW & KJ), the Shanahan Foundation Fellowship (LM & MSB), the Paul G. Allen Foundation (MR, KS, KD MDG), NIH award R00-MH121533 (MDG), NSF CCF award 2007036 (KJ), and NSF CAREER award 2141511 (KJ).
References
[1] Saurabh Vyas, Matthew D Golub, David Sussillo, and Krishna V Shenoy. Computation through neural population dynamics. Anmual Review of Neuroscience, 43:249-275 2020.
[2] Mark M Churchland, John P Cunningham, Matthew T Kaufman, Justin D Foster, Paul Nuyu- jukian, Stephen Ryu, and Krishna V Shenoy. Neural population dynamics during reaching. Nature, 487(7405):51-565 2012.
[3] David Sussillo, Mark M Churchland, Matthew T Kaufman, and Krishna V Shenoy. A neural net- work that finds naturalistic solution for the production of muscle activity. Nature Neuroscience,
18(7): 2015.
[4] Nuo Li, Kayvon Daie, Karel Svoboda, and Shaul Druckmann. Robust neuronal dynamics in premotor cortex during motor planning. Nature, 532(7600):459-465 2016.
[5] Evan D Remington, Devika Narain, Eghbal A Hosseini, and Mehrdad Jazayeri. Flexible sensorimotor computations through rapid reconfiguration of cortical dynamics. Neuron, 98(5): 1019, 2018.
[6] Hidehiko K Inagaki, Susu Chen, Margreet C Ridder, Pankaj Sah, Nuo Li, Zidan Yang, Hana Hasanbegovic, Zhenyu Gao, Charles R Gerfen, and Karel Svohoda. A midbrain-thallmussoom circuit reorganizes cortical dynamics to initiate movement. Cell, 185(6):1065-1085 2022.
[7] Valerio Mante, David Sussillo, Krishna V Shenoy, and William T Newsome. Context-dependent computation by recurrent dynamics in prefrontal cortex. Nature, 503(77744:78-844 2013.
[8] Federico Carnevale, Victor de Lafuente, Ranulfo Romo, Omri Barak, and Néstor Parga. Dynamic control of response criterion in premotor cortex during perceptual detection under temporal uncertainty. Neuron, (6(4))1067-1777 2015.
[9] Krithika Mohan, Ou Zhu, and David J Freedman. Interaction between neuronal encoding and population dynamics during categorization task switching in parietal cortex. Neuron, 109(4): 700-712, 2021.
[10] Arseny Finkelstein, Lorenzo Fontolan, Michael N Economo, Nuo Li, Sandro Romani, and Karel Svoboda. Attractor dynamics gate cortical information flow during madiiiinnmmaiig.... Nature Neuroscience, 24(6):843 850, 2021.
[11] Warasinee Chaisangmongkon, Sruthi Swaminathan, David J Freedman, and Xiao-Jing Wang. Computing by robust transience: how the fronto-parietal network performs sequential, category- based decisions. Neuron, 93(6):
2017.
[12] Aditya Nair, Tomomi Karigo, Bin Yang, Surya Ganguli, Mark J Schnitzer, Scott W Linderman, David J Anderson, and Ann Kennedy. An approximate line attractor in the hypothalamus encodes an aggressive state. Cell, 186(1):178 193, 2023.
[13] Saurabh Vyas, Nir Even-Chen. Sergey D Stavisky, Stephen I Ryu, Paul Nuyujukian, and Krishna V Shenoy. Neural population dynamics underlying motor learning transfer. Neuron, 97 157:1177-1186, 2018
[14] Britton A Sauerbrei, Jian-Zhong Guo, Jeremy D Cohen, Matteo Mischiati, Wendy Guo, Mayank Kabra, Nakul Verma, Brett Mensh, Kristin Branson, and Adam W Hantman. Cortical pattern generation during dexterous movement input-driven. Nature, 577(7790):386339) 2020.
[15] Xulu Sun, Daniel J O' Shea, Matthew D Golub, Eric M Trautmann, Saurabh Vyas, Stephen I Ryu, and Krishna V Shenoy. Cortical preparatory activity indexes learned motor memories. Nature, 602(7896):24--79, 2022.
[16] Emily R Oby, Alan D Degenhart, Erinn M Grigsby, Asma Motiwala, Nicole T McClain, Patrick J Marino, Byron Yu, and Aaron Batista. Dynamical constraints on neural population activity. bioRxiv, pages 2024-01, 2024.
11

---

[17] Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen I Ryu, Krishna V Shenoy, and Maneesh Sahani. Gaussian- prrscesspp factor analysis for low-dimensional single-trial analysis of neural population activity. Journal of Neurophysiology, 102(1):6144-6) 635, 2009.
[18] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu. Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances in Neural
Information Processing Systems, 24, 2011. 1.
[19] Evan W Archer, Urs Koster, Jonathan W Pillow, and Jakob H Macke. Low-dimensional models of neural population activity in sensory cortical circuits. Advances in Neural Information Processing Systems. 27, 2014.
[20] Scott Linderman, Ryan P Adams, and Jonathan W Pillow. Bayesian latent structure discovery from multi-neuron recordings. Advances in Neural Information Processing Systems, 29, 2016.
[21] Yuanjun Gao, Evan W Archer, Liam Paninski, and John P Cunningham. Linear dynamical neural population models through nonlinear embeddings. Advances in Neural Information Processing Systems, 29, 2016.
[22] Chethan Pandarinath, Daniel J O' O'Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew T Kaufman, Stephen I Ryu, Leigh R Hochberg, Jaimie M Henderson, Krishna V Shenoy, Larry F Abbott, and David Sussillo. Inferring single- trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10): 805-815, 2018.
[23] Joshua Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott Linderman. Recurrent switching dynamical systems models for multiple interacting neural populations. Advances in Neural Information Processing Systems, 33: 33:14867-14878, 2020.
[24] Timothy D Kim, Thomas 7. Luo, Jonathan W Pillow, and Carlos D) Brody. Inferring latent dynamics underlying neural population activity via neural differential equations. In International Conference Machine Learning, pages 5551-5561. PMLR, 2021.
[25] Orren Hamiol--Tabour., David M Zoltowski, E Mika Diamanti, Lucas Pinto, David W Tank, Carlos D Brody, and Jonathan W Pillow. Modeling communication and switching nonlinear dynamics in multi-region neural activity. bioRxiv, pages 2022-09, 2022.
[26] Daniel O'Shea, Lea Duncker, Werapong Goo, Xulu Sun, Saurabh Vyas, Eric M Trautmann, Ilka Diester, Charu Ramakrishnan, Karl Deisseroth, Maneesh Sahani, et al. Direct neural perturbations reveal dynamical mechanism for robust computation. bioRxiv, pages 2022-12, 2022.
[27] Mohammad Reza Keshtkaran, Andrew R Sedler, Raeed H Chowdhury, Raghav Tandon, Diya Basrai, Sarah I. Nguyen, Hansem Sohn, Mehrdad Jazayeri, Lee E Miller, and Chethan Pandari- nath. A large-scale neural network training framework for generalized estimation of single-trial population dynamics. Nature Methods, 19(12):152211577 2022.
[28] Adrian Valente, Jonathan W. Pillow, and Srdjan Ostojic. Extracting computational mechanisms from neural data using low-rank RNNs. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Prucessing Systems, 2022.
[29] Arthur Pellegrino, N Alex Cayco Gajic, and Angus Chadwick. Low tensor rank learning of neural dynamics. Advances in Neural Information Processing Systems, 36:1 16711-11-1,72, 2023.
[30] Daniel Durstewitz, Georgia Koppe, and Max Ingo Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24 (11):693-710, 2023.
[31] Hyun Dong Lee, Andrew Warrington, Joshua Glaser, and Scott Linderman. Switching autore- gressive low-rank tensor models. Advances in Neural Information Processing Systems, 36: 57976-58011 2023.
[32] Aniruddh R Galgali, Maneesh Sahani, and Valerio Mante. Residual dynamics resolves recurrent contributions to neural computation. Nature Neuroscience, 26(2):326-338, 2023.
[33] Omid G Sani, Bijan Pesaran, and Maryam M Shanechi. Dissociative and prioritized modeling of behaviorally relevant neural dynamies using recurrent neural networks. Nature Neuroscience, pages 1-13, 2024.
12

---

[34] Adam M Packer, Lloyd E Russell, Henry WP Dalgleish, and Michael Häusser. Simultaneous all-optical manipulation and recording of neural circuit activity with cellular resolution vivo. Nature methods, 12(2):140-146, 2015.
[35]
Zihui Zhang, Lloyd E Russell, Adam M Packer, Oliver M Gauld, and Michael Häusser. Closed- loop all-optical interrogation of neural circuits in vivo. Nature methods, 15(12):0037-1040 2018.
[36] Selmaan N Chettih and Christopher D Harvey. Single-neuron perturbations reveal feature- specific competition in vl. Nature, 537 2019.
[37] James H Marshel, Yoon Seok Kim, Timothy A Machado, Sean Quirin, Brandon Benson, Jonathan Kadmon, Cephra Raja, Adelaida Chibukhchyan, Charu Ramakrishnan, Masatoshi Inoue, et al. Cortical layer-specific critical dynamics triggering perception. Science, 365(6453): eaaw5202, 2019.
[38] Luis Carrillo-Reid, Shuting Han, Weijian Yang, Alejandro Akrouh, and Rafael Yuste. Con- trolling visually guided behavior by holographic recalling of cortical ensembles. Cell, 178(2): 447-457, 2019.
[39] Kayvon Daie, Karel Svoboda, and Shaul Druckmann. Targeted photostimulation uncovers circuit motifs supporting short-term memory. Nature Neuroscience, (4(2):25922633 2021.
[40] Hillel Adesník and Lamiae Abdeladim. Probing neural codes with two-photon holographic optogenetics. Nature Neuroscience, 24(10):136--366, 2021.
[41] Marcus Triplett, Marta Gajowa, Hillel Adesnik, and Liam Paninski. Bayesian target optimisation for high-precision holographic optogenetics. Advances in Neural Information Processing Systems, 36, 2024.
[42] Christopher A Baker, Yishai M Elyada, Andres Parra, and M McLean Bolton. Cellular resolution circuit mapping with temporal-focusee excitation of soma-targeted channelrhodopsin. Elife, 5: e14193, 2016.
[43] Laurence Aitchison, Lloyd Russell, Adam M Packer, Jinyao Yan, Philippe Castonguay, Michael Hausser, and Srinivas C Turaga. Model-based bayesian inference of neural activity and con- nectivity from all-optical interrogation of a neural circuit. Advances in Neural Information Processing Systems, 30, 2017.
[44] Anne Draelos and John Pearson. Online neural connectivity estimation with noisy group testing. Advances in Neural Information Processing Systems, 33:7437-7448, 2020.
[45] Travis A Hage, Alice Bosma-Moody Christopher A Baker, Megan B Kratz, Luke Campagnola, Tim Jarsky, Hongkui Zeng, and Gabe Murphy. Distribution and strength of interlaminar synaptic connectivity in mouse primary visual cortex revealed by two-photon optogenetic stimulation. BioRxiv, pages 2019-12, 2019.
[46] Arseny Finkelstein, Kayvon Daie, Marton Rozsa, Ran Darshan, and Karel Svoboda. Connectiv- ity underlying motor cortex activity during naturalistic goal-directed behavior. bioRxiv, pages 2023-11, 11, 2023.
[47] Mark M Churchland and Krishna Shenoy. Delay of movement caused by disruption of cortical preparatory activity. Journal of neurophysiology, 97(1):348-359 2007.
[48] Nishal Shah, Sasidhar Madugula, Pawel Hottowy, Alexander Sher, Alan Litke, Liam Paninski, and EJ Chichilnisky. Efficient characterization of electrically evoked responses for neural interfaces. Advances in Neural Information Processing Systems, 32, 2019.
[49] Yuxiao Yang, Shaoyu Qiao, Omid Sani, Isaac Sedillo, Breonna Ferrentino, Bijan Pesaran, and Maryam M Shanechi. Modelling and prediction of the dynamic responses of large-scale brain networks during direct electrical stimulation. Nature biomedical engineering, 5(4):
324-345, 2021.
[50] Emmanuel Candes and Benjamin Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):111-11,, 2012.
[51] Roman Vershynin. High-dimensiooal probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
[52] Martin J Wainwright. High-dimensional statistics: non-asymptotic viewpoint, volume 48.
Cambridge university press, 2019.
13

---

[53]
Lennart Ljung. System identification. Springer, 1998.
[54]
Raman Mehra. Optimal input signals for parameter estimation in dynamic systems-suree and new results. IEEE Transactions on Automatic Control, 19(6)753-7668, 1974.
[55] László Gerencsér and Håkan Hjalmarsson. Adaptive input design in system identification. In Proceedings of the 44th IEEE Conference on Decision and Control, pages 4988-4993. IEEE, 2005.
[56] Dimitrios Katselis, Cristian R Rojas, Häkan Hjalmarsson, and Mats Bengtsson. Application- oriented fimite sample experiment design: A semidefinite relaxation approach. IFAC Proceedings Volumes, 45(16):1635-166) 1640, 2012.
[57] Ian Manchester. Input design for system identification via convex relaxation. In 49th IEEE Conference Decision and Control (CDC), pages 2041-2046. IEEE, 2010.
[58] Cristian R Rojas, James S Welsh, Graham C Goodwin, and Arie Feuer. Robust optimal experiment design for system identification. Automatica, 43(6);993-1008, 2007.
[59] Graham Clifford Goodwin and Robert L Payne. Dynamic system identification: experiment design and data analysis. Academic press, 1977.
[60] Kristian Lindq and Hákan Hjalmarsson. Identification for control:
Adaptive input design using convex optimization. In Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No. 01CH37228), volume 5, pages
IEEE, 2001.
[61]
László Gerencsér, Jonas Mårtensson, and Häkan Hjalmarsson. Adaptive input design for arx systems. In 2007 European Control Conference (ECC), pages 5707-5714, IEEE, 2007.
[62] Andrew Wagenmaker and Kevin Jamieson. Active learning for identification of linear dynamical systems. In Conference on Learning Theory, pages 3487-3582. PMLR, 2020.
[63] Andrew J Wagenmaker, Max Simchowitz, and Kevin Jamieson. Task-optimal exploration in linear dynamical systems. In International Conference on Machine Learning, pages 10641- 10652. PMLR, 2021.
[64] Horia Mania, Michael Jordan, and Benjamin Recht. Active learning for nonlinear system identification with guarantees. Mach. Learn. Res., 23:32-1, 2022.
[65] Andrew Wagenmaker, Guanya Shi, and Kevin G Jamieson. Optimal exploration for model-based
rl in nonlinear systems. Advances in Neural Information Processing Systems, 36, 2024.
[66] Aditi Jha, Zoe Ashwood, and Jonathan W Pillow. Active learning for discrete latent variable models. Neural Computation, 36(3):437-474, 2024.
[67] Sumeet Katariya, Branislav Kveton, Csaba Szepesvari, Claire Vernade, and Zheng Wen. Stochas- tic rank-1 bandits. In Artificial Intelligence and Statistics, pages 392-401. PMLR, 2017.
[68] K Kwang-Sung Jun, Rebecca Willett, Stephen Wright, and Robert Nowak. Bilinear bandits with low-rank structure. In International Conference on Machine Learning, pages 3163-3172. PMLR, 2019.
[69] Yangyi Lu, Amirhossein Meisami, and Ambuj Tewari. Low-rank generalized linear bandit problems. In International Conference on Artificial Intelligence and Statistics, pages 468.
PMLR, 2021.
[70] Yue Kang, Cho-Jui Hsieh, and Thomas Chun Man Lee. Efficient frameworks for generalized low-rank matrix bandit problems. Advances in Neural Information Processing Systems, 35: 19971-19983, 2022.
[71] Ery Arits-Castro, Emmanuel Candes, and Mark Davenport. On the fundamental limits of adaptive sensing. IEEE Transactions Information Theory, 59(1):772-481, 2012.
[72] John P Cumningham and Byron M Yu. Dimensionality reduction for large-scale neural record- 11:110001110099 2014.
ings. Nature Neuroscience,
[73] Patrick T Sadtler, Kristin M Quick, Matthew D Golub, Steven M Chase, Stephen I Ryu, Elizabeth C TyberaKabara, Byron M Yu, and Aaron P Batista. Neural constraints on learning. Nature, 512:423-453
2014.
[74] Saul Kato, Harris S Kaplan, Tina Schrödel, Susanne Skora, Theodore H Lindsay, Eviatar Yemini, Shawn Lockery, and Manuel Zimmer. Global brain dynamies embed the motor command sequence of cacnorhabditis elegans. Cell. 163(3):656-669, 2015.
14

---

[75] Matthew D Golub, Patrick T Sadtler, Emily R Oby, Kristin M Quick, Stephen Ryu, Eliza- beth C Tyler-Kabara, Aaron P Batista, Steven M Chase, and Byron M Yu. Learning by neural reassociation. Nature Neuroscience, 21(4):607-616, 2018.
[76] Juan Gallego, Matthew G Perich, Stephanie N Naufel, Christian Ethier, Sara A Solla, and Lee E Miller. Cortical population activity within preserved neural manifold underlies multiple motor behaviors. Nature Communications, 9()):4233, 2018.
[77] Rishidev Chaudhuri, Berk Gerçek, Biraj Pandey, Adrien Peyrache, and Ila Fiete. The intrinsic attractor manifold and population dynamics of canonical cognitive circuit across waking and sleep. Nature Neuroscience, 2219):112-1520, 2019.
[78] Gongguo Tang and Arye Nehorai. Lower bounds on the bean-squared error of low-rank matrix reconstruction. IEEK Transactions Signal Processing, 59 (10444559455555 2011.
[79] Friedrich Pukelsheim. Optimal design of experiments. SIAM, 2006.
[80] Marguerite Frank and Philip Wolfe. An algorithm for quadratic programming. Naval research
logistics quarterly, 3(1-2):55-10, 1956.
[81] Maryam Fazel. Matrix rank minimizution with applications. PhD thesis, PhD thesis, Stanford University, 2002.
[82] Marius Pachitariu, Carsen Stringer, Sylvia Schröder, Mario Dipoppa, L Federico Rossi, Matteo Carandini, and Kenneth D Harris. Suite2p: beyond 10,000 neurons with standard two-photon microscopy. BioRviv, page 061507, 2016.
[83] Diederik P Kingma. Adam: A method for stochastic optimization. 1rXiv:1412.665 2014.
arXiv preprint
15

---

A Proof of Theorem 1
As θ is minimizer, we have that |$ =113 ||
=11 Consequently,
[ 2112 <
min G.tAck:AAl22)
| A) 21 = θ ş (e. + E | 2 p)
|6 θ || < p.
Thus, our strategy will attempt to find minimum p that makes the first expression true. Equivalently, we show that the following quantity is non-positive
(0.)) =1
min +AEK:|AA
| + + 21
AA (A); ((O.) -) (}(()
min
2(T(A), 11) |
M+AEX:| Al
2(t(A), n) (
max (2p M+ |Al|st
max
Intuitively, when is large, the quadratic term will dominate the inner product making the entire expression non- non-positive. We will employ geometric fact about the nuclear norm ball.
Lemma A.1. Ler K = (0: ||||= I [[.1l-]. If + E K then [| S -AA,UVT) <
(A)]ee
Proof. If ||0, + All. > |0,||- then O + A & K. By the convexity of the nuclear norm ball, we have that
||0 + All. > ||-||. . (A. UV (A, w)
VW W P1(W), [[]] s 1.
Consequently, as the dual norm to ||2 is the nuclear norm
||. we have
||, + || + 2 ||0||| 1|+ + (A, UV + (A)]l
Thus. if (A, TVT) + (A) > then θ, A & K. Consequently, if e. + A IC then AA + |P (A)Il. < 0.
Recalling that ||| S ||M for any matrix M, an interesting consequence of the above lemma is that [|P) (A)IIF || < (A) which implies (A) < || < 2| (A) That is, the total error dominated by the error in the tanget space of +..
Applying this lemma, we have that
max ((0(A), - I <
max (
-A ((T(A), n) (
M M+AEK:||A|||=
max
2 2(0(P) + h (A)), m) (()
A:IIP SI PAA)IL -==
max
2( (A)), m) (
+
P1 S1 (1) 2(4(P1 m)
max
where the equality uses the fact that A = P1 + (A) and the linearity of 4. For any v R we have
N (A,
Σ
X,vn) =: (A, (v))
n=1
where Φ* denotes the adjoint of Φ. We also recognize that the operator g dd + Rd xda
also linear, defined as (MM) = Φ (L(X,, M) = Σ (X, M). Consequently ( well-defimed and the same operator as o* after taking the square root its eigenval- ues.
The next three lemmas bound the two terms of above. Combining them yields the result of the theorem.
16

---

Lemma A.2. With probability at least 1 8 we have
max ( n) 2t
+ ddz + log((116))
Proof. Computing this term amounts to bounding a Gaussian width. Begin by recognizing that by the non-expansive property of projections, || ( |F S | S which results in the simplification:
max (()(P)(A)) n) S max PL
())(()(A))
N 2(
max (A))les n=1
2
Tn P1(A))
max ( 2(P
mm PL(A)
==1
S 2|
Σ nn Xn)lop
using the fact that the operator norm and nuclear norm are dual to each other, and that the projection Pi non-expansive. Note that
X, mat Th vec(Xx))) mat vec (xnoXX))))((,,) 22 =: (0'2)
where v/ E IR xd2 with This ~N(0, 1). We then observe that
I Σ Th Xn)
[| (mat( vec vec(Xx)vec(Xx, T ) vec(t) ))
< P| ( lopllll
which completes the proof of the first claim. Recognizing that Il7'llpp is just the largest singular value of a Gaussian matrix, we find that |/sup//ulls [1lllla (wv (uy], I vdi + vdz by Exercise 5.14 of [52]. Applying a sub-Gaussian tail bound completes the proof.
Lemma A.3. Define μ = ( $ ())
)/2 Then 2 ( m) 1
max (MIP
2( (P)() - (1 -1- (A))
max
(A)
Proof. Recognizing that Φ(A) RN we have
( [| + P (3
||((P) + ((P: (A))
= ( + (| + (A)) (PLL()))
To aid in readability, we make number of notational modifications. First, we drop parentheses so that de (A)) is just notated as & A. Second, we define M1 == (A where M is the pseudoin- verse. If th & is invertible restricted to the range of Pp.. then P A = (P) x x t/2 1) 12 /
for all A. Thus,
= (P
-KAAA,
"" (T A)
S 1" P/ || 3) 1/2P A//
< (PP 1/2 lop 1 1) 1/2 1//1////444 A/p |1 3) 1/2 AlF
17

---

where the last two lines follow from Cauchy-Schwarte Observe that for μ < 1 we have
a² + +2 - 2aby = (1 - u1)a² + (1 - -1)b²² + na² + 12b² 2ab
= (1 1-)aa² + (1 -1)b² + uuaa b)²
(1 -))aa
Thus, jf : ((0000)) 1/2 lop then
( |0 + (A(L) (
( 3 + (A) - 0) (0'0) P1 Alp
Z (1 - (
Thus.
max (I) ((0()) - ( -1(A))
2(0(P)) n) (1 1)) (PP)
max SII
P:I)P SIP 11 (PP() (). m7) (1))))((P((A))1
max
(P 1), 17 - 111-pp1 (
max A:IIM 10/2/22//
where we've used the facts that || |1. S 11. ||I and t² = =11 = II + (I)L < (A)lE
Lemma A.4. Let K = min(log, (N), Then for any a 0. if
t I- 1- 16tr (P))))))))) + (P ( log (K)) 2a]|( 4- P1) P
1 1
then with probability at least 1 8 we have
al +
max
2( ()) - (1 - (P)))
Proof. The linear operator PP:: Rdixd2 + R.N can be decomposed as &P En Bn Wn Un where {w,) In are orthonormal on RT. (th).. n are orthonormal linear operators on Rdd xdz and 3a 2 0 are decreasing. For = 0, 1, 1...., min (lo,, d2d2) - 1 let Wt- = [Wye.... W2 so that
834
max |A|
u W. (A)
min u TW &P B2811
| |
Then
2(0(P - (1-
max
A:
K WT (P)) ( (A)
max
2A
k=0
2(W (P) WT. m) - (1 -(1-1) (A))
max
(A)Il 21/V2
W 7ll2 -1 (1 ( ( (PI
max
AIIl
K
(2 - 11)
t 282er++
W 7l12
B2
2
0=0
18

---

where the first inequality holds by Cauchy-Schwartz and the second holds for all 2 W n||22v2 Moreover, ty W 7/12 ((1-H) 32 I -at if
m
(1-1)
1 n|222v2
2a
8 W nl12
2a
max 0
(1 -1- 12
O...,KK K (1 - 11)3
max
max O
(1 122
+ max Am O...., (1 - 11)
Note that for K = min (NV), d1d2) we have
P(URET nll 2 2 + / 1og((K1)))) - PPURE nll 2 Bll nll2] + / log( K/ K/)}} < 6.
On this good event, we have that
IWT ml1
(V2 + / (2log(K/8))
max 0= K
B2
max k.....,
2
4 log(K/0)
k= K 32h+1
B2
N
4 4log( K/6)
- 2
+
max
82 1=
N
$2
==1
= 20r
+ 4 (P(0-()))
log(K/o)
where we use the fact that
N
K 2
K
2
2
32
32
K K-
02
k=D
n= 24
k=U
B Additional Details on Experiments
B.1 Further Details on Dataset
The photostimulation data were collected from transgenic reporter mice Ai229, which express Cre- recombinase-dependent cytosolic GCaMP6m and soma-targeted ChRmine, crossed with the Vglut1 cre mouse line. Imaging and photostimulation experiments were performed on Bergamo (Thorlabs) microscope equipped with a 16x (0.8 NA) Nikon objective. Post-hoc motion correction and neuron segmentation were performed with the Suite2p package [82] https//ggtttbbcommMoossLm//sssst
B.2 Further Details on Experiment of Section 3.1
We split each of our photostimulation datasets into con-ovellapping training and test datasets. All models were trained exclusively using the training dataset and were then evaluated (as shown in Figure 2) using the test dataset. To build our test datasets, we randomly chose (out the 100 total) unique photostimulation patterns and then included all 70-timestep windows about each of the 20 instances of those unique photostimuli. The resulting test set amounted to ~20% of each dataset, In Figure 2, all models were evaluated using these 70-timestep test sequences of the form (M,, ) where Y/ Rd is the recorded neural activity and E Rd is the photostimulation delivered at time t. During evaluation on a given test window, all models were provided {y, ) and {ue 7 {ve 170
predict
Autoregressive-k models: We fit the full-rank AR- models training datasets via linear regression by expressing
h-1
Vt++
Σ
Astle +
B Baut-s +
(B.1)
==0
s=0
19

---

as Y XW, where
1/1
yo
y-1
y1-k
U0
u-
u1
17
1/2
1/1
yo
1/-
01
12
Y
X
:
1
+1++1
L/T
MT-1
Mr+1-K
WT
uT-1
WT+1-k
1
Ao
A1
A
W
Bo
B
Bx-1
and the closed-form solution is W- -X For the low-rank AR-k models, we fit all parameters via gradient descent using Adam [83] over 100 training epochs with a learning rate of 0.01. Gradient descent was implemented in Py PyTorch and ran single NVIDIA Tesla T4 GPU.
During evaluation of the AR-k models, for each test window we first computed 7x++1 given {3/t ue ue using (B.1). Then for all subsequent predictions, on the right-hand side of (B.1) we replaced all instances of yt with De for > In this manner, each entire roll-out prediction of {ye J=kk++1 used all photostimulation inputs {u, ) but only the first & timesteps of neural activity {yr 1E All AR-k models in this paper used k- = 4.
Gated recurrent unit (GRU) networks: GRU networks were loosely based on the sequential variational autoencoders of [22]. Each model consisted of an encoder GRU network that encodes initial timesteps recorded neural activity into a bottlenecked initial state for decoder GRU network. The decoder then unrolls an entire predicted timeseries of recorded neural activity given (as input) all photostimulation that was delivered over that time period. Model fitting proceeded by optimizing the evidence lower bound (ELBO) with respect all model parameters. Both encoder and decoder GRUs had 512 hidden units. We used Adam optimization with a learning rate of 0.001 over 4000 training epochs of batch size 100. Models were implemented with Py PyTorch, and optimized on a single NVIDIA Tesla T4 GPU.
We evaluated all models using roll-out predictions on held-out test windows. We quantified performance with mean squared error between recorded and predicted neural activity for each neuron. We also performed thresholded response detection, whereby detections were defined as timesteps at which given neuron's 's measured calcium fluorescence exceeded predefined threshold. To calculate a receiver operator characteristic (ROC) curve, we enumerated range of thresholds, normalized by the standard deviation of each neuron's 's empirical activity distribution, and performed threshold detection separately on the real and mrdel-predicted neural activity traces. We then compute the overall false-positive rate and true-pose noe-positive rate at each threshold level to trace out an ROC curve. We calculate area under the ROC curve (AUROC) to quantify the accuracy of each model.
Evaluation metries:
Longer roll-out evaluations: To assess AR-k models' ability to predict over longer time horizons, we implemented another train-test strategy, where the first 80% of timesteps in a recording are used for training, and the last 20% timesteps (6736 steps) are used for testing. During the test phase, we use the same procedure described above, providing only the k = 4 initial timesteps of neural activity and then unrolling predictions over the remainder this long test window. We report these results in Figure 5.
20

---

direct stimulation
stimulation inputs
indirect stimulation
ARK model (full rank) r- 663
measured activity prediction
AR- model (low rank) r
imaged sciivity neuron)
(ARFF)
AR-K model (low rank)
AR-K model (low rank
100
25
time (8)
50
100
0.85
ROC curve indirect responses
ROC ourve for responses
AUROC
resporses
troe postive rate (TPR)
true postiu rate (TPR)
ricired restonset
0.55
5.8e-41
MSE
A2-K Bnk
0.0
0.0
false positive rate (FPR)
1.0
0.0
0.0
false positive rate (FPR)
1.0
4.8e-4
663 100 50 35 15 10 5
AR-K model rank dim
Figure 5: Longer roll-out evaluations. Same format as in Figure 2.
B.3 Further Details on Experiment of Section 5.1
To fit the H parameter in this experiment, we generate observations described in Section 5.1. We then estimate H as:
arg min
Σ
|| - Hul
HER
(u,=)ED
for u our input, and 2 = E 2e the observed response, where here 2; are the observations generated from playing input u, and 15.
As K is defined with respect the nuclear norm of the true parameter, which we do not assume is known, we run each method with range of possible values for the nuclear-norm constraint, and plot the performance of each method for the constraint value that has minimum error. We state the value of the nuclear-norm constraint used for each plot below:
Active
Random
Uniform
Mouse rank
10
10
10
Mouse 1, rank 35
10
10
10
Mouse 2, rank 15
5
5
Mouse 2, rank
5
Mouse 3 (FoV A), rank 15
25
25
25
Mouse 3 (FoV A), rank 35
25
25
25
Mouse 3 (FoV B). rank 15
50
100
100
Mouse 3 (FoV B), rank 35
100
100
100
Table 1:
Nuclear-Norm Constraint Settings for Results of Section 5.1
To choose the input rank of Algorithm 1, we ran our experiment with several different ranks and provide results for the best perforreronnnng rank. We found, however, that results are typically robust to the setting of the rank parameter of Algorithm 1, and our choice did not significantly impact performance. Furthermore, believe this could effectively be chosen adaptively. We state chosen values below.
For all experiments, we add observation noise distributed as (V0 to Σt= It.
21

---

Input Rank
Mouse 1, rank 15
10
Mouse 1, rank 35
10
Mouse 2, rank 15
5
Mouse 2, rank 35
5
Mouse 3 (FoV A), rank 15
25
Mouse 3 (Fo' A), rank 35
25
Mouse 3 (FoV B), rank 15
25
Mouse 3 (FoV B), rank 35
50
Table 2: Input Rank r for Results of Section 5.1
(a) Ground Truth
(b) Error 0.1
Error 0.2
(d) Error 0.3
Error 0.4
(f) Error = 0.5
(g) Error = 0.6
(h) Error 0.7
(i) Error 0.8
(j) Error = 0.9
Figure 6: Causal connectivity matrix for Mouse FoV with different levels of estimation error (corresponding to Figure 3).
To ground the estimation error values shown in Figure 3, in Figure 6 will illustrate the causal connectivity matrix for Mouse 3 FoV with different levels of estimation error.
B.4 Further Details on Experiment of Section 5.2
For this experiment, on the data D we observed thus far, we fit the AR-k model described in Section 3.1 with & We found that for this experiment, simply using the least squares estimator with no low-rank penalty produced the best results. We use the same estimation method for both our method and the bascline method.
Given input response trajectory in the test set, (x1....,,ay with input u, to compute the test MSE, provide our learned dynamics model with the initial state æ1 and input u, and then roll this out for 15 timesteps to generate predictions .2...., 115. Precisely, if À and B are our estimated parameters, we let
T2 Ar Bu,
Ad,, t 1.
We the compute the MSE on this segment
15
14
2
not immediately obvious how to apply Algorithm to this setting, since we must choose each trajectory sequentially, and once we have observed trajectory it can no longer chosen again. Rather than solving the optimization of Algorithm I to find the best inputs, instead seek to
It is
22

---

(a) Neuron 0, MSE 0.175 (b) Neuron 0, MSE 0.150 (c) Neuron 0, MSE= 0.140 (d) Neuron 0, MSE = 0.138
(c) Neuron 3, MSE 0.175 (f) Neuron 3, MSE 0.150 (g) Neuron 3, MSE (h) Neuron 3, MSF 0.138
(i) Neuron 95. MSE =() Neuron 95, MSE 0.175
(k) Neuron 95, MSE (1) Neuron 95, MSE
0.150
0.140
0.138
Figure 7: Estimated neural activity vs true neural activity on heldout trials for Mouse 2, Neurons 0, 3, and 95, at different levels of overall MSE on heldout trials (corresponding to Figure 4).
iteratively choose the next input that would maximize "information gain" some sense. In particular, note that applying the Frank-Wolfe algorithm [80] to the objective, if we have inputs U available:
min tr((V) A(A)V)-),
λE
the update is given by:
0 min u T V(V A(A.)V) -V TT well
A++1 n)) + xlffu u++13
for learning rate %.-
In this experiment, we simply choose u. as above, with 2/ the set of remaining active inputs in D Derain and A (A,,) replaced with Σx=1 u.u. This therefore approximates the solution to the experiment design of Algorithm 1, and has the advantage of being very computationally efficient. Furthermore, set V to be the right singular vectors of B. We believe this is reasonable in dynamical system settings with fast decay.
The primary hyperparameter for this experiment is the choice of r, the rank of V. As in the previous section, we not find the results particularly sensitive to setting of r. For each dataset, we ran with E [25, 50, 75, 100, 125, 150], and include results for the best-performing setting.
For both of experiments in Section 5, we ran on 56 Intel(R) Xeon(R) CPU E5-2690 v4 @ 2.60GHz CPUs.
To ground the MSE values shown in Figure 4, in Figure 7 plot the predictions from the estimated model at different MSE values on heldout trials for Mouse 2.
23

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
We validate all claims with theoretical results or experiments.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Please see our discussion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]

---

Justification:
Proofs are given in the supplemental for stated theorems (or citations given to works where they proved).
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: To the extent possible, we provide as much information as we can on all results to ensure they are reproducible.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
25

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: We have not yet released our code but plan to in the future. We also hope to release the data we used in the future.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
cetaomissssooooooooooooooooooooo cesssonooliooo olicyyyy))llo for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / //nips. cc/ wwcppbllc///aeeecceeeeceecceelee http//uudee//eeeee Sttps///ooSeumiieicc cessonnoolcyyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
To the extent possible, we state all hyperparameters.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Where possible, we state confidence intervals in our results.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi-
dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
26

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
Computational resources are listed in the supplemental.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics httpss//neurs: /ttp///eeeii../.eeeee cc. cc/ Cpuubic//thics//te csGuidel ines?
Answer: [Yes]
Justification: Our work does not raise any ethical concerns, and conforms with all ethical guidelines.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work studies fundamental concepts in machine learning and neuroscience, and we do not believe it has any immediate societal consequences.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
27

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We are not releasing high-risk models or data with this work.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We have credited all those involved with the collection of the data we utilize.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in package should be provided. For popular datasets, paperswi wwwode.hhoe...c. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
28

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets introduced in this work.
Guidelines:
The answer NA means that the paper does not release new assets. Researchers should communicate the details the dataset/code/model as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: No human subjects were involved in this work.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: No human subjects were involved with this work.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
29


Reviews:
[{"review_id": "PQd5X8Ag0N", "reviewer": "Reviewer_YQWs", "summary": "The paper proposes an active learning framework for choosing the next set of neurons to stimulate to best inform a dynamical model of the neural population activity. The active learning procedure takes advantage of the low-rank structure of the dynamical systems model. With synthetic and real datasets, the authors demonstrate that the approach can obtain as much as a 2x reduction in the amount of data required to reach a given predictive power.", "questions": "- Minor comment:\nThe authors could cite more recent uses of data-driven dynamical models for neural data [1, 2].\n\n[1] Karniol-Tambour, O., Zoltowski, D. M., Diamanti, E. M., Pinto, L., Tank, D. W., Brody, C. D., & Pillow, J. W. (2022). Modeling communication and switching nonlinear dynamics in multi-region neural activity. bioRxiv, 2022-09.\n\n[2] Lee, H. D., Warrington, A., Glaser, J., & Linderman, S. (2023). Switching autoregressive low-rank tensor models. Advances in Neural Information Processing Systems, 36, 57976-58010.\n\nIn particular, the authors could consider replacing the low-rank AR models with the model from [2] when fitting models with more lags. While low-rank AR models can be overparameterized easily when incorporating more lags, the model from [2] has fewer parameters than standard low-rank AR models and draws a connection to LDSs.\n\n- How is the lag hyperparameter chosen for the AR models?\n\n- 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.", "limitations": "- As the authors noted, one limitation could be that the low-rank AR-k models may not be effective enough to capture potential nonlinear neural dynamics. Another limitation is that the experiments of the paper are offline. It will be important to test the method in real-time during closed-loop experiments.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Clarity: The paper is well-written, with appropriate equations and figures. Sections 3 and 4 especially present complex mathematical derivations and theorems with details, allowing readers with little expertise to follow the paper.\n\n- Significance: Due to the enormous space of potential stimulations and the time-consuming nature of neuroscience experiments, a reduction in the number of experiments via active learning is crucial for fast advancement in neuroscience.", "weaknesses": "- 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?\n\n- 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.\n\n[1] Jha, A., Ashwood, Z. C., & Pillow, J. W. (2024). Active Learning for Discrete Latent Variable Models. Neural computation, 36(3), 437-474.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 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?\n\nIn 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.\n\nWe 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.\n\n> 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.\n\nWe 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.\n\nSecond, 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. \n\nIf the reviewer has additional benchmarks that would be appropriate to consider with that we have missed, we would be happy to consider.\n\nPukelsheim, Friedrich. Optimal design of experiments. Society for Industrial and Applied Mathematics, 2006.\n\nSimchowitz, Max, et al. \"Learning without mixing: Towards a sharp analysis of linear system identification.\" Conference On Learning Theory. PMLR, 2018.\n\nChaloner, Kathryn, and Isabella Verdinelli. \"Bayesian experimental design: A review.\" Statistical science (1995): 273-304.\n\n> Minor comment: The authors could cite more recent uses of data-driven dynamical models for neural data…\n\nWe 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.\n\n> How is the lag hyperparameter chosen for the AR models?\n\nWe 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.\n\n> 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.\n\nPlease 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "dAf7OCL4h3", "reviewer": "Reviewer_cE26", "summary": "The advent of holographic optogenetics has brought about an unprecedented level of specificity in the way we stimulate and measure the activity across the neural population. The authors propose methods for efficiently determining effective photostimulation patterns to study neural population dynamics using two-photon holographic optogenetics. They propose a novel active learning approach that leverages a low-rank linear dynamical systems model to determine patterns of stimulation with the highest information capacity aiming to characterize the neural population dynamics by minimizing the necessity for data collection.", "questions": "Part of Sector 4.1 seems more like a methods sector. The authors could try to make the algorithm explanation a bit more intuitive and justify a bit more the rationale behind each equation.\n\nThe 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. It would be interesting to try different causal connectivity matrices with more complex data. I understand though that such datasets might not be available.", "limitations": "The authors provide 3 limitations of their approach, namely the uniformity in which the causal connectivity matrix is being approached across all neurons, the applicability of their approach to non-linear dynamics and whether their algorithm can be used online during the experiment.\n\n 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. For example, the neural activity structure of sensory areas could differ significantly from that of motor areas. \n\nFinally, they discuss that the don’t have a clear validation test with new recordings in which the optimal stimulation patterns could be utilized. \n\nIn any case, the direction of this paper is very important for future experiments.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors use holographic 2-photon calcium imaging in order to optimize the pattern of photo-stimulation, something which will be useful for many subsequent studies and experiments.\nThe way in which the authors recover low-rank matrices with the novel use of the nuclear norm and also bound those low rank matrices seems to work well for linearly approachable problems and can also be applied further to other scientific domains.\nFurthermore, they explain the reduced autoregressive model fairly well and provide proofs for 2 new theorems that will ultimately bound the low rank matrices.\nStudying the aforementioned low rank matrices through the lens of active learning gives a novel perspective on how to optimize data acquisition and handling. The paper also makes clear how active learning is used to effectively opt for those stimulations that accelerate the dynamics estimation.                                                                                \nIt appears that this method can significantly reduce the data required for accurate estimation, which is important given the time-limits that these experiments enforce.", "weaknesses": "The dataset is very simple and does not contain anything to enrich the dynamics of the neural population except for random noise. There is not a task in which the animal is involved, at least as presented in the text. The modality of choice is movement which has fairly straightened out dynamics and well characterized, low dimensional neural trajectories. Will the approach work the same in a sensory modality? This could significantly cut off from the model’s wide applicability.\n\nIf the model is used in a more complex modality, the structure of the causal connectivity matrix should be significantly altered, thus current assumptions might not hold true. \n\nSome of the figures are not adequately explained in their legends and some of the ranks are obsolete (Figure 3 and Figure 4)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The dataset is very simple and does not contain anything to enrich the dynamics of the neural population except for random noise…\n\n> 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…\n\n> 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…\n\nWe 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.\n\n> Some of the figures are not adequately explained in their legends and some of the ranks are obsolete (Figure 3 and Figure 4)\n\nWe 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.\n\n> 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.\n\nWe 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "W9LNuhX19w", "reviewer": "Reviewer_d6y7", "summary": "The authors develop active learning techniques to design photostimulation experiments combined with 2p imaging to uncover dynamical systems model of brain activity. To this end, the authors employ low rank matrix recovery techniques. They demonstrate their approach on a dataset from mouse motor cortex.", "questions": "* 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.\n* 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. \n* Line 179: Why is it so important to have an error bound for active learning as compared to other settings? \n* 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.\n* Line 205: What is meant by “played isotropically”? \nAlgorithm 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?\n* Line 248: Is this the same T as in Algorithm 1?", "limitations": "The limitations are very briefly discussed but sufficient.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "* Develops a low-rank matrix estimation framework for dynamical systems of neural activity upon photostimulation.\n* Develops an active learning approach to effectively identify components of the system.\n* The authors demonstrate the usefulness of their approach using simulated experimental settings, which make use of real data.", "weaknesses": "* Fig. 2: It is unclear why the GRU is shown if it is not used later.\n* Line 21: I feel like have read the first sentence many times, almost every year at NeurIPS, in papers on computational neuroscience, where dynamics of neural populations is replaced by whatever the paper studies. You could simply start by the third or the forth sentence without any loss in information.\n* 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. \n* No errorbars in Figure 3 and 4.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Fig. 2: It is unclear why the GRU is shown if it is not used later.\n\nWe 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.\n\n> Line 21: I feel like have read the first sentence many times, almost every year at NeurIPS, in papers on computational neuroscience…\n\nWe thank the reviewer for pointing this out and will revise for the final version accordingly.\n\n> 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. \n\nWe 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.\n\nIntuitively, 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.\n\n> No errorbars in Figure 3 and 4.\n\nFigures 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.\n\n> 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.\n\nWe 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.\n\n> 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. \n\nTheorem 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.\n\nFor 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.\n\n> 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.\n\nBy 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.\n\n> 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?\n\nBy “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.\n\n> Line 248: Is this the same T as in Algorithm 1?\n\nYes.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank each of the reviewers for their helpful feedback, and will work to incorporate all suggestions in the final version. We have addressed specific questions in the following, and are also attaching a pdf with additional visualizations, as requested in some of the reviews.

Meta-review: This paper describes novel active learning methods for optogenetic perturbation experiments aimed at reducing the amount of data required to infer accurate models of neural dynamics.  All three reviewers felt that it was above the threshold for acceptance, and I'm pleased to report that it has been accepted to NeurIPS.  Congratulations!  Please revise the manuscript according to the reviewer comments and discussion points.
Decision: Accept (poster)