Active design of two-photon holographic stimulation for identifying neural population dynamics

Authors: andrew  wagenmaker
lu  mi
georgia  tech
marton  rozsa
matthew s bull
karel  svoboda
kayvon  daie
matthew d golub
kevin  jamieson
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 22:55:30

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.

Main paper: 
 introduction: 
 Neural population dynamics describe how the activities across a population of neurons evolve over time due to local recurrent connectivity and inputs to the population from other neurons or 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 a multitude of hypotheses about how specific neural populations support the neural computations that underlie, for example, motor control [2] [3] [4] , motor timing [5, 6] , decision making [7] [8] [9] [10] , working memory [11] , social behavior [12] , and learning [13] [14] [15] [16] .The traditional approach to data-driven modeling of a neural population typically involves two separate stages. First, neural population activity is recorded while an animal performs a task of interest. Then, a dynamical systems model is fit to the recorded neural responses [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [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 collection-oversampling in some parts of neural activity space while altogether missing others. Given constraints on time and resources in neurophysiological 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 [34] [35] [36] [37] [38] [39] [40] [41] . 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] [43] [44] [45] [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 connectivity. First, we introduce a low-rank autoregressive model that captures low-dimensional structure in neural population dynamics 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 as a simulator of the true dynamics, and 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 on a novel analysis of nuclear-norm 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. 
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] [48] [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 a half [50] [51] [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 measurements we can take is constrained by the physical limits of the photostimulation system. Second, as we aim to adapt and actively learn these 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 a quantification of the estimation error in terms of the precise individual measurements taken (i.e., in contrast to a more global property like RIP).Active Learning and Low-Rank Estimation. The active learning literature is vast, and a 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 system identification problem-is central to many areas of engineering and science [53] . The problem of actively designing inputs to effectively estimate the parameters of a dynamical system has been studied extensively for decades [54] [55] [56] [57] [58] [59] [60] [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 structure, in particular works on low-rank bandits [67] [68] [69] [70] . While the setting considered in these works is somewhat different-they 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-only applying to certain "hard" problems-and do not address the matrix recovery problem.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 are not discussed in the paper. • The authors are encouraged to create a 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 on 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 a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not 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 to address problems of privacy and fairness. • 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
preliminaries: 
 Dataset Details. Neural population activity was recorded in mouse motor cortex using two-photon calcium imaging at 20Hz of a 1mm×1mm field of view (FoV) containing 500-700 neurons. Each recording spanned approximately 25 minutes and 2000 photostimulation trials. In each trial, a 150ms photostimulus was delivered and was followed by a 600ms response period before the next trial began. Each photostimulus targeted a group of 10-20 randomly selected neurons, and a total of 100 unique photostimulation groups were defined for each experiment (≈ 20 trials per group). We evaluate our techniques on four such datasets.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 2 . We describe the model details next.At discrete time t ∈ N, we denote the true neural activity across the d imaged neurons as x t ∈ R d , the noisy, measured activity as y t ∈ R d , and the photostimulus intensity applied across those same d neurons as u t ∈ R d . Applying stimulus u t influences the measured neural activity at the next timestep y t+1 . However, just the snapshot y t 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: where A s ∈ R d×d and B s ∈ R d×d describe the coupling between neurons and stimulus at the time lag of s timesteps, s = 0, . . . , k -1, and offset v ∈ R d accounts for baseline neural activity. Given input-observation pairs {(u t , y t )} t , the coefficients {(A s , B s ) k-1 s=0 , v} 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 ).EQUATIONNeural population dynamics are frequently reported as residing in a subspace of lower dimension than the total number of recorded neurons [2, 7, 17, [72] [73] [74] [75] [76] [77] . The population dynamics in our datasets are consistent with such low-dimensional structure, as indicated by the singular value spectrum in Figure 1(c ). Inspired by this observation, we introduce a set of low-rank dynamical models, where each matrix of {(A s , B s ) k-1 s=0 } is re-defined as diagonal plus low-rank. Explicitly, we parameterize A s = D As + U As V ⊤ As and B s = D Bs + U Bs V ⊤ Bs , where D ∈ R d×d with D ij = 0 for all i ̸ = j, U ∈ R d×r , and V ∈ R d×r for predefined rank r. The diagonal matrices account for substantial autocorrelation in each neuron's activity (D As ) and for the reliable response of each neuron to direct photostimulation (D Bs ), whereas the low-rank matrices (U V ⊤ ) confer coupling between neurons. To fit these parameters, we optimize the following objective function with gradient descent over all parameters:minimize As,Bs∈R d×d ,s=0,...,k-1,v∈R d T t=1 y t+1 - k-1 s=0 A s y t-s - k-1 s=0 B s u t-s -v) 2 .(3.2)Figure 2 shows that these low-rank models perform comparably to the full rank versions in terms of predictive performance; indeed the rank r = 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.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 a causal connectivity matrix, H ∈ R d×d , to be the mapping such that Hu ∈ R d quantifies the total response (across time) of each neuron to a single-timestep photostimulus u. That is, ∞ t=1 x t = Hu, where x 1 , x 2 , x 3 , . . . are the neural activities generated by the population dynamics if u 0 = u, u t≥1 = 0, and x t≤0 = x ∞ is the steady state or resting state of the system subject to no photostimulation. If the dynamics are linear, or more specifically follow (3.1), 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, if it is obtained from a low-rank AR-k model, it too will exhibit low-rank structure.In our experimental paradigm, photostimulation acts as a causal perturbation to the population dynamics, 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 dynamics, 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 our setting.To fit H, we could first fit { A s , B s } k-1 s=0 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 u asymptotes after some finite number of timesteps τ . Thus, we can apply some photostimulus u ∈ R d at time t = 0 and then measure the total response z = τ t=1 y t , where y t ∈ R d is the noisy measurement of the true neural response x t . If we repeat this for many pairs {(u n , z n )} n then we can approximate H asH := arg min H ′ n ∥z n -H ′ u n ∥ 2 2 .In this work we adopt this latter approach. Since we believe H 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 {u n } n 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 dynamics as well.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Question: Does the paper provide open access to the data and code, with sufficient instructions to 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. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
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 ∈ R d and observe the neural population response z ≈ 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 we present novel results characterizing the estimation error of the nuclear norm regression estimator, and in Section 4.2 present an algorithm motivated by these results which seeks to actively estimate low-rank matrices. These results will directly motivate a procedure for designing photostimulation inputs.To demonstrate the generality of our results, in Section 4.1 we consider a general matrix regression setting. In particular, let Θ ⋆ ∈ R d1×d2 be a rank r (potentially non-square) matrix, φ n ∈ R d1×d2 some input matrix, and assume scalar observations:EQUATIONwhere ⟨Θ ⋆ , φ n ⟩ = tr(Θ ⊤ ⋆ φ n ) for tr(•) the trace of a matrix. Note that the setting considered in Section 3.2 is a special case of this observation model with Θ ⋆ ← H and, for each input stimulation u, measuring the response of (4.1) to d inputs φ j of the form φ j ≡ e j u ⊤ for j = 1, . . . , d.Matrix Notation. We let ∥ • ∥ F , ∥ • ∥ op , ∥ • ∥ * denote the Frobenius, operator, and nuclear norm of a matrix, respectively. † denotes the pseudo-inverse of a matrix. vec(•) denotes the vectorization of a matrix, and mat(•) the inverse of the vectorization. We also let △ U denote the simplex-the set of distributions-over a set U.We are interested in understanding how we can effectively take into account the low-rank structure of Θ ⋆ , if our goal is to estimate Θ ⋆ from the observations of (4.1). To this end, we consider the following nuclear-norm constrained least-squares estimator for Θ ⋆ :EQUATIONwhere here we let Φ(Θ ⋆ ) ∈ R N denote the vector where the nth element is ⟨φ n , Θ ⋆ ⟩, and z = Φ(Θ ⋆ ) + η the vector of observations, for η the vector with elements η n . Define Θ ⋆ = U ΣV ⊤ as the skinny SVD such that U ∈ R d1×r , V ∈ R d2×r , and consider the linear projection operators P ⊥ , P ∥ : R d1×d2 → R d1×d2 defined as:P ⊥ (M ) := (I -U U ⊤ )M (I -V V ⊤ ) and P ∥ (M ) := M -P ⊥ (M ),for any M ∈ R d1×d2 . We call P ∥ the projection onto the tangent space of Θ ⋆ . Note that the dimension of the range ofP ∥ is equal to just r(d 1 + d 2 ) -r 2 ≪ d 1 d 2 .We are now ready to state our main result on the estimation error of Θ, for Θ as defined in (4.2). Theorem 1.Define µ := ∥(Φ * Φ) 1/2 ((P ∥ Φ * ΦP ∥ ) † ) 1/2 ∥ op .Then with probability at least 1 -2δ:∥ Θ -Θ ⋆ ∥ F ≤ 4 1 -µ tr (P ∥ Φ * ΦP ∥ ) † + 2∥(P ∥ Φ * ΦP ∥ ) † ∥ op log d1d2 δ + 4∥P ⊥ (Φ * Φ) 1/2 P ⊥ ∥ op ∥(P ∥ Φ * ΦP ∥ ) † ∥ op ( d 1 + d 2 + 2 log 1 δ ),where here Φ * Φ(M ) := n φ n ⟨φ n , M ⟩ and tr(•) describes the sum of the eigenvalues of the linear operator(P ∥ Φ * ΦP ∥ ) † : R d1×d2 → R d1×d2 .Theorem 1 provides a precise bound on the estimation error of the nuclear norm estimator under arbitrary inputs {φ n } 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 problem 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 a special case of the setting considered in (4.1) with Θ ⋆ ← H, Theorem 1 provides a 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 φ n in the tangent space of Θ ⋆ . Indeed, if [w 1 , . . . , w d1 ] and [v 1 , . . . , v d2 ] are the left and right singular vectors of the full SVD of Θ ⋆ , and L ∈ R d1d2×r(d1+d2)-r 2 is a matrix with orthonormal columns vec(w i v ⊤ j ) for (i, j) : {i ≤ r} ∪ {j ≤ r}, thentr (P ∥ Φ * ΦP ∥ ) † = tr L ⊤ N n=1 vec(φ n )vec(φ n ) ⊤ L † ,so we see that the estimation error depends only on the scaling of N n=1 vec(φ n )vec(φ n ) ⊤ in the space spanned by vec(u i v ⊤ j ) for i ≤ r or j ≤ r-the tangent space to Θ ⋆ . As an example of how this scales, assume that for n = 1, . . . , N the entries of each φ n are IID N (0, 1) andN ≥ r(d 1 + d 2 ) -r 2 . Then µ ≈ 0, tr (P ∥ Φ * ΦP ∥ ) † ≈ r(d1+d2)-r 2 N , ∥(P ∥ Φ * ΦP ∥ ) † ∥ op ≈ 1 N , and ∥P ⊥ (Φ * Φ) 1/2 P ⊥ ∥ op ≈ √ N . This translates to a bound of ∥ Θ -Θ ⋆ ∥ 2 F ≤ r(d1+d2)-r 2 +log(1/δ) N. Critically, we see that this does not scale with the total number of parameters, d 1 d 2 , but instead with r(d 1 + d 2 ), 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 Θ,E[∥ Θ -Θ ⋆ ∥ 2 F ] ≥ tr (P ∥ Φ * ΦP ∥ ) † .Given the above characterization, we turn now to the active learning problem: how can we best choose our inputs φ n to speed up estimation error of Θ ⋆ ? 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 u n ∈ R d2 and observe z n = Θ ⋆ u n + η n , for η n ∼ N (0, I d1 ). A single vector observation corresponds to observing d 1 observations from (4.1), the responses to the matrix inputs φ j ≡ e j u ⊤ n for j = 1, . . . , d 1 . Assume that Θ ⋆ is rank r and let V 0 := [v 1 , . . . , v r ] denote the first r right singular vectors of the full SVD of Θ ⋆ . Then we have that:tr (P ∥ Φ * ΦP ∥ ) † = (d 1 -r) • tr (V ⊤ 0 Σ N V 0 ) † + r • tr (Σ N ) † , Σ N := N n=1 u n u ⊤ n . (4.3) This calculation, combined with Theorem 1, shows that the estimation error of Θ ⋆ scales with a weighting of two terms: one quantifying the amount of input energy we put into directions spanned by the topr 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 0 is weighted by a factor of d 1 -r ≈ d 1 , much larger weight than the weight of r given to the term quantifying the isotropic input energy. This suggests that, to minimize the estimation error of Θ ⋆ , we should focus a large portion of our sampling budget to target the directions spanned by the top-r right singular vectors of Θ ⋆ . This strategy admits a transparent intuition. If Θ ⋆ is rank-r and some vector u is orthogonal to the topr right singular vectors of Θ ⋆ , then Θ ⋆ u = 0. Thus, if we know what subspace the top-r right singular vectors of Θ ⋆ span, playing u orthogonal to this subspace gives us no additional information about Θ ⋆ ; in this case we should instead play u 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 Θ ⋆ , for which playing inputs isotropically is optimal.In general, as we do not know Θ ⋆ , we do not know V 0 , and so cannot directly compute inputs minimizing (4.3). To circumvent this, we consider the following iterative procedure, which alternates between obtaining an estimate of Θ ⋆ , Θ, and then playing the inputs that would minimize the estimation error-minimize (4.3)-if Θ were the true parameter. We present this procedure in Algorithm 1.1: input: horizon N , feasible inputs U, rank r, feasible set K 2: Θ 1 ← I, D ← ∅ 3: for ℓ = 1, 2, 3, . . . , ⌈log 2 N ⌉ do 4:Let V 0 denote the top-r right singular vectors of Θ ℓ and Λ(λ) := u∈U λ u uu ⊤ , solve:λ V ℓ ← arg min λ∈△ U tr ( V ⊤ 0 Λ(λ) V 0 ) † , λ unif ℓ ← arg min λ∈△ U tr Λ(λ) † 5:For 2 ℓ steps, play inputu n ∼ 1 2 λ V ℓ + 1 2 λ unif ℓ , add observations to D 6: Update estimate of Θ ⋆ : Θ ℓ+1 ← arg min Θ∈K (u,z)∈D ∥z -Θu∥ 2 F 7: return Θ ℓ+1 .At every iteration ℓ, Algorithm 1 computes two distributions over inputs: λ V ℓ , which targets the top-r right singular vectors of our current estimate of Θ ⋆ , and λ unif ℓ , which plays inputs isotropically, covering all directions. Rather than playing these distributions according to the precise weighting given in (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 Θ ⋆ , λ V ℓ is not guaranteed to target the correct directions, especially in early iterations. λ unif ℓ plays inputs in every direction, however, and thus, even if λ V ℓ is not aligned to the top-r right singular vectors of Θ ⋆ , will ensure sufficient energy is still being played in the correct directions to allow for learning. Given this, we increase the weight of playing λ unif ℓ relative to that prescribed by (4.3) .Note that the computation of the optimal inputs is a form of A-optimal experiment design [79] , which in general can be efficiently solved by, for example, the Frank-Wolfe algorithm [80] . Furthermore, efficient procedures for solving nuclear-norm regression problems exist, allowing us to estimate Θ ℓ+1 on line 6 efficiently [81] . We remark that Algorithm 1 takes as input r, the rank of Θ ⋆ , and K, which requires knowledge of ∥Θ ⋆ ∥ * . In general, when these quantities are unknown, they can be chosen via standard cross-validation procedures.We emphasize again that the setting considered here corresponds precisely to the setting considered in Section 3.2 with Θ ⋆ ← H, u n the input stimulation patterns, and z n the observed neural response to input u n . As such, if the causal connectivity H is low rank, Algorithm 1 and the preceding results provide a methodology to select input stimuli to most efficiently estimate H. In the following section, we will apply this to our photostimulation datasets. 
active learning for estimating neural population dynamics: 
 We return now to the problem of photostimulus design for learning neural population dynamics, and seek to apply the insights of Section 4 to 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 a simulator for the true dynamics, and then demonstrate that we can learn the causal connectivity matrix H 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 3) 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 to observe, training a 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 non-adaptive approaches.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 we treat them as a simulated representation of our true dynamics, allowing us to query them arbitrarily as a 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-k model to each dataset as described in Section 3.1. In all cases we use an AR-k model with order k = 4. We do one run of the experiments using low-rank model parameters U V ⊤ with rank 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 at each new trial.In practice, both the magnitude of the stimuli and number of neurons stimulated at each timestep are constrained by the photostimulation 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 U in Algorithm 1 to be U := {u ∈ [0, 1] d : ∥u∥ 1 ≤ γ}, for some value γ > 0 (which we set to γ = 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 2γ-sparse. As baseline methods, we consider the following:• Random Stimulation: At each trial n, choose γ neurons at random, and set corresponding elements of u n to 1.• Uniform Stimulation: Compute λ unif as in Algorithm 1 and play inputs u n ∼ λ unif for all n. Treating this as a simulator of the true dynamics, we compare our active stimulation design procedure (Active, Algorithm 1) 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 1 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 ∥M ⊙(H-H)∥F, for H our estimate of H, M a matrix with all entries 1 except its diagonal, which is 0, and ⊙ element-wise multiplication.Experiment Results. We present our 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, on Mouse 1 and both datasets for Mouse 3, we observe a gain of between 1.5-2× over baselines-that is, to achieve a given estimation error, our approach requires between 1.5-2× fewer samples than baseline methods. This demonstrates the effectiveness of our active learning procedure for estimating low-rank matrices-our method is able to exploit the low-rank structure present in the underlying dynamics to 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 to speed up the estimation of the dynamics.As described in Section 3, 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 to simulate this process using real experimental data, but offline, by choosing the ordering of the trials available in our pre-collected datasets. This serves as a testbed for active learning procedures: if we can more efficiently learn models in this offline setting, that is a 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 to 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 set aside a test set containing all 20 repeated trials of those photostimuli. This creates an 80%/20% train-test split of non-overlapping stimulus patterns. For D train and D test our train and test datasets, respectively, we consider the following query model: Each 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 of Algorithm 1 (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, as well as splits where Active has the largest and smallest improvement over Random. We plot error bars denoting 1 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 of 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 1 to learning a full dynamical system, we choose our inputs to target the right singular vectors of B s in (3.1). As a baseline method, we consider the procedure which randomly chooses an unobserved segment from D train at each iteration. We run the above experiment for 20 different randomly generated train-test splits on each dataset, and present our results in Figure 4 , providing the results for the average performance over the train-test splits, as 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 it is able to give a gain of up to a factor of 2× in the number of samples required over the random baseline, 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 can improve estimation rates, and, we 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 inputs 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 or changes in connectivity due to learning. Second, while we found that linear dynamics fit our data remarkably well, this may not always be the case. Does our methodology effectively scale to nonlinear dynamics? Finally, our real-data experiments were performed offline. Future work may explore running our algorithm online during closed-loop photostimulation experiments.As Θ is a minimizer, we have that ∥Φ( Θ) -z∥ 2 2 ≤ ∥Φ(Θ ⋆ ) -z∥ 2 2 . Consequently, ∥Φ(Θ ⋆ ) -z∥ 2 2 < min Θ⋆+∆∈K:∥∆∥F≥ρ ∥Φ(Θ ⋆ + ∆) -z∥ 2 2 =⇒ Θ ̸ ∈ {Θ ⋆ + ∆ ∈ K : ∥∆∥ F ≥ ρ} =⇒ ∥ Θ -Θ ⋆ ∥ F < ρ.Thus, our strategy will attempt to find a minimum ρ > 0 that makes the first expression true. Equivalently, we show that the following quantity is non-positive∥Φ(Θ ⋆ ) -z∥ 2 2 - min Θ⋆+∆∈K:∥∆∥F≥ρ ∥Φ(Θ ⋆ + ∆) -z∥ 2 2 = - min Θ⋆+∆∈K:∥∆∥F≥ρ 2⟨Φ(∆), Φ(Θ ⋆ ) -z⟩ + ∥Φ(∆)∥ 2 2 = max M +∆∈K:∥∆∥F≥ρ 2⟨Φ(∆), η⟩ -∥Φ(∆)∥ 2 2 = max t≥ρ max M +∆∈K:∥∆∥F=t 2⟨Φ(∆), η⟩ -∥Φ(∆)∥ 2 2 .Intuitively, when t is large, the quadratic term will dominate the inner product making the entire expression non-positive. We will employ a geometric fact about the nuclear norm ball. Lemma A.1.Let K = {Θ : ∥Θ∥ * ≤ ∥Θ ⋆ ∥ * }. If Θ ⋆ + ∆ ∈ K then ∥P ⊥ (∆)∥ * ≤ -⟨∆, U V ⊤ ⟩ ≤ ∥P ∥ (∆)∥ F . Proof. If ∥Θ ⋆ + ∆∥ * > ∥Θ ⋆ ∥ * then Θ ⋆ + ∆ ̸ ∈ K.By the convexity of the nuclear norm ball, we have that∥Θ ⋆ + ∆∥ * ≥ ∥Θ ⋆ ∥ * + ⟨∆, U V ⊤ ⟩ + ⟨∆, W ⟩ ∀W : W = P ⊥ (W ) , ∥W ∥ 2 ≤ 1.Consequently, as the dual norm to ∥ • ∥ 2 is the nuclear norm ∥ • ∥ * we have∥Θ ⋆ + ∆∥ * ≥ ∥Θ ⋆ ∥ * + ⟨∆, U V ⊤ ⟩ + ∥P ⊥ (∆)∥ * . Thus, if ⟨∆, U V ⊤ ⟩ + ∥P ⊥ (∆)∥ * > 0 then Θ ⋆ + ∆ ̸ ∈ K. Consequently, if Θ ⋆ + ∆ ∈ K then ⟨∆, U V ⊤ ⟩ + ∥P ⊥ (∆)∥ * ≤ 0.Recalling that ∥M ∥ F ≤ ∥M ∥ * for any matrix M , an interesting consequence of the above lemma is that∥P ⊥ (∆)∥ F ≤ ∥P ∥ (∆)∥ F which implies ∥P ∥ (∆)∥ 2 F ≤ ∥∆∥ 2 F ≤ 2∥P ∥ (∆)∥ 2 F .That is, the total error is dominated by the error in the tanget space of Θ ⋆ .Applying this lemma, we have that maxM +∆∈K:∥∆∥F=t 2⟨Φ(∆), η⟩-∥Φ(∆)∥ 2 2 ≤ max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(∆), η⟩ -∥Φ(∆)∥ 2 2 = max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)) + Φ(P ⊥ (∆)), η⟩ -∥Φ(∆)∥ 2 2 ≤ max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)), η⟩ -∥Φ(∆)∥ 2 2 + max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ⊥ (∆)), η⟩where the equality uses the fact that ∆ = P ∥ (∆) + P ⊥ (∆) and the linearity of Φ. For any v ∈ R N we have⟨Φ(∆), v⟩ = ⟨∆, N n=1 X n v n ⟩ =: ⟨∆, Φ * (v)⟩where Φ * denotes the adjoint of Φ. We also recognize that the operator Φ * Φ : R d1×d2 → R d1×d2 is also linear, defined asΦ * Φ(M ) = Φ * ({⟨X n , M ⟩} n ) = N n=1 X n ⟨X n , M ⟩. Consequently (Φ * Φ) 1/2 iswell-defined and is the same operator as Φ * Φ after taking the square root of its eigenvalues.The next three lemmas bound the two terms of above. Combining them yields the result of the theorem.Lemma A.2. With probability at least 1 -δ we have max∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ⊥ (∆)), η⟩ ≤ 2t∥P ⊥ (Φ * Φ) 1/2 P ⊥ ∥ op ( d 1 + d 2 + 2 log(1/δ)).Proof. Computing this term amounts to bounding a Gaussian width. Begin by recognizing that by the non-expansive property of projections, ∥P ∥ (∆)∥ F ≤ ∥∆∥ F ≤ t which results in the simplification:max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ⊥ (∆)), η⟩ ≤ max ∆:∥P ⊥ (∆)∥ * ≤t 2⟨Φ(P ⊥ (∆)), η⟩ = max ∆:∥P ⊥ (∆)∥ * ≤t 2⟨ N n=1 η n X n , P ⊥ (∆)⟩ = max ∆:∥P ⊥ (∆)∥ * ≤t 2⟨P ⊥ ( N n=1 η n X n ), P ⊥ (∆)⟩ ≤ 2t∥P ⊥ ( N n=1 η n X n )P ⊥ ∥ opusing the fact that the operator norm and nuclear norm are dual to each other. Note thatN n=1 η n X n = mat( N n=1 η n vec(X i )) = mat(( N n=1 vec(X n )vec(X n ) ⊤ ) 1/2 vec(η ′ )) =: (Φ * Φ) 1/2 (η)where η ′ ∈ R d1×d2 with η ′ i,j ∼ N (0, 1). We then observe that∥P ⊥ ( N n=1 η n X n )P ⊥ ∥ op = ∥P ⊥ (mat(( N n=1 vec(X n )vec(X n ) ⊤ ) 1/2 vec(η ′ )))P ⊥ ∥ op ≤ ∥P ⊥ (Φ * Φ) 1/2 P ⊥ ∥ op ∥η ′ ∥ opwhich completes the proof of the first claim. Recognizing that ∥η ′ ∥ op is just the largest singular value of a Gaussian matrix, we find that E[sup ∥u∥2≤1,∥v∥2≤1 ⟨uv ⊤ , η⟩] ≤ √ d 1 + √ d 2 by Exercise 5.14 of [52] . Applying a sub-Gaussian tail bound completes the proof.Lemma A.3. Define µ := ∥(Φ * Φ) 1/2 ((P ∥ Φ * ΦP ∥ ) † ) 1/2 ∥ op . Then max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)), η⟩ -∥Φ(∆)∥ 2 2 ≤ max ∆:∥P ∥ (∆)∥F≥t/ √ 2 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2 Proof. Recognizing that Φ(∆) ∈ R N we have ∥Φ(∆)∥ 2 2 = ∥Φ(P ∥ (∆) + P ⊥ (∆))∥ 2 2 = ∥Φ(P ∥ (∆)) + Φ(P ⊥ (∆))∥ 2 2 = ∥Φ(P ∥ (∆))∥ 2 2 + ∥Φ(P ⊥ (∆))∥ 2 2 + 2⟨Φ(P ∥ (∆)), Φ(P ⊥ (∆))⟩.To aid in readability, we make a number of notational modifications. First, we drop parentheses so that Φ(P ∥ (∆)) is just notated as ΦP ∥ ∆. Second, we define M †/2 := (M † ) 1/2 where M † is the pseudoinverse. If Φ * Φ is invertible restricted to the range of P ∥ , then P ∥ ∆ = (P ∥ Φ * ΦP ∥ ) †/2 (Φ * Φ) 1/2 P ∥ ∆ for all ∆. Thus,|⟨Φ(P ∥ (∆)),Φ(P ⊥ (∆))⟩| = |⟨ΦP ∥ ∆, ΦP ⊥ ∆⟩| = |⟨P ∥ ∆, (Φ * Φ)P ⊥ ∆⟩| = |⟨(Φ * Φ) 1/2 P ∥ ∆, (Φ * Φ) 1/2 P ⊥ ∆⟩| = |⟨(Φ * Φ) 1/2 P ∥ (P ∥ Φ * ΦP ∥ ) †/2 (Φ * Φ) 1/2 P ∥ ∆, (Φ * Φ) 1/2 P ⊥ ∆⟩| ≤ ∥(Φ * Φ) 1/2 (P ∥ Φ * ΦP ∥ ) †/2 (Φ * Φ) 1/2 P ∥ ∆∥ F ∥(Φ * Φ) 1/2 P ⊥ ∆∥ F ≤ ∥(Φ * Φ) 1/2 (P ∥ Φ * ΦP ∥ ) †/2 ∥ op ∥(Φ * Φ) 1/2 P ∥ ∆∥ F ∥(Φ * Φ) 1/2 P ⊥ ∆∥ Fwhere the last two lines follow from Cauchy-Schwartz. Observe that for µ < 1 we havea 2 + b 2 -2abµ = (1 -µ)a 2 + (1 -µ)b 2 + µa 2 + µb 2 -2abµ = (1 -µ)a 2 + (1 -µ)b 2 + µ(a -b) 2 ≥ (1 -µ)a 2 . Thus, if µ := ∥(Φ * Φ) 1/2 (P ∥ Φ * ΦP ∥ ) †/2 ∥ op then ∥Φ(∆)∥ 2 2 ≥ ∥Φ(P ∥ (∆))∥ 2 2 + ∥Φ(P ⊥ (∆))∥ 2 2 -2|⟨Φ(P ∥ (∆)), Φ(P ⊥ (∆))⟩| ≥ ∥Φ(P ∥ (∆))∥ 2 2 + ∥Φ(P ⊥ (∆))∥ 2 2 -2µ∥(Φ * Φ) 1/2 P ∥ ∆∥ F ∥(Φ * Φ) 1/2 P ⊥ ∆∥ F ≥ (1 -µ)∥Φ(P ∥ (∆))∥ 2 2 Thus, max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)), η⟩ -∥Φ(∆)∥ 2 2 ≤ max ∆:∥P ⊥ (∆)∥ * ≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2 ≤ max ∆:∥P ⊥ (∆)∥F≤∥P ∥ (∆)∥F,∥∆∥F=t 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2 ≤ max ∆:∥P ∥ (∆)∥F≥t/ √ 2 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2where we've used the facts that∥ • ∥ * ≤ ∥ • ∥ F and t 2 = ∥∆∥ 2 F = ∥P ∥ (∆)∥ 2 F + ∥P ⊥ (∆)∥ 2 F ≤ 2∥P ∥ (∆)∥ 2 F . Lemma A.4. Let K = min{log 2 (N ), d 1 d 2 }. Then for any α > 0, if t ≥ 1 1 -µ 16tr (P ∥ Φ * ΦP ∥ ) † + 32∥(P ∥ Φ * ΦP ∥ ) † ∥ op log(K/δ) + 2α∥(P ∥ Φ * ΦP ∥ ) † ∥ op 1 -µthen with probability at least 1 -δ we have αt + max∆:∥P ∥ (∆)∥F≥t/ √ 2 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2 ≤ 0.Proof. The linear operator ΦP ∥ : R d1×d2 → R N can be decomposed as ΦP ∥ = N n=1 β n w n ψ n where {w n } n are orthonormal on R T , {ψ n } n are orthonormal linear operators on R d1×d2 , andβ n ≥ 0 are decreasing. For k = 0, 1, . . . , min{log 2 (N ), d 1 d 2 } -1 let W k = [w 2 k , . . . , w 2 k+1 -1 ] so that β 2 k = max ∥∆∥F=1,∥u∥2=1 u ⊤ W ⊤ k ΦP ∥ (∆) ≥ min ∥∆∥F=1,∥u∥2=1 u ⊤ W ⊤ k ΦP ∥ (∆) ≥ β 2 k+1Then max∆:∥P ∥ (∆)∥F≥t/ √ 2 2⟨Φ(P ∥ (∆)), η⟩ -(1 -µ)∥Φ(P ∥ (∆))∥ 2 2 = max ∆:∥P ∥ (∆)∥F≥t/ √ 2 K k=0 2⟨W k W ⊤ k Φ(P ∥ (∆)), η⟩ -(1 -µ)∥W k W ⊤ k Φ(P ∥ (∆))∥ 2 2 = max ∆:∥P ∥ (∆)∥F≥t/ √ 2 K k=0 2⟨W ⊤ k Φ(P ∥ (∆)), W ⊤ k η⟩ -(1 -µ)∥W ⊤ k Φ(P ∥ (∆))∥ 2 2 ≤ K k=0 max ∆:∥P ∥ (∆)∥F≥t/ √ 2 2∥W ⊤ k Φ(P ∥ (∆))∥ 2 ∥W ⊤ k η∥ 2 -(1 -µ)∥W ⊤ k Φ(P ∥ (∆))∥ 2 ≤ K k=0 t √ 2β 2 k+1 ∥W ⊤ k η∥ 2 - t 2 (1 -µ) 2 β 2 2 k+1where the first inequality holds by Cauchy-Schwartz, and the second holds for all t ≥ max k=0,...,K∥W ⊤ k η∥22 √ 2 (1-µ)β 2 k+1 . Moreover, t √ 2β 2 k+1 ∥W ⊤ k η∥ 2 -t 2 (1-µ) 2 β 2 2 k+1 ≤ -αt if t ≥ max k=0,...,K ∥W ⊤ k η∥ 2 2 √ 2 (1 -µ)β 2 k+1 + max k=0,...,K 2α (1 -µ)β 2 2 k+1 = max k=0,...,K 8∥W ⊤ k η∥ 2 2 (1 -µ) 2 β 2 2 k+1 + max k=0,...,K 2α (1 -µ)β 2 2 k+1Note that for K = min{log 2 (N ), d 1 d 2 } we haveP(∪ K k=1 {∥W ⊤ k η∥ 2 ≥ √ 2 k + 2 log(K/δ)}) ≤ P(∪ K k=1 {∥W ⊤ k η∥ 2 ≥ E[∥W ⊤ k η∥ 2 ] + 2 log(K/δ)}) ≤ δ. On this good event, we have that max k=0,...,K ∥W ⊤ k η∥ 2 2 β 2 2 k+1 ≤ max k=0,...,K ( √ 2 k + 2 log(K/δ)) 2 β 2 2 k+1 ≤ max k=0,...,K 2 k+1 β 2 2 k+1 + 4 log(K/δ) β 2 2 k+1 ≤ 2 N n=1 1 β 2 n + max n=1,...,N 4 log(K/δ) β 2 n = 2tr (P ∥ Φ * ΦP ∥ ) † + 4∥(P ∥ Φ * ΦP ∥ ) † ∥ op log(K/δ)where we use the fact thatN n=1 1 β 2 n = K k=0 2 k+1 -1 n=2 k 1 β 2 n ≥ K-1 k=0 2 k β 2 2 k+1 ≥ max k=0,...,K-1 2 k β 2 2 k+1.The photostimulation data were collected from transgenic reporter mice Ai229, which express Crerecombinase-dependent cytosolic GCaMP6m and soma-targeted ChRmine, crossed with the Vglut1cre mouse line. Imaging and photostimulation experiments were performed on a 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://github.com/MouseLand/suite2p).We split each of our photostimulation datasets into non-overlapping 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 5 (out of the 100 total) unique photostimulation patterns and then included all 70-timestep windows about each of the 20 instances of those 5 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 {y t , u t } 70 t=1 , where y t ∈ R d is the recorded neural activity and u t ∈ R d is the photostimulation delivered at time t. During evaluation on a given test window, all models were provided {y t } 4 t=1 and {u t } 70 t=1 to predict {y t } 70 t=5 . Autoregressive-k models: We fit the full-rank AR-k models to training datasets via linear regression by expressingy t+1 = k-1 s=0 A s y t-s + k-1 s=0 B s u t-s + v (B.1)as Y = XW , whereY =     y 1 y 2 . . . y T +1     X =     y 0 y -1 . . . y 1-k u 0 u -1 . . . u 1-k 1 y 1 y 0 . . . y 2-k u 1 u 0 . . . u 2-k 1 . . . . . . . . . . . . . . . . . . 1 y T y T -1 . . . y T +1-k u T u T -1 . . . u T +1-k 1     W =                A 0 A 1 . . . A k-1 B 0 B 1 . . . B k-1 v               and the closed-form solution is W = (X T X) -1 X T Y . 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 PyTorch and ran on a single NVIDIA Tesla T4 GPU.During evaluation of the AR-k models, for each test window we first computed y k+1 given {y t , u t } k t=1 using (B.1). Then for all subsequent predictions, on the right-hand side of (B.1) we replaced all instances of y t with y t for t > k. In this manner, each entire roll-out prediction of {y t } 70 t=k+1 used all photostimulation inputs {u t } 70t=1 , but only the first k timesteps of neural activity {y t } k t=1 . 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 k = 4 initial timesteps of recorded neural activity into a bottlenecked initial state for a 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 to 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 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 a given neuron's measured calcium fluorescence exceeded a predefined threshold. To calculate a receiver operator characteristic (ROC) curve, we enumerated a range of thresholds, normalized by the standard deviation of each neuron's empirical activity distribution, and performed threshold detection separately on the real and model-predicted neural activity traces. We then compute the overall false-positive rate and true-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.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% of 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 of this long test window. We report these results in Figure 5 . As K is defined with respect to the nuclear norm of the true parameter, which we do not assume is known, we run each method with a 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: To choose the input rank of Algorithm 1, we ran our experiment with several different ranks and provide results for the best-performing rank. We found, however, that results are typically robust to the setting of the rank parameter of Algorithm 1, and our choice of r did not significantly impact performance. Furthermore, we believe this could effectively be chosen adaptively. We state our chosen values of r below.ActiveFor all experiments, we add observation noise distributed as N (0, 0. For this experiment, on the data D we observed thus far, we fit the AR-k model described in Section 3.1 with k = 1. 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 baseline method.Given a input response trajectory in the test set, (x 1 , . . . , x 15 ), with input u, to compute the test MSE, we provide our learned dynamics model with the initial state x 1 and input u, and then roll this out for 15 timesteps to generate predictions x 2 , . . . , x 15 . Precisely, if A and B are our estimated parameters, we letx 2 = Ax 1 + Bu, x t+1 = A x t , t ≥ 1.We the compute the MSE on this segment as:1 14 15 t=2 ∥ x t -x t ∥ 2 2It is not immediately obvious how to apply Algorithm 1 to this setting, since we must choose each trajectory sequentially, and once we have observed a trajectory it can no longer be chosen again. Rather than solving the optimization of Algorithm 1 to find the best inputs, we instead seek to iteratively choose the next input that would maximize "information gain" in some sense. In particular, note that applying the Frank-Wolfe algorithm [80] to the objective, if we have inputs U available:min λ∈△ U tr((V ⊤ Λ(λ)V ) -1 ),the update is given by:u i+1 = min u∈U u ⊤ V (V ⊤ Λ(λ i )V ) -2 V ⊤ u λ i+1 ← (1 -γ i )λ i + γ i I{u = u i+1 }for learning rate γ i .In this experiment, we simply choose u n as above, with U the set of remaining active inputs in D train , and Λ(λ n ) replaced with n-1 s=1 u s u ⊤ s . This therefore approximates the solution to the experiment design of Algorithm 1, and has the advantage of being very computationally efficient. Furthermore, we 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 did not find the results particularly sensitive to setting of r. For each dataset, we ran with r ∈ [25, 50, 75, 100, 125, 150] , and include results for the best-performing setting.For both sets 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 we plot the predictions from the estimated model at different MSE values on heldout trials for Mouse 2. 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 be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?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 to a 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, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?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, confidence intervals, or statistical significance tests, at 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 run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to 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. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 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 out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?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, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as 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 it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?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 circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
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 is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 a 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 or 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 to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe 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)? 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 a 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 a best faith effort. 12. Licenses for existing assets Question: 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? 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 or dataset.• The authors should state which version of the asset is used and, if possible, include a 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 the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: - 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?

- 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.

[1] Jha, A., Ashwood, Z. C., & Pillow, J. W. (2024). Active Learning for Discrete Latent Variable Models. Neural computation, 36(3), 437-474.
questions: questions: - Minor comment:
The authors could cite more recent uses of data-driven dynamical models for neural data [1, 2].

[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.

[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.

In 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.

- How is the lag hyperparameter chosen for the AR models?

- 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: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

If 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. 

Some of the figures are not adequately explained in their legends and some of the ranks are obsolete (Figure 3 and Figure 4)
questions: 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.

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. 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: 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.

 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. 

Finally, they discuss that the don’t have a clear validation test with new recordings in which the optimal stimulation patterns could be utilized. 

In any case, the direction of this paper is very important for future experiments.
paper_summary: paper_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.
strength_weakness: strength_weakness: * Fig. 2: It is unclear why the GRU is shown if it is not used later.
* 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.
* 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. 
* No errorbars in Figure 3 and 4.
questions: 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.
* 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. 
* Line 179: Why is it so important to have an error bound for active learning as compared to other settings? 
* 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.
* 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?
* Line 248: Is this the same T as in Algorithm 1?
limitations: limitations: The limitations are very briefly discussed but sufficient.

Comments: comment: I would like to thank the authors for their response to my comments and the clarifications and modifications that will appear in the manuscript. I would like to increase the score from 5 (Borderline accept) to 6 (Weak accept).
I would like to thank the authors for the response to my comments and the clarifications that will appear in the manuscript.