Bayesian target optimisation for high-precision holographic optogenetics
Authors: Marcus Triplett
Marta Agnieszka Gajowa
Hillel Adesnik
Liam Paninski
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Neuroscience
neural stimulation
optogenetics
calcium imaging
TL;DR: A Bayesian approach to optimising stimulation parameters significantly improves the precision of two-photon optogenetics

Abstract: Two-photon optogenetics has transformed our ability to probe the structure and function of neural circuits. However, achieving precise optogenetic control of neural ensemble activity has remained fundamentally constrained by the problem of off-target stimulation (OTS): the inadvertent activation of nearby non-target neurons due to imperfect confinement of light onto target neurons. Here we propose a novel computational approach to this problem called Bayesian target optimisation. Our approach uses nonparametric Bayesian inference to model neural responses to optogenetic stimulation, and then optimises the laser powers and optical target locations needed to achieve a desired activity pattern with minimal OTS. We validate our approach in simulations and using data from in vitro experiments, showing that Bayesian target optimisation considerably reduces OTS across all conditions we test. Together, these results establish our ability to overcome OTS, enabling optogenetic stimulation with substantially improved precision.

Bayesian target optimisation for high-precision holographic optogenetics
Marcus A. Triplett 12..t
Marta Gajowa²
Hillel Adesnik
Liam Paninski 1.2
Department of Statisties, Columbia University Zuckerman Mind Brain Behavior Institute, Columbia University Department of Molecular and Cell Biology, UC Berkeley
marcus. triplett@columiia edu
Abstract
Two-photon optogenetics has transformed our ability to probe the structure and function of neural circuits. However, achieving precise optogenetic control of neural ensemble activity has remained fundamentally constrained by the problem off-target stimulation (OTS): the inadvertent activation of nearby non-target neurons due to imperfect confinement of light onto target neurons. Here propose a novel computational approach to this problem called Bayesian target optimisation. Our approach uses noupaameerrcc Bayesian inference to model neural responses to optogenetic stimulation, and then optimises the laser powers and optical target locations needed to achieve desired activity pattern with minimal OTS. We validate our approach in simulations and using data from in vitro experiments, showing that Bayesian target optimisation considerably reduces OTS across all conditions we test. Together, these results establish our ability to overcome OTS, enabling optogenetic stimulation with substantially improved precision.
Introduction
key technological goal in neuroscience is to gain precise control over the activity of neurons. Cur- rently, of the most promising approaches for achieving such control is two-photon optogenetics [14. This technique relies on the use of two-photon excitation to activate opsin molecules expressed in the somatic membrane, allowing individual neurons be targeted for stimulation. Holographic optogenetics further extends this technique by focusing two-photon excitation into ~ 10 yam disks of light that illuminate all of neuron's opsin molecules parallel [.)). Replicating this illumi- nation pattern across many neurons at once then grants simultancous optogenetic control of entire neural ensembles [0 [2]. However, the spatial precision of holographic optogenetics has remained fundamentally limited by the problem of off-target stimulation (OTS): if multiple opsin-expresse neurons lie in close proximity (in the most extreme case, their membranes are physical contact), two-photon excitation will frequently activate opsin molecules in each neuron when attempting to stimulate only one. This can inadvertently elicit spikes non-target neurons, compromising the precision and specificity of the optogenetic manipulation [3]
Previous attempts overcoming OTS have been either optical or molecular in nature (e.g. by using temporal focusing [2]
6, 7 [4] or improving the soma-targeting of the opsin [15 16 10)), but
have not yet achieved "true" single-cell precision under many experimental conditions. Frequently, experimenters will simply express opsin sparsely in population to avoid OTS, as this reduces the number of nearby oxxin-exp neurons that could be mistakenly activated [17]. However, sparse opsin expression reduces the number of neurons that can potentially be probed during an
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

experiment, and therefore limits the kinds of scientific questions that can be addressed using two- photon optogenetics. Conversely, in an ideal experiment, one could precisely stimulate any neuron despite a moderate or high density of opsin-expressing cells.
Here, develop novel computational strategy for all-optical experiments [1-21 that reduces (and in some cases entirely eliminates) OTS (Figure La. b). The key insight is that heterogeneity in opsin expression and intrinsic excitability causes neurons have different sensitivities to optical stimulation. This implies that in some cases laser power could be safely turned down to prevent the spiking of nearby non-target neuron that is weakly driven by stimulation while still spiking a target neuron that is strongly driven by stimulation. Further, the spatial arrangement of neurons can be exploited by moving holographic target off of neuron's 's nucleus (the typical stimulation point) towards a less dense area (e.g. towards neuropil or less photosensitive neurons) while still exciting a sufficient number of the target neuron's is opsin molecules to elicit a spike. While some prior experimental research has considered cell-specific tuning of laser powers [22], our advance is to identify both the optimal laser powers and target locations automatically using efficient Bayesian techniques.
To achieve highly precise optogenetic ensemble control, we propose two-phase process called Bayesian target optimisation (Bataro). In the first phase, we map the "optogenetic receptive field" (ORF) of each neuron by stimulating at various locations around cell nuclei and a range of laser powers to test whether they elicit spikes. To infer these ORFs without exhaustively testing every combination of location and power, we use Gaussian processes (GPs) to encode prior knowledge about how neurons respond to optogenetic stimulation. Further, we propose to use holographic stimulation to test multiple locations around one or many neurons simultaneously, ensuring that the mapping phase completes quickly. This also allows to model how neurons integrate two-photon excitation from multiple nearby holograms. In the second phase, we then exploit knowledge of the ORFs to optimise holographic targets. We use properties of the GP to infer gradients of objective function at a spatial resolution that exceeds the original sampling resolution, enabling us to optimise for the exact holographic targets and laser powers needed to evoke specific neural activity pattern.
2 Related work
Our work is related to prior research in three areas: GP-based receptive field inference, optimal stimulus design, and statistical methods for estimating connectivity using optogenetic stimulation. Since ref. 123), many studies have used GPs to infer the relationship between neural spiking and multi-dimensional covariates. This includes using GPs to infer place fields 24, 25 and orientation preference maps (26). as well as with Bayesian active learning techniques for inferring sensory receptive fields 27 28 However, none to our knowledge have used GP inference to model responses to optogenetic stimulation. Closed-loop methods for optimal control of neural activity using one- photon or electrical stimulation have been recently explored (2 34], though these approaches lack the necessary spatial specificity for highly precise ensemble control. Ref. (35] considers optimising spike-timing using two-photon stimulation (in addition to electrical stimulation), but focused on single neurons and did not GP estimation methods. Finally, number of papers have developed statistical models of optogenetic data to infer functional or synaptic connectivity [36 but do not consider identifying the exact stimulation parameters needed to evoke specific activity patterns. Here, we unify these three approaches to develop a novel computational framework for optimal holographic stimulation of neural ensembles.
3 Methods
3.1
Optogenetic receptive field model
Given population of N neurons, we first wish to learn ORF functions In that model how each neuron responds to optogenetic stimulation at different laser powers and locations near the cell soma. Ideally, ORF mapping should be completed quickly, stimulating as few times as possible, so that stimulus optimisation and the desired experiment can begin. Therefore, an ORF model should leverage prior knowledge of how neurons respond to stimulation. Importantly, spike probability should increase monotonically with laser power until saturation (4), and the ORF should approximately match the shape of the somatic membrane (though enlarged to account for the ~0 pm diameter of the

---

Holographic ensemble stimulation
sampled
60% sampled
100% sampled
miw datence ((m)
o spike
no
response
opain condenmaation
-20
Target neuran
Naive holo
Optimised holo
distansa (um)
20
godd soyds
-00
distance (jm)
X distance (um)
distance (um)
Figure I: Off-target stimulation with two-thoton optogenetics. (a) Two-photon holographic opto- genetics can be used to elicit spikes in specific ensembles of experimenter-seleee OTS arises due to an inability to confine two-photon excitation to the soma of a target neuron. By repositioning the hologram away from the soma, OTS could be avoided while still activating the target neuron. (c) Data from real two-photon optogenetics experiment showing that at high power (e.g. 40 mW) a neuron can be activated from 15-20 jmm away, though this depends on the specific pattern of opsin expression at the soma and proximal dendrites. Red and gray circles indicate locations where stimulation resulted in successful or unsuccessful spikes. Colour map shows the inferred probability of spiking (i.e., the ORF) using the log- barrier Newton method from (Sbseecttn 3.2] As the number of sampled locations and laser powers increases, the GP model adapts to the particular ORF shape (see supplementary material for the prior mean). This shape can then be exploited to precisely optimise holographic stimuli. Data from PV-neuron in 1,/3 of V1 expressing the soma-targeted, excitatory opsin ChroME2f [4]
neurons. (b)
holographic disk 22)). However, any residual expression of opsin molecules in the proximal dendrites of a neuron creates unique differences in its ORF shape compared to other neurons (Figure Te). Further, neurons can express opsin in variable concentrations and vary their intrinsic excitability, leading to an unpredictable dependence of spiking on power. Thus, we need an approach to modelling an ORF that roughly describes the typical shape and photosensitivity of a cell, but that can flexibly adapt variation across neurons.
To simultaneously account for these factors, we a novel variant of the GP-Bernoulli model. Let Yns E {0, 1} denote the response of neuron on trial to a holographic stimulus X; 1x3 As our approach designed for all-optical experiments, the response stimulation is assumed to be observed for all neurons at once using calcium or voltage imaging. Each stimulus delivers two-photon excitation to J 2 1 different target locations. with x/ = (cl,, S I/) E R representing the two-dimensional coordinates and laser power of the jth target (though note that we can also handle three-dimensional hologram coordinates without any substantial changes to the model). While some cases the targets will far enough apart space to not interact with each other, frequently a neuron will integrate two-photon excitation from multiple nearby targets, increasing the risk of OTS. To account for this, we model the probability of evoking a spike in neuron n by summing across all J points that neuron's 's ORF,
Ont Bernoulli - 9,)), Yn (x,) ) = Σ In (x/),
where is the logistic sigmoid function and 0,, IR is scalar threshold. We model each ORF using a three-dimensional GP prior that describes the effect of stimulating at a given location and power: In (P(m,,,) (:), k(-, -)), where m,, and & are respectively the mean and covariance functions of the GP. During inference (as discussed below) we constrain 9n. to be non-negative, which ensures that stimulation cannot have an inhibitory effect. Further, using non-negativity constraint rather than (e.g.) applying a rectifying nonlinearity to In guarantees that the posterior distribution remains log-concave, facilitating tractable identification of the global optimum.
To encode prior knowledge about optogenetic stimulation, the mean function should have the property that the probability of evoking spike increases with laser power but decays quickly the hologram
3

---

is positioned further away from the nucleus. To this end, we set
my (x) = pl -x /o
-, (c1, 02) denotes the coordinates of
(2)
for = L,.... N (note that we suppress reference to the target dimension j for notational clarity here and in the definition of the covariance function below). Here the holographic target, L, is the location of neuron n. and can be loosely interpreted as the average opsin expression level, determining how laser power affects the probability of spiking. While treat p and o as static hyperparameters, in future work one could consider a hierarchical model where the mean function is also learned, with its accuracy improving additional ORFs are probed.
We model the covariance between points on the ORF using the radial basis function (RBF) kernel k; which regularises the ORF by encouraging it to vary smoothly through space and power. However, we note that other covariance kernels could be used here provided that they are differentiable with respect x. We define
x(x,x') = a² exp (-x -oxx x)) + o2
(3)
where = diag(As ()s, A. Ar) is diagonal matrix with each diagonal term giving the characteristic lengthscale for the corresponding GP dimension. We assume that the lengthscales X, in the two spatial dimensions are equal (a typical assumption in holography, see e.g. ref. [6], but take a different lengthscale x for the laser power dimension. Finally, we denote the hyperparameters of the ORF prior by o = (a², oy As, Al, P. o.)
3.2 Inference
During ORF mapping experiment we first collect calibration data by performing stimulation in series of trials X1,..., Xy,, where each X₃ targets J different locations and laser powers near random set of neurons. Then we combine these stimuli with the corresponding observations of any evoked activity N.,....,,3 Y,, to estimate the ORFs. Here each y.,. E 13, describes the population response to the stimulus X,. Due to the non-negativity constraints p(gn Yn, {x, ,/= On, @) is not itself GP, and hence we perform maximum a posteriori (MAP) inference rather than attempting to obtain complete description of posterior uncertainty (though see the supplementary material for full treatment of posterior uncertainty in the = I case). To this end, the MAP estimates are obtained by maximising the log-posterior, which, via Bayes rule, is equivalent to computing
Jn. the posterior
gn. 8, argmax On
Σ In p(yne
Xe. 0%) + In (xI), In (x(x)) (0)
(4)
such that 9n (xx) 2 for I....,
Note that we have abbreviated 9. (X) (the MAP version of the ORF) as In- To solve Equation 4 we alternate between using Newton's 's method with log-barrier update 9n (such that the optimised ORFs respect the non-negativity constraints), and performing gradient descent steps for updating O,, We also adaptively set the Newton step-size using a standard backtracking line-search method.
3.3 Optimisation approach
Given the estimated ORFs, we then aim to identify holographic stimuli that evoke specific neural activity patterns as accurately possible (i.e., minimising off- off-target activation). To do this, introduce a novel optimisation approach that exploits properties of the GP to compute gradients of an eтror function between target activity pattern and the (predicted) evoked activity. To this end, let (ân, 0n I be the MAP-estimated ORFs, and define the predicted evoked activity as
(x. 9) (o((1 ".).....o. (0....,,(xx(( 0N)) RR,,
(5)
where (x) = E/1 = Hn (x'). Let Ω E {0, 1. denote a target activity pattern. The optimal holographic stimulus is then
Xaptimal argmin || () such that 0 I < I
(6)

---

Note that we include Amax as an upper bound on the laser power, which can be set to match the power practically deliverable by the microscopy system or to prevent tissue damage due to excess heating.
We propose to solve Equation Juuinng a projected gradient descent algorithm. While for the ORF mapping phase only need compute the ORFs In small set of stimulation points x,, for the optimisation phase we must evaluate the gradients of of[Euutio 6 (and therefore of the ORFs) at series of new, unobserved points constituting an optimisation path. However, is immediately obvious how to evaluate these gradients for a general nonparametrie function In- The key idea for our approach that each step, gradient descent will provide an updated estimate x" of the optimal stimulus, which will act as a "test point" commonly used in GP regression. We then leverage properties of the GP to perform inference of the ORF gradients at x*, and subsequently perform gradient descent update that reduces the value of the objective function in [lEEuuatio 6
Inferring the ORF gradients relies the fact that GP and its derivative are jointly GP-distributed, and hence the derivative can be inferred just from small number observations the GP itself. To this end, note that the covariance between GP and its derivative is given by [45] See 9.4]
0k(x,, x*)
COv (xx)) Duy Dra (x(x))
(7)
where d is any dimension of x" and where we have suppressed reference to the hologram target for notational clarity here and in the following equation. The expression in Eqqation 7 can be used to obtain a predictive distribution over derivative functions consistent with the observed data points, which we to evaluate the derivatives of 9n the test point x" via
(x*) (x') + cov (2). Ox²
0gn(x*)) Ox²s
(x) mx,
(8)
Or²
where X' E IB-Tx3 is the collection of JT unique points on the ORF that were probed during ORF mapping. 9r. (X), mn (X)= RJ give the value of the GP and mean function evaluated at each such point, and where K is the corresponding GP covariance matrix obtained by evaluating the covariance kernel & at every pair of elements of X. Further details are provided in the supplementary material.
Equation 8 alllows us to define closed-form gradient of the objective function, which we use in a projected gradient descent algorithm (Algorithm Tp. Note, however, that the objective function is not convex in x, and therefore such updates are only guaranteed to converge to a locally optimal solution. We therefore typically run Algorithm with several random initialisations and select the optimised target with the lowest predicted error.
Algorithm 1: Bayesian target optimisation (Bataro).
Compute MAP estimates of ORFs {ân, 0, 1N from calibration data {ym (x1 K using Newton's 's method with log- lor-berrier.
2 Initialise targets x E R. p x3 to random locations near the somas the target neurons and with random laser powers.
3 while target not converged do Construct gradient vectors (x) for 1...., N using inference of ORF derivatives (Equation 8).
n=1 Set 6- = -2 (2,, - ofn (x) (n))o'(An, - V (x). Perform gradient descent update x x x 55xx with step-size 3. Project laser power onto feasible domain, I, +- min(I,, I for 1, I....,J.
8 end
4
Results
4.1
Simulations
We first tested Bataro by simulating an optical "write-in" experiment, where we attempted to write specific activity patterns into a hypothetical neural population. To do this, we positioned 50 neurons randomly in 250 pm 250 μm field of view (FOV), sampled each neuron's 's ORF from its GP
5

---

Nuclear stimulation
b
Optimised stimulation
Spika probabiity 1.0 0.6
0.1
200
Power (mW)
Y (um)
150
60
50
250
500
Grad descent iters
50
8
random 6-target ensembles
Write- error
50
100
150
200
250
100
150
200
250
distance (um)
distance (um)
Spike probability
Write-in
Write-in error 1.61
2
1.0
0.5
Nuclear
Optimised
A
Target nourons
Target nourons
Figure 2: Minimising off- stimulation using Bayesian target optimisation. (a) Direct nuclear stimulation 70 mW successfully activates the target neurons with high probability, but also activates nearby non-target neurons due to OTS. Triangles indicate target neurons. Shading indicates probability of spiking. Optical write-in error (bottom) given as the sum of squared errors between the evoked and desired activity patterns. (b) Optimised stimulation using Algorithm I repositions the hologram locations away from the nuclei non-target neurons, resulting in a substantial reduction off-target activation. (c) Optimisation trajectory of the 6 different target laser powers. Initial laser powers selected randomly between 50 and 70 mW. (d) Optimising holographic targets for 20 diffcrent random ensembles shows robust reduction of optical write-in error (average reduction, 75%).
prior, and generated responses to ensemble stimulation. To quantify the optical write-in error, we used the sum of squared errors Ex,,,, Yn (x))², where (x) the ground truth probability of neuron spiking in response to stimulus x. This error takes the interpretable value of ~-m if m non-target neurons are inadvertently recruited due to OTS. In this simulation, holographic stimulation of an example ensemble consisting of 6 selected neurons resulted in a substantial recruitment of activity from non-target neurons due to OTS, especially when multiple holograms converged on nearby non-target neurons Figure 2); optical write-in error 7.84).
To optimally reposition the holographic targets. we probed the population's ORFs using 10-target ensemble stimulation, and then used Bataro (Algorithm I.. This successfully repositioned holographic targets an appropriate distance from the nuclei of the target neurons, such that while minimal excitation was applied to the non-target neurons, the target neurons were still driven to spike with high probability (Figure 2b,c; write-in error = 1.61). To confirm that the improvement write-in precision was robust, we repeated the optimisation for 20 additional random 6-target ensembles and found that off-targed effects were substantially reduced in every case, with the write-in error reducing by 75% on average compared to naively stimulating the nucleus of each neuron (Figure 2t.).
Next, we investigated two variables known critically constrain the ability evoke specific neural activity patterns: (1) the density of ppsin-expressing neurons in the stimulation FOV, and (2) the size of the ensemble that one is attempting activate. Figure 3/ illustrates how increasing the number of apsin-expressing be inadvertently activated. This causes errors induced by direct nuclear stimulation of (e.g.) 5-neuron ensembles to grow rapidly (Figure 3p, grey curve), though note that the magnitude of the error depends on many experimental factors including microscope, and the soma-targeting efficacy. By optimising holographic target locations and laser powers, the average write-in error was reduced by 76% relative to nuclear stimulation (Figure 3p, red curve). Similarly, in these simulations, nuclear stimulation resulted in error that grew almost
neurons introduces many more potential non-target neurons that could
size of the ORF, numerical aperture of

---

neurons, error-2.88
nourons,
250
16
Nuclear
16-
14
Y distanco (jmm)
Write-in erro
12
Write-in
10
8
6
Space pionablty
Optimised
4
2
125
250
125
250
10
X distance (jm)
distance (μm)
Population size
Stimulated ensemble size
Figure 3: Performance of Bayesian target optimisation in increasingly difficult contexts. (a) Two example scenarios with low (left) and high (right) density opsin-expressing neurons. Triangles indicate example 5-target ensemble to be stimulated. At low density, the risk of OTS is low because neurons are often spaced far apart. However, high density, OTS arising from direct nuclear stimulation at high power is unavoidable. (b) Optical write-in error for nuclear and optimised stimulation of 5-target ensembles. Error bars show the mean error s.d. over 10 different simulations. For each simulation averaged the write-in error over 20 random ensembles. (c) Same as (b), but for fixed population size of 50 and varying ensemble size.
400
40
1.2
Laser power (mW)
0.9
ddistance (um)
300
0000
Write-in
20
200
100
30 jm
Intial
Optimised
Nuclear
Optimised
100
200
400
distance (uim)
Figure 4: Performance of Bayesian target optimisation using simulations based two-photon holographic optogenetics data. (a) Left: field view from a real in vitro optogenetics experiment showing every optimised single-target hologram. Opsin fused to red fluorescent protein mRuby3 so that oxssnnexpr neurons can be visualised. Unfilled white circles represent putative neurons detected by automated cell segmentation method. Smaller filled white circles represent optimised holographic targets, shown in relation to the cell nucleus by a straight white line. Right: zoomed view of optimised targets, corresponding to dashed region the left panel. (b) Optimised laser powers relative to their initialised values show how Bataro exploits differences in photosensitivity to avoid OTS. Each circle represents the power delivered single holographic target. (c) Reduction of optical write-in error using optimised holographic targets. Each circle represents the error when attempting to stimulate single neuron.
directly proportionally to the size the stimulated ensemble, indicating that almost as many non- target neurons activated as target neurons. However, Bataro reduced the average write-in error by 69% across ensembles size 1-15 (Figure 3b).
Application to holographic optogenetics experiments
To test whether Bataro could eliminate OTS in realistic settings, we created synthetic optogenetics experiments involving 100 neurons from a small number (n=4) of detailed cell-attached recordings Briefly, each slice experiment was performed by first establishing cell-attached clectrophysiological recording of a single 1.2/3 interneuron from mouse V1 expressing the soma- targeted opsin ChroME2f. Then, the recorded neuron was optogenetically stimulated at dense grid
in slice 42

---

75 jm
Lm
50 JM
Z- Z-projection
50
Any
65
harget
neuron
jocloop
25 jim
neuron
pm
Axy um pm
jm
60
Power (mV)
35
Writo in arror
?
30
80 00 000 00
Depth
25
Nuclear (2 03)
Nuclear (1.71)
20
Spite grobil
-20
Nuolear
Optimined
Inaial
Optimised
Nucioa
Optimised
Optimised (x0.01)
Optimised (0.25)
target
25
50
25
50
Neuron
Neuron
Figure 5: Optimisation of holographic targets in threc-dimensional space. (a) Example target neuron and optimised holographic stimulus (plane 3, 25 μm). By repositioning the hologram (primarily in the x/y dimensions), off-target activation entirely eliminated (bottom spike probability plots, inset numbers show write-in error). Deepest plane labelled as 0 yum by convention. Solid white circle indicates target neuron. Dashed white circles indicate nearby non-targed neurons that must be avoided. (b) Similar to (a), but for neuron in plane 2 (50 jm) and with repositioning of the hologram primarily in the z dimension. (c) All optimised single-target holograms over four stimulation planes. Targets shown their original depth for visualisation. Average displacement of optimised holographic targets relative to nuclei, 8.3 Jum. Scale bar, 20 jm. (d) Max-projection over z-planes simultaneously showing locations of all segmented opsin-exp neurons. (e) Optimised depths of holographic targets corresponding to neurons in (c). Depths are shown in comparison to one of four depths that the target neuron was segmented during the experiment. Optimised laser powers relative their initialised values. (g) Reduction of optical write-in error for all neurons, across multiple depths.
of locations surrounding the cell and at range of laser powers to comprehensively map its ORF (see Figure Ip. "100% sampled" column, also see supplementary material for additional examples). We used this data create a set of four "ground truth" ORFs by fitting the EP--Bennuulli model from Equation I where the hyperparameters of the GP covariance kernel were selected using the cross- validated predictive likelihoodd. Next, we used fluorescence image from separate experiment to extract the locations of 16 putative opsin-expressing neurons, and randomly assigned each putative neuron one of four ground-truth ORFs. Together, this enabled us to simulate responses to optogenetic stimulation arbitrary locations and laser powers and for a large number of neurons that were realistically distributed in space.
We sampled responses to optogenetic stimulation at sparse grid of locations surrounding each neuron to map their ORFs hypothetical experiment (see supplementary methods for details). Then, we used Bataro to compute single-target stimuli that activated each neuron individually and minimised OTS (Figure 4b). We found that, the one hand, subset of neurons required no change from nuclear stimulation at high power due to being spatially isolated. However, regions with large numbers of closely-packed, opsin-expressing neurons required repositioning of the holographic targets (Figure 4p, right) and adjustment to the laser powers (Figure 4b) eliminate OTS (average distance between nuclear and optimised stimuli, 4.9 jm). Across all single-target holograms, target optimisation reduced the average write-in error by 85% (Figure 4c). We also performed control

---

analysis by matching the nuclear stimulation laser power to the average optimised laser power and obtained similar results (Figure S3).
Some holographic microscopy systems cannot position holograms at arbitrary continuous depths, instead either requiring holograms to be positioned in two dimensions on preselected number of discrete stimulation planes. Hence, thus far, we focused our efforts on the much more common case of two-dimensional optogeneties experiments. However, OTS is generally observed be even more prevalent in three dimensions, as the optical point spread function is elongated in the "z"-axis 1 17.. We therefore tested Bataro in the more general setting of three spatial dimensions and laser power (i.e., four dimensions in total), under the assumption that some microscopy systems may be able to flexibly reposition holograms continuously across depths.
15
We used the fact that the cell-attached recordings were repeated at multiple depths to generate synthetic optogenetics experiments with four-dimensional ORFs (Figure 5p.b). Optimisation of holographic targets successfully reduced OTS by exploiting both the x/y dimensions (Figure 5p), as well the z dimension (Figure 5p), in addition automatically adjusting the laser power where needed, Repeating the optimisation for every single-target hologram reduced the average optical write-in error by 89% (from 0.7 0.08. Figure 5p-g; average distance between nuclear and optimised targets, 8.3 jm), demonstrating the feasibility of largely eliminating OTS in optogenetics experiments involving all three spatial dimensions.
5 Conclusions
We developed Bataro, a novel computational framework for two-photon optogenetics that optimises holographic stimuli to evoke neural activity patterns with minimal off-target activation. Our pre- liminary in vitro experiments predict that specific adjustments to holographic target placements of 5-10 μam (on average) could substantially reduce OTS, even in the more challenging regime of three-dimensional optogenetic stimulation. Our key idea was to simultaneously exploit variability in photoexcitabilityy residual opsin expression in proximal dendrites (due to imperfections in opsin toraatarge,ing), and the specific spatial arrangement of neurons in the stimulation FOV. Optimisation using these three factors ultimately led to reductions in the optical write-in error of 85-89% for single-target holograms using data from in vitro two-photon optogenetics experiments.
Whether achieving such stimulus precision is worth spending valuable experiment time mapping ORFs is, of course, highly dependent on the application. While for some experiments minimising off-target activation may not be critical, there are many cases where stimulating with as close to "true" single-cell precision as possible is essential to testing a scientific hypothesis. For example, "hub" neurons in the hippocampus are believed to orchestrate network-level function through their unilateral activity [46] 47 Thus, using two-photon optogenetics to test whether an individual neuron represents functional hub by stimulating that neuron alone requires true single-cell precision. In such cases, we believe devoting experiment time to ORF mapping and stimulus optimisation is a necessary trade-off.
Throughout our analysis we have assumed that ORFs are stable over time (i.e., are stationary). However, desensitisation of opsins following prolonged illumination has previously been observed M. indicating that some nonstationarity could be encountered in real experiments. To account for this, we note that after the ORFs have been mapped and holographic stimuli optimised, it is straightforward to periodically recompute the MAP estimate of the ORF as the experiment proceeds. This would provide an updated estimate of the laser power required to maintain optimal precision. We hypothesise that this would not require repeating the entire ORF mapping phase as we do not expect the shape of the ORF to change drastically beyond desensitisation.
Standard GP regression methods scale cubically with the number of locations and laser powers probed, and therefore estimation of each neuron's 's ORF could encroach on experiment time if ORFs are mapped at high detail just due to required computation time. However, computational efficiency could be improved by appealing to modern GP techniques such as inducing point methods 448) or by adaptively mapping ORFs using Bayesian active learning to minimise the number probed locations and powers (4 50 While in this work we have focused primarily on demonstrating the feasibility of overcoming OTS, in future work we plan to minimise both required experiment and computation time in order to deploy Bataro online.

---

Acknowledgements
We thank Darcy Peterka, Benjamin Antin, Cole Hurwitz, Shizhe Chen, and Ben Shababo for helpful discussions and suggestions. This work was funded by NIH awards 1RF1MH120680 and IU19NS10761344 to НA and LP. MAT and LP were supported by the Gatsby Charitable Foundation and NSF NeuroNex award 1707398.
References
[1] John Peter Rickgauer and David W Tank. Two-photon excitation of channelrhodopsinn2 saturation. Proceedings of the National Academy of Sciences, 106(35): 1503,, 2009.
[2] Eirini Papagiakoumou, Francesca Anselmi, Aurélien Bègue, Vincent De Sars, Jesper Glückstad, Ehud Y Isacoff, and Valentina Emiliani. Scanless two-photon excitation of channeirhoddopsin-2 Nature Methods, 7(10):4488855 2010.
[3] Adam M Packer, Darcy S Peterka, Jan Hirtz, Rohit Prakash, Karl Deisseroth, and Rafael Yuste. Two-photon optogenetics of dendritic spines and neural eircuits. Nature Methods, 9(22):1222-1255, 2012.
[4] Rohit Prakash, Ofer Yizhar, Benjamin Grewe, Charu Ramakrishnan, Nancy Wang, Inbal Goshen, Adam M Packer, Darcy S Peterka, Rafael Yuste, Mark J Schnitzer, et al. Two-photon optogenetic toolbox for fast inhibition, excitation and bistable modulation. Nature Methods, 9(12):11111119, 2012.
[5] Oscar Hernandez, Eirini Papagiakoumou, Dimitrii Tanese, Kevin Fidelin, Claire Wyart, and Valentina Emiliani. Three-dimensional spatiotemporal focusing of holographic patterns. Nature Communications, 7(1):1-11, 2016.
[6] Nicolas C Pégard, Alan R Mardinly, Ian Antón Oldenburg. Savitha Sridharan, Laura Waller, and Hillel Adesnik. dimee-dimeniional scanless holographic optogenetics with temporal focusing (3d-shot). Nature Communications, 8(1):1-14, 2017.
[7] Emiliano Ronzitti, Cathie Ventalon, Marco Canepari, Benoit Forget, Eirini Papagiakoumou, and Valentina Emiliani. Recent advances in patterned photostimulation for optogenetics. Journal of Optics, 19(1 1):1 13001, 2017.
[8] Nicolò Accanto, Clément Molinier, Dimitrii Tanese, Emiliano Ronzitti, Zachary L Newman, Claire Wyart, Ehud Isacoff, Eirini Papagiakoumou, and Valentina Emiliani. Multiplexed temporally focused light shaping for high-resolution multi-cell targeting. Optica, 5(11):1478- 1491, 2018
[9] Weijian Yang, Luis Carrilo--Reid Yuki Bando, Darcy S Peterka, and Rafael Yuste. Simultaneous two-photon imaging and two-photon optogeneties of cortical circuits in three dimensions. elife, 7:e32671, 2018
[10] Alan R Mardinly, Ian Antón Oldenburg, Nicolas C Pégard, Savitha Sridharan, Evan H Lyall, Kirill Chesnov, Stephen G Brohawn, Laura Waller, and Hillel Adesnik. Precise multimodal optical control of neural ensemble activity. Nature Neuroscience, 21(6)8881-8) 18933, 2018.
[11] 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)eeaww5200 2019.
[12] Luis Carrillo-Rede Shuting Han, Weijian Yang, Alejandro Akrouh, and Rafael Yuste. Control- ling visually guided behavior by holographic recalling of cortical ensembles. Cell, 178(2)447- 457, 2019.
[13) Hillel Adesnik and Lamiae Abdeladim. Probing neural codes with two-photon holographic optogenetics. Nature Neuroscience, 24(10):1356-136) 2021.
10

---

[14] Eirini Papagiakoumou, Emiliano Ronzitti, and Valentina Emiliani. Scanless two-photon excita- tion with temporal focusing. Nature Methods, 17(6):571-581, 2020.
[15] Christopher A Baker, Yishai M Elyada, Andres Parra, and M McLean Bolton. Cellular resolution circuit mapping with temporal-focused excitation of soma-targeted channelrhodopsin. Elife, 3:e14, 2016.
[16] Or A Shemesh, Dimitrii Tanese, Valeria Zampini, Changyang Linghu, Kiryl Piatkevich, Emil- iano Ronzitti, Eirini Papagiakoumou, Edward S Boyden, and Valentina Emiliani. Temporally precise single-cell-resolutio optogenetics. Nature Neuroscience, 10612):176--1806 2017.
[17] Travis A Hage, Alice Bosma-Moody. Christopher A Baker, Megan B Kratz, Luke Campagnola, Tim Jarsky, Hongkui Zeng, and Gabe J Murphy, Synaptic connectivity to 12/3 of primary visual cortex measured by two-photon optogenetic stimulation. Elife, 11:e7 103, 2022.
[18] Valentina Emiliani, Adam Cohen, Karl Deisseroth, and Michael Häusser. All-optical interro- gation of neural circuits. Journal of Neuroscience, 35(41)::3917-
2015.
[19] Adam Packer, Lloyd E Russell, Henry WP Dalgleish, and Michael Häusser. Simultancous all-optical manipulation and recording of neural circuit activity with cellular resolution in vivo. Nature Methods, 12(2): 140-146, 2015.
[20] Zihui Zhang, Lloyd E Russell, Adam M Packer, Oliver M Gauld, and Michael Häusser. Closed- loop all-optical interrogation neural circuits in vivo. Nature Methods, 15(12)10077-1000 2018.
[21] Lloyd E Russell, Henry WP Dalgleish, Rebecca Nuthrown, Oliver M Gauld, Dustin Herrmann, Mehmet Fişek, Adam M Packer, and Michael Häusser. All-optical interrogation of neural circuits in behaving mice. Nature Protocols, 17(7):1579-1620 2022.
[22] Hayley A Bounds, Masato Sadahiro, William D Hendricks, Marta Gajowa, Ian Antón Oldenburg, Karthika Gopakumar, Daniel Quintana, Tanya Daigle, Hongkui Zeng, and Hillel Adesnik. Multifunctional cre-dependent transgenic mice for high-precision all-optical interrogation of neural circuits. bioRxiv, page 463223, 2021.
[23] Kamiar Rahnama Rad and Liam Paninski. Efficient, adaptive estimation of two-dimensional firing rate surfaces via gaussian process methods. Network: Computation in Neural Systems, 21(3-4):142 168, 2010.
[24] Cristina Savin and Gasper Tkacik. Estimating nonlinear neural response functions using gp priors and kronecker methods. Advances in Neural Information Processing Systems, 29, 2016.
[25] M Rule, P Chaudhuri- Chaudhuri-Vayallmbreee Vayalambrone, Marino Krstulovic, Marius Bauza, Julija Krupic, and Timothy O'Leary. Variational gaussiann point-process methods for grid cells. bioRxiv, pages 2023 03, 2023.
[26] Jakob H Macke, Sebastian Gerwinn, Leonard E White, Matthias Kaschube, and Matthias Bethge. Gaussian process methods for estimating cortical maps. Neurolmage, (6(2))570-588) 2011.
[27] Mijung Park, Greg Horwitz, and Jonathan Pillow. Active learning of neural response functions with gaussian processes. Advances in Neural Information Processing Systems, 24, 2011.
[28] Mijung Park, J Patrick Weller. Gregory Horwitz, and Jonathan W Pillow. Bayesian ac- tive learning of neural firing rate maps with transformed gaussian process priors. Neural Computation, 26(8): 2014.
[29] Jonathan P Newman, Ming-fai Fong, Daniel C Millard, Clarissa J Whitmire, Garrett B Stanley, and Steve M Potter. Optogenetic feedback control of neural activity. Elife, 4:e07 2015.
[30] Yuxiao Yang, Allison T Connolly, and Maryam M Shanechi. A control-theoretic system identification framework and a real-time closed-loop clinical simulation testbed for electrical brain stimulation. Journal of Neural Engineering, 15(6):066077, 2018,

---

[31] Michael F Bolus, Adam A Willats, Clarissa J Whitmire, Christopher J Rozell, and Garrett B Stanley. Design strategies for dynamic closed-loop optogenetic neurocontrol in vivo. Journal of Neural Engineering, 15(2):026011 2018.
[32] 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.
[33] Michael Bolus, Adam A Willats, Christopher J Rozell, and Garrett B Stanley. State-space op- timal feedback control of optogenetically driven neural activity. Journal of Neural Engineering. 18(3):036006, 2021.
[34] 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)3224 345, 2021.
[35] Yashar Ahmadian, Adam M Packer, Rafael Yuste, and Liam Paninski. Designing optimal stimuli to control neuronal spike timing. Journal of Neurophysiology, 106(2):1038--003, 2011.
[36] Tao Hu and Dmitri Chklovskii. Reconstruction of sparse circuits using multi-neuronal excitation (rescume). Advances in Neural Information Processing Systems, 22:790-798, 2009.
[37] Alyson K Fletcher, Sundeep Rangan, Lav R Varshney, and Aniruddha Bhargava. Neural reconstruction with approximate message passing (neuramp). Advances in Neural Information Processing Systems, 2011.
[38] Ben Shababo, Brooks Paige, Ari Pakman, and Liam Paninski. Bayesian inference and on- line experimental design for mapping neural microcircuits. Advances in Neural Information Processing Systems, 26: 134-11233,, 2013.
[39] Alyson K Fletcher and Sundeep Rangan. Scalable inference for neuronal connectivity from calcium imaging. Advances in Neural Information Processing Systems, 2014.
[40] Laurence Aitchison, Lloyd Russell, Adam Packer, Jinyao Yan, Philippe Castonguay, Michael Häusser, 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, 2017.
[41] Anne Draelos and John Pearson. Online neural connectivity estimation with noisy group testing. Advances in Neural Information Processing Systems, 33, 2020.
[42] Marcus A Triplett, Marta Gajowa, Benjamin Antin, Masato Sadahiro, Hillel Adesnik, and Liam Paninski. Rapid learning of neural circuitry from holographic ensemble stimulation enabled by basel-based compressed sensing. bioRxiv, page 507926, 2022.
[43] Yoav Printz, Pritish Patil, Mathias Mahn, Asaf Benjamin, Anna Litvin, Rivka Levy, Max Bringmann, and Ofer Yizhar. Determinants of functional synaptic connectivity among amygdala- projecting prefrontal cortical neurons in male mice. Nature Communications, 14(1):1667, 2023.
[44] Savitha Sridharan, Marta A Gajowa, Mora B Ogando, Uday K Jagadisan, Lamiae Abdeladim. Masato Sadahiro, Hayley A Bounds, William D Hendricks, Toby S Turney, Ian Tayler, et al. High-performance microbial opsins for spatially and temporally precise perturbations of large neuronal networks. Neuron, 2022.
[45] Christopher KI Williams and Carl Edward Rasmussen. Gaussian Processes for Machine Learning, volume 2. MIT Press Cambridge, MA, 2006.
[46] Paolo Bonifazi, Miri Goldin, Michel A Picardo, Isabel Jorquera, A Cattani, Gregory Bianconi, Alfonso Represa, Yehezkel Ben-Ari, and Rosa Cossart. Gabaergic hub neurons orchestrate synchrony in developing hippocampal networks. Science, 126(5958):4419-1444 2009.
[47] Marco Bocchio, Claire Gouny, David Angulo-Garcia, Tom Toulat, Thomas Tressard, Eleonora Quiroli, Agnès Baude, and Rosa Cossart. Hippocampal hub neurons maintain distinct connec- tivity throughout their lifetime. Nature Communications, 11(1):4559, 2020.
12

---

[48] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial Intelligence and Statisties (AISTATS), pages 567-574. PMLR, 2009.
[49] Jeremy Lewi, Robert Butera, and Liam Paninski. Sequential optimal design of neurophysiology experiments. Neural Computation, 21(3):19-687, 2009.
[50] Jonathan W Pillow and Mijung Park. Adaptive bayesian methods for closed-loop neurophysiol- ogy. Closed Loop Newroscience, pages 3-18, 2016.
[51] Andrew James McHutchon et al. Nonlinear modelling and control using Gaussian processes. PhD thesis, Citeseer, 2015.
13

---

Supplementary material
6.1 Animal ethics statement
All experiments on animals were conducted with approval of the Animal Care and Use Committee of the University of California, Berkeley.
6.2 Compute
All computational procedures were performed either desktop workstation running Ubuntu 18.04 with an Intel Xeon E5-2620 v4 CPU, four GTX 1080 Ti GPUs, and 112GB RAM, or on the Axon computer cluster based at the Zuckerman Institute (Columbia University) using nodes comprised of two Xeon E5-2660 v4 CPUs, eight GTX 1080 Ti GPUs, and 125GB RAM.
6.3 Broader societal impact
Our work is significant for interventional approaches to studying the brain and its connection to disease. By minimising off-target activation, Bayesian target optimisation could enable (e.g.) more precise synaptic connectivity mapping, improving our understanding of neural circuitry. This advancement has potential implications for understanding brain disorders like epilepsy, where abnormal synaptic connections are central to seizure generation and propagation. Deepening our understanding of these diseases can lead to enhanced targeted interventions and more effective therapeutic strategies, benefiting individuals with neurological disorders.
6.4 Code availability
An open-source implementattinnnnnnnnnnnn of Bayesian target optimisation is available in Python at https htt.abbb.ccmmmmrruuutpppett/bbba
6.5
Single-target holographic stimulus optimisation with posterior uncertainty
Here we provide further mathematical details for optimising holographic stimuli. First, we develop the approach for single optogenetic targets, as this is most closely related to existing GP-based receptive field inference techniques. The single-target case also allows us to have a full treatment of posterior uncertainty (unlike for optimising ensemble stimuli) which may be desired in certain applications.
Optogenetic receptive field model. We use GP-Bernoulli approach to model the response you of neuron n on trial a single-target stimulus X1,
Yat (emnouli(o(g,.11555
where the stimulus x, = (Cie, C2t, I E I represents the two-dimensional coordinates and laser power of the I-th hologram. Each ORF follows three-dimensional GP prior Jn. ~ GPPPmm (m, where m, and i again are the mean and covariance functions of the GP.
Posterior inference. Unlike for ensemble stimulation, for single-target stimulation do not require that the ORF 9n is non-negative. This is because now a point on the ORF becomes inhibitory (by taking a negative value). it will not conflict with excitation from any other hologram. Consequently, the posterior of 9n. is a GP, which allows us to work with a full description of posterior uncertainty. To compute the posterior, approximating the posterior using a multivariate normal g(gn Mn: E.) Normal pm, Σn) RV plyn Yn, X, Φ)). The mean pnn is obtained by maximising the log-posterior, given by the expression
the conventional Laplace approximation. Briefly, this consists of
T
Imp(gn Yn, X, Φ)
In p(y( xx, + In ((). 9n (xr) | () + const, Σ
(10)
where X = (x XT) and where const does not depend on 9n. Since the posterior is log-concave
Ga, we use Newton's method identify the global optimum of Equation 10] and adaptively set the Newton step-size using standard backtracking line-search method. Letting H VV g.. In p(g.
14

---

.n. X. σ) be the Hessian of the log-posterion the posterior covariance matrix is obtained by setting En
5n Pm
Target optimisation. Let G = (g1,..., .aal.....g....... gN). and define the predicted evoked activity for single holographic targets i(x, G) (o(g) (egg--)))( To minimise the error between target binary activity pattern $2 E {0, 1]. and the predicted evoked activity, we solve an optimisation problem that accounts for the uncertainty in the ORF estimates:
Xoptimal = argmin E [1 i(x, G]]|[] such that 0 < S
Iaxx
(11)
x
where g(G
- (9n Mu: En) gives the joint posterior across all ORFs. To solve we first sample ORFs gú (for = S)..., from their posterior distributions to approximate the expected error the current estimate x
p.. [1
Equation 11
1 S
[) 11-0((x,
G)||
3E #=1
On ((nb
(12)
Then, we compute the partial derivative (in dimension d) of the expected error by differentiating through the Monte Carlo approximation.
8 N >ΣΣ 8=1n=1
① Dzy E [| O(x' G)1|
S2n (9. (x) (ol)"(()))))
Oz -.! ( (x*).
~
(13)
Next must evaluate the partial derivative the right-hand side of Equuation 13 We use the fact that GP and its derivative are jointly GP-distributed, and hence infer the derivative from observations of the ORF The covariance between GP and its derivative is given by [45] Sec 9.4]
cov (B.(xx))
ak(x,, x*)
a²
| 2
x
In (xxx))
(a.t x5)
(14)
Ozyy
dzy
Xa
where the second equality specific to the RBF covariance. Thus, we can use Equation 14 obtain the posterior predictive mean for the derivative GPs in closed form as [5]] Sec 2.7]
[0g((x)) Oxy
dmn(x") cov (x), 0gn(x²) dr²
E Eg(g lnnn E.n)
K (Mn - mn (X)). (15)
Ōr²
Here X (x., ........ Xy. is the collection of unique points the ORF probed during calibration.
Equation 15 is combined with an expression for the posterior predictive variance, one obtains a full predictive distribution over derivative functions consistent with the observed neural responses. However, rather than working with this full distribution, we instead use Equation 15) to approximate the derivatives of the Monte Carlo samples by replacing the posterior mean f4, with Monte Carlo sample,
0gh (x*)
mx")
0gn(x²)) dr²
cov (x)),
K (gl -m, (X)).
(16)
Ox
dr
Equation 16 then allows to define closed-fome approximate gradient Vx gh" at test point x*, defined as
Oggh (4) Oxi
(x')
0gh"
(x*) Ox
Vx-
((()
(17)
which we use in the single-target projected gradient descent algorithm (Algorithm 2). Note that one could also consider quadrature approach to solving Eqqation 12, which may be more efficient than Monte Carlo sampling. However, the presentation of the Monte Carlo approach is instructive for deriving the optimisation of ensemble stimuli below.
15

---

Algorithm 2: Projected Monte Carlo gradient descent algorithm for optimising single-target holograms
Infer ORF posterior q(G μ.. Σ) from calibration data {y,, J Jn X using the Laplace
approximation.
Precompute the negative of the Hessian W.. VV pyy, X, gn)) Ia. -,, for each nt.
3 Initialise x to random location near soma of target neuron and with random laser power. while targer not converged do for 7 .,...,... .V do
Compute mean and variance of posterior predictive distribution at current target estimate
x via pm (x) = (x) + k(x, x) K ( m, (X)), and
(x) k(xx
le X
) (XX x).
Sample ORFs at the current target estimate, gn (x) Nu (x), (x)) for 1..... S.
Construct approximate gradients V xgnn for == 1...., S using Equation 17
(s)
end
(9,, (ggh (x))o'(g V (e) (x) per Equation 13 Perform gradient descent update, 35x with step-size B.
Set 8x
12
Project laser power onto feasible domain, t min (I, Imax).
13
end
6.6
Additional details on ensemble stimulus optimisation approach
The approach for optimising holographic ensemble stimuli is based the approach for single-target optimisation, but modified to account for differences in the ORF model and inference. In particular, we again seek to minimise the error between a target activity pattern S2 and the predicted evoked activity, but now using the MAP estimates G = {ūn; 0n ) in place of the full posterior distributions.
Let 9(x,9) = (a(9) (x) - i (ffN((( - 0x)) be the predicted population response to an ensemble stimulus, where in (x) Ej=1 In (x)). The optimal ensemble stimulus is now
Xoptimal = argmin || j(xx 9)1 argmin Σ
2n -(fm(x)) - 8n))
(18)
n=1
such that max Evaluating the partial derivative of Equation 18 with respect dimension d of test point x* yields,
0
N
1 0(x 9)12 -2 dry
Σ (on - o0(n(xx)) - (Sm(("))) - 0m) Oz
(x*).
(19)
n=1
The derivative the right-hand side of Equation 19 is given by a (x)) Σ/= J. 9n (x²),
J27
which requires computing the derivative of 9. (x²). ). To evaluate this derivative, to Equation 16 but substituting the MAP estimate in place of the posterior mean or Monte Carlo sample,
similar trick
0 dr
T
In (x') =
ar: mn (x") + cov (x((X), dxy "n((x))
K (g - m, (X)).
(20)
This expression can also arrived at by first evaluating the posterior predictive mean of gn (x*). and then differentiating with respect to z²
We use Equation 20 to define closed-form gradient Vx at test point x" via
(x"
29
Fn =
(21)
Oxt
Oxy
Finally, Equation 21 used in the projected gradient descent algorithm for optimising ensemble stimuli Algorithm
16

---

High coverage
Low coverage
b
18
C
1.2
20
20
16
Nuclear
Y (Ja))
Y distarcc (Jm)
Wiitl-in emor
14
12
Low coverage
High coverage
Wiife-inn arror
10-
R6
0,3
-10
-10
-20
-20
distance (um)
ddsstanc (um))
15
Nuciear
Optmieed
Stimulated ensemble size
Figure S1: Effect of reducing the number of points at which each ORF is probed. (a) In the high coverage case (left), each ORF is probed by stimulating at a 5xx5 5 grid of points near the soma (grid points separated by 10 ju), at three different laser powers. In the low coverage case (right), this reduces to stimulating just 3 grid (grid points separated by 12 ju) three powers. However, as the density of opsin-expressing neurons increases, ORFs are probed at high density even in the low coverage case as the grids from different neurons increasingly overlap. (b) Minimal performance difference between the high and low coverage cases in simulations with 50 neurons. (c) Reduction in optical write-in error using cell- attached recordings as in Figure 4 but with low coverage. Reduction in average write-in error, 74% (c.f. 85% with high coverage, Figure 4b).
6.7 Further details on simulations and "synthetic" optogenetics experiments
Simulations consisted of both ORF mapping and stimulus optimisation phases. ORF mapping required probing responses to stimulation at range of laser powers and stimulus locations. We defined a grid of stimulation points surrounding each neuron. the spatial dimensions, the grid ranged from -20 yum to 20 jam relative to the centroid of the neuron in steps of 10 μam, and powers ranged from 30 mW to 70 mW steps of 20 mW. The complete grid was thus given by the Cartesian product 202 10, 0, 10, 20} x { 20, 100, 0, 10, 20 {30, 50, 70). For opsin-expressing neurons that were spaced far apart, this coasse-resolution grid was sufficient because risk of OTS was low, and therefore ORF mapping
not needed at high detail. On the other hand, as the density of ppsin-expressing neurons increased, the grids surrounding each neuron increasingly overlapped with each other, resulting in much denser sampling of the ORFs.
For the synthetic optogeneties experiments (based on the cell-attached recordings), we used the same spatial grid spacing but used laser powers of 10, 25, and 40 mW to match the range of powers used in the underlying slice experiment, though note that the slice experiment had denser spacing than our chosen 15 mW (see example loose-patee patch recordings below), which we chose to reduce the ORF mapping time. For the optogenetics experiments involving three spatial dimensions, we extended the grid sampling to include depths of -60 μam to 60 um in steps of 30 μu. We also explored the effect of reducing the number of probed grid points to further reduce the time spent mapping ORFs, and found that Bayesian target optimisation maintained high performance when probing with 3 x 3 spatial grid of (-12, 0, 12} x 1-12, 12} (Figure SIp.
We selected the parameters of the GP covariance kernel using 5-fold cross-validation on separate set of recordings that were made on the same set of four cells, ensuring the hyperparameter selection was using out-of-sample data. Cross- Cross-validation validation was performed using grid search over set of possible hy perparaameers:: the possible radial lengthscales were 2, 4, 8, 16, the power lengthscales were 2, 4, 8, 16, and the amplitudes were 1, 2, 4, 8, 16. For each hy hyperarrameter combination 0 and for each cell, we used Newton's method to fit the GP-Bernoulli model to 80% of the trials in the loose- patch data, yielding ORF estimate ge. On the remaining 20% of the trials (denoted as Theld-out). we evaluated the log-likelinn likeeloodoo E EigTageop {yr In(ar (xx)) + (1 - ye) In - o (ğe (xx })) We averaged the likeliood across all five folds and across all four cells, and chose the hyperparameter combination& 9 that yielded the largest average log- likeliiooodd resulting in a radial lengthscale of 8, a power lengthscale of 16, and a kernel amplitude of 8.
The GP parameters for generating the simulations in Figure 3
inferring the resulting ORFs, and generating synthetic optogenetics experiments with two and three spatial dimensions are given in Table ST. For reference, a typical ORF mean function given in Figuue $2
17

---

GP
70 mw
Sample
Sample 2
Sample 3
Sample
25
-25
Y distance (Jm)
50 mW
25
0
-25
30 mW
1.0
25
goid Spiee
0.5
-25
0.0
-25
25
-25
25
distance (μm)
Figure S2: Example mean function (shown at three powers) used for simulations (left column). Also shown are four samples from the ORF prior corresponding to this mean function (right four columns). Parameters given in Table SI
Average performance (5 random repeats, matched powers)
1.5
40
000 00
1.2
Laser power (mW)
38
Write-in error
0.9
0
0.6
34
0.3
32
30
Nuclear
Optimised
Nuclear
Optimised
(average)
Figure S3: Performance of Bayesian target optimisation compared nuclear stimulation when average laser powers are matched. Each repeat consists of randomly reassigning neurons different optogenetic receptive fields while keeping their spatial positions fixed, and performing full mapping and target optimisation sequence. Within each repeat, nuclear stimulation was performed at the average optimised power. Left: average write-in error for nuclear stimulation remains substantially higher than with optimised stimulation despite having matched average powers. Each circle represents a single neuron. Right: powers used within each repetition. Note that these powers are lower than 40 m W used to benchmark performance in Figures 4 and 5 in the main text.
18

---

Parameter
Symbol
Value
Simulations (data generation)
Mean function excitability
p
0.125
Mean function width
3 x 10? pmm
Spike threshold
0
3.5
Kernel radial lengthscale
A
pm
Kernel power lengthscale
A
20 mW
Kernel amplitude
a²
0.2
Kernel marginal variance
o²
10-5
Simultaneously stimulated neurons during ORF mapping
J
Simulations (ORF inference)
Mean function excitability
P
0.125
Mean function width
3 x
102 pm
Kernel radial lengthscale
X
jm
Kernel power lengthscale
A
mW
Kernel amplitude
a²
Kernel marginal variance
o²
10-5
Learning rate for spike thresholds (fo, N.N
Number of random initialisations
Synthetic optogenetics experiments (two spatial dimensions)
Mean function excitability
P
0.175
Mean function width
x
10 μm
Kernel radial lengthscale
X
jum
Kernel power lengthscale
A
mW
Kernel amplitude
a²
Kernel marginal variance
o
10 10-5
Learning rate for spike thresholds ((On),) ).
Number of random initialisations
Synthetic optogenetics experiments (three spatial dimensions)
Mean function excitability
P
0.175
Mean function width (x/y dimensions)
22
x 102 Jm
Mean function width (z dimension)
3 3 103 pm
Kernel radial lengthscale (x/y dimensions)
X
pm
Kernel axial lengthscale (z dimension)
X
32 pm
Kernel power lengthscale
X
16 mW
Kernel amplitude
a
Kernel marginal variance Learning rate for spike thresholds ((O,) J.
o3
10 5
Number of random initialisations
Table S1: Parameters used for simulations and generating synthetic optogenetics experiments.
19

---

6.8 Additional examples of optogenetic receptive fields from cell-attached recordings
S7 show examples of four ORFs that have been comprehensively mapped using two- photon optogenetic stimulation and cell-attached recordings of evoked spikes. Note the unpredictable differences in ORF shape across laser powers and depths, motivating nonparametric approach.
Figures
S4
mW
mW
20 mw
30 mW
40 mW
Depth = 75 um
Y distance (um)
0000000000
000000000000
0000000000000
000000000
00000000000 0000000
25
00000000000000000000000000000000
-25
50 pm Y distance (um)
255
Y distance (um)
25
25 um
23
Y distance (um)
25
O pm
255
Y distance (um)
25
-25 um
25
Y distance (um)
25
-50 um
255
Y distance (um)
-55 um
-25
-25 O
25
-25
25
25 X distance (um)
-25
-25
0
25
-25 O
25
X distance (μm)
X distance (Hm)
X distance (um)
X distance (Hm)
Figure S4: Loose-patch recording and inferred ORF (experiment 1/4).
20

---

10 mW
15 mW
20 mW
30 mW
40 mW
Depth = 75 um Y distance (um)
000000
Y distance (um)
25
50 μm
92
Y distance (wm)
25 pm
25
Y distance (um)
O pm
255
Y distance (um)
25
-25 jum
25
Y distance (um)
25
-50 jm
-25
Y distance (um)
25
-75 um
25
-25
0
25
25
25
-25
25 distance (um)
25 X distance (um)
-25
-25 0 25
X distance (um)
X distance (Hm)
X distance (um)
Figure S5: Moose-patch recording and inferred ORF (experiment 2/4)
21

---

10 mW
15 mW
20 mW
30 mW
40 mW
Depth = 75 um Y distance (um)
000000
25
Y distance (um)
25
50 μm
Y distance (wm)
25 pm
25
Y distance (mm)
O pm
255
Y distance (um)
25
-25 jum
Y distance (um)
25
-50 jm
-25
Y distance (um)
25
-75 um
-25
0
25
25
25
-25
25 distance (um)
-25
25
-25
0
25
X distance (um)
X distance (Hm)
X distance (um)
X distance (um)
Figure S6: Loose-patch recording and inferred ORF (experiment 3/4)
22

---

10 mW
15 mW
20 mW
30 mW
40 mW
Depth = 75 um Y distance (um)
000000
Y distance (um)
25
50 μm
92
Y distance (wm)
25 pm
25
Y distance (um)
O μm
255
Y distance (um)
25
-25 jum
Y distance (um)
25
-50 jm
-25
Y distance (um)
25
-75 um
25
-25
0
25
25 X distance (Hm)
25
-25
25 distance (um)
-25
25
-25
0
25
X distance (um)
X distance (um)
X distance (um)
Figure S7:
Loose-patchs recording and inferred ORF (experiment 4/4)
23


Reviews:
[{"review_id": "rYzggRqOt5", "reviewer": "Reviewer_FU6t", "summary": "The authors present a set of methods to efficiently characterize the activation field under optogenetics and optimize the optogenetic stimulation patterns to target certain cells while avoiding the others. ", "questions": "• How is linear summation of activation from different optogenetic stimulation sites justified, especially if the sites are very close to each other?\n• Could a simple stimulation optimization approach, where the location is at the intersection of high g(x) for one cell and low g(x) for other cells, work?\n\nA characterization of the efficiency (reduction in the number of measurements for the same estimation quality) is missing.", "limitations": "Limitations have been adequately addressed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper tackles an important problem for reproducing the neural responses with high resolution in a neural circuit. The presented approach is direct, principled and effective.", "weaknesses": " The experimental validation of the method is limited. It is unclear how the real-data was used in analysis (Sec 4.2), and if it represents what would happen in a real experiment. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks again for spending the time to review our submission and for providing your comments.\n\n**The experimental validation of the method is limited. It is unclear how real-data was used in analysis (Sec 4.2) and if it represents what would happen in a real experiment.**\n\nThanks for noting that there is room to improve the clarity in the text about how the real data was used in the analysis. As we commented on in the global rebuttal, the ideal experimental data for validating these techniques are not currently available to us (i.e., an _in vivo_ demonstration of the technique, though we are actively collaborating to acquire such data). Therefore, we are limited to working with the data that we have, which come from slice experiments.\n\nWe made use of two kinds of experimental data: (1) detailed optogenetic receptive field (ORF) maps for single neurons, obtained by making a loose-patch recording of a single neuron, stimulating at regularly spaced locations surrounding the patched neuron and at multiple powers, and then fitting the GP-Bernoulli model; and (2) a fluorescence image of opsin expression from a separate experiment showing a typical distribution of opsin-expressing neurons in space. To use these data to simulate an experiment, we randomly assigned each opsin-expressing neuron in (2) an ORF from (1), and used these ORFs to sample responses to neural stimulation at arbitrary locations and powers.\n\n**How is linear summation from multiple nearby sites justified?**\n\nTo our knowledge, an experimental characterisation of the photocurrent evoked by multiple closely positioned holograms has not yet been performed. However, the photocurrent directly depends on the number of opsin molecules illuminated [4], and we do not expect an interaction effect between two-photon laser pulses that would nonlinearly recruit opsin molecules beyond more quickly approaching a saturation point (that is currently accounted for by the sigmoid nonlinearity). \n\nThat being said, some prior experiments using two-photon glutamate uncaging have found cases of nonlinear summation when stimulating multiple (_distal_) dendritic spines on the same branch [5]. However, in the experiments that our technique is designed for, opsin expression is restricted to the soma and (at most) the _proximal_ dendrites. We expect that instances of simultaneously stimulating multiple points along the same segment of proximal dendrite immediately adjacent to a target neuron's soma will be very rare, and we do not currently know whether this would measurably change the ability to optogenetically stimulate the neuron.\n\n**What about a simpler approach of finding high $g(x)$ for one cell and low $g(x)$ for other cells?**\n\nWe thank the reviewer for this intuitive suggestion, and note that for the single-target case this is similar to what the proposed gradient descent solution finds (though the supplementary algorithm for single-target optimisation further accounts for posterior uncertainty in the inferred ORFs). However, this does not generalise to the ensemble stimulation case because it does not account for the contribution from multiple holograms at once.\n\n**A characterisation of the efficiency is missing**\n\nWhile we show the effect of probing the ORFs with fewer stimuli (i.e. using downsampled grids) in the supplement, we will provide further data related to statistical efficiency in the final version of the manuscript and make an explicit reference to the result in the main text.\n\nThanks again for your comments!"}, {"review_id": "mEyo57Y5B9", "reviewer": "Reviewer_4tSm", "summary": "The authors present a method for reducing off-target stimulation of neurons during photostimulation by modifying the laser power and target locations. They use a Bayesian optimization approach to determine neuron responses (ORFs) to stimulations at different targets and laser powers and then choose the optimal target parameters to attain a desired neural response to stimulation. In contrast to typical stimulation targeted to the desired neuron’s soma, stimulating somewhat ‘off-target’ can still produce spikes in the desired neurons and avoid stimulating its nearby neighbors. ", "questions": "1. How well does this scale from an optimization (time) perspective? Can this be done for thousands of neurons? Ensembles of 50?\n\n2. How does the exact neuron expression matter (to the ORF shape)? Can a simplified version be used during the 2nd step of target optimization?\n\n3. Similarly, what do ORFs look like? Are there other references used for establishing the prior beyond the few shown here?\n\n4. Would it be possible to integrate spatial information (from e.g. calcium imaging) into the ORF prior in a useful way? To minimize the number of required mapping stimulations before the target optimization step. \n\n5. It appears that the method uses random locations to first map the ORFs. Would it be possible to optimize the selection of locations for the ORF mapping based on the desired targets for later stimulation?\n\n6. If Algorithm 1 needs to be run many times with different seeds, how long does it take to determine the optimal stimulation pattern? How practical is this in an experimental setup?\n\n7. Clarification on Figure 3: Line 199 says increasing density of neurons, and Fig. 3a, shows 50 - 150 neurons (in the same 250 x 250 um plane) – would be nice to explicitly state that is what was done in simulation, and state what the density therefore is. Why were these density ranges chosen (for different example brain regions)? Are the number of neuron-neuron connections a potential factor here? \n\n8. Is the performance increase the same for lower laser power in the nuclear condition? Ie if we used nuclear stimulation at lower power would we do just as well as optimized? Figure 4 & 5 indicate the optimization chose overall lower laser powers than initial settings.  \n", "limitations": "The authors discussed some trade-offs between the time spent mapping the ORFs and the need for single-cell precision stimulations. They mention potential future work to minimize the time needed for Bayesian target optimization, but fail to detail exact time requirements for the current method. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Originality: This work is well-motivated by the growing number of methods used for direct neuronal stimulation using photostimulation to determine e.g. causal effects of single neuron activity or the functional structure of a neural circuit. Here, rather than address the inference of neural function or connectivity itself, the authors present a solution to the problem of off-target stimulation (OTS), an experimental reality where neural targets not selected for stimulation are stimulated by the laser anyway (due to e.g. poor spatial resolution). The idea of using Bayesian optimization (BO) and finding locations for off-soma stimulations that reduce OTS is very appealing and the method proposed here with BO, GPs to model neural responses, and optimization of stimuli location shows a nice original combination of techniques to address this problem.\n\nQuality: The methods used (GPs, BO, gradient descent for optimization) are all well-established and robust. The inference method for the ORFs with the additional non-negativity constraint is nicely explained. The results clearly demonstrate an advantage of this method over traditional nuclear stimulation for the single-neural-target case in particular. \n\nClarity: Overall the paper is very well written. Ideas from the literature are well-sourced and the authors make it clear how their approach uses some similar techniques and where their method is novel. The figures are good quality as well. The code is relatively clear, though it could use additional documentation or commenting for ease of use. \n\nSignificance: As more neuroscience researchers turn to photostimulation for causal testing of hypotheses, this method should prove useful for the cases where laser accuracy is insufficient for selective stimulations. ", "weaknesses": "One weakness is the lack of experimental data, though it is understandably still difficult to obtain. While the authors mentioned slice data, this was only used as input to their simulations. Looking at the various ORFs mentioned in the supplement, is it reasonable to think that only 4 types of ORFs in the simulation is enough? And if these ORFs that were experimentally obtained were done so via stimulating a grid nearby, would that change the estimation of 1 neuron's ORF if another neuron nearby was also being stimulated (in an ensemble)? In vivo, one would expect that confounding neural activity could occur due to connectivity in addition to laser power spillover. \n\nIt is unclear if, at the end, this method is useful for ensemble stimulation in addition to single target paradigms. While the initial motivation and explanation of the method appears to be geared towards both cases (including Algorithm 1), it is not clear if in Figures 4-5 the optimization of each 'single-target hologram' is the end result of an optimization across multiple targets, or if the optimization was made per target. The different factors examined here ('density', size of ensemble, xy plus z dimensions) are all nicely presented but showing the method works in the 'worst case' (high density, larger ensembles, fully 3D, ...) would make this work stronger. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks very much for your insightful comments and feedback! Unfortunately we had to cut much of our response to meet the character limit -- we would have liked to address every point and with more detail.\n\n**Is it reasonable to think that 4 ORFs are enough?**\n\nWe do not have a sense of the variability of optogenetic receptive field (ORF) shapes beyond the slice recordings that have been made. However, we have attempted to account for this by sampling over many different GPs in our simulations, which should extensively explore the space of ORFs.\n\n**It is not clear if in Figures 4-5 the optimization of each single-target hologram is the end result of an optimization across multiple targets, or if the optimization was made per target**\n\nThank you for highlighting this ambiguity. The optimisation was made for each target individually, and we will further clarify this in the figure caption and main text.\n\n**How does this scale? Can this be done for thousands of neurons? Ensembles of 50?**\n\nThank you for raising this important point. We first note that experiments of such scale (with thousands of neurons being stimulated) are beyond what we expect the current implementation of Bayesian target optimisation can handle. Such scales could become feasible computationally if more advanced GP estimation techniques were employed, though it's not clear to us that biological experiments of this size (where a single SLM can flexibly deliver 2p excitation to thousands of neurons individually, within a single stimulation field, and with high spatiotemporal fidelity) are feasible with existing technology. \n\nWhile our objective in this submission was to establish the feasibility of overcoming off-target stimulation computationally rather than achieving the fastest possible implementation, we have provided some preliminary runtime data in the attached supplement (see **Table 1** in the rebuttal pdf), and we plan to extend the characterisation for the camera-ready copy. \n\nAlso note that both ORF inference and target optimisation are performed sequentially, and thus major speed gains could be made by parallelising ORF inference over different neurons, and target optimisation over different desired ensembles. \n\nThe runtimes given in **Table 1** demonstrate suitability for typical all-optical experiments using GCaMP6 indicators [6], the most widely used class of calcium sensors. We believe there is room for further optimisation (e.g. with more parallelisation), and in future work we plan to push the efficiency of the technique for use with the class of faster GCaMP8 indicators. We hope this also answers Question 6 from the reviewer, which we do not have space to respond to in a separate comment.\n\n**How does the exact neuron expression matter to the ORF shape? Can a simplified version be used during the 2nd step of target optimization?**\n\nWhile at a fixed power and in two dimensions some ORFs might be reasonably well-described by a simpler parametric model (e.g. a circle), in three dimensions the shape of the ORF can change so drastically and with unpredictable asymmetries that a single parametric model is unlikely to capture their varying shape with much accuracy at all (see supplemental figures S3-S6). Thus, we believe we have taken the simplest nonparametric approach (an RBF Gaussian process) that adequately describes the data.\n\nDuring the second step of target optimisation, a simple parametric ORF model might go some way in better repositioning holograms compared to naive nuclear stimulation, but for achieving the highest performance possible our experiments have found that it is certainly advantageous to exploit asymmetries in the ORF shapes, especially in three dimensions. For example, with the neuron shown in supplemental figure S4, it is clearly better to learn during an experiment that one should stimulate this neuron deeper in the tissue (if trying to avoid other neurons) than shallower.\n\n**Are there other references used for establishing priors beyond the few shown here?**\n\nThe convention in the two-photon optogenetics field has been to characterise the ability to stimulate a neuron by moving only in one dimension at a time (radially or axially, i.e. \"left and right\" or \"up and down\"; see e.g. [7-9]). The only existing reference of ORFs profiled at this level of detail (i.e. a complete grid) is [10, Figure 1G], which is similar data to what we present in the supplement.\n\n**Would it be possible to integrate spatial information (from e.g. calcium imaging) into the ORF prior?**\n\nWe thank the reviewer for making this interesting suggestion. In most cases, calcium sensors and opsins are expressed through separate viral vectors, and neurons are not guaranteed to express both proteins at once. However, some research groups have developed constructs that fuse the two proteins together directly [10], in which case the relative brightness of a calcium transient across the cell soma and proximal dendrites might provide some hint at different locations that the cell could be better activated, though it remains to be seen how useful this would actually be in practice.\n\n**Would it be possible to optimize the selection of locations for the ORF mapping based on the desired targets for later stimulation?**\n\nYes! We definitely think this can and should be done in practice. If one only wishes to stimulate a desired subset of neurons, then ORF mapping should only be performed in neighbourhoods surrounding those neurons.\n\n**If we used nuclear stimulation at lower power would we do just as well as optimized?**\n\nWe performed a control analysis by matching the nuclear stimulation laser power to the average optimised laser power (see **Figure 1** in the rebuttal supplement) and reran the mapping and target optimisation phases multiple times while randomly reassigning neurons different ORFs. This process showed that we maintain a very similar improvement in performance as what was shown in the main text.\n\nThanks again!"}, {"review_id": "C7jxWg0Tcm", "reviewer": "Reviewer_5v9K", "summary": "The problem the authors tackle in their manuscript is the problem of target selection in holographic optogenetics. Briefly, of the many neurons in a field of view, many experiments require the selective stimulation of a small subset of cells, while minimizing off-target stimulation that may muddy the interpretation of the data. This question is of considerable importance to the field of systems neuroscience, and several experimental methods have been pursued trying to solve it, such as the use of soma-localized opsins or sparse expression of optogenetic protein. The authors present in this manuscript a complementary, Bayesian computational approach that takes into account differences in sensitivity of stimulation of different neurons in a given FOV to optimize laser intensity of position to minimize off-target stimulation while maintaining sufficient target stimulation. The authors use gaussian processes to GPs to infer the response properties of neurons based on stimulation parameters. The authors validate this approach in simulation and in real optogenetic experiments.", "questions": "Are the real-world validations provided using soma-localized channelrhodopsin, or is there channelrhodopsin in the processes? Does the method provide a greater benefit in one situation versus the other?\n\nAs stated in weaknesses section:\nHow sensitive is an OTS-minimized stimulation protocol to sample motion? How does this compare to soma-targeting? How about to changes in neuron sensitivity to stimulation? Can these be incorporated into optimization?", "limitations": "The authors have adequately addressed the limitations.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "Originality: The authors introduce a novel computation algorithm to address a pressing issue in systems neuroscience; that of off-target stimulation. While the authors note that similar methods have been used to infer receptive fields for other types of stimulation, the authors cleverly employ these methods in a new way for inferring the responses of neurons to optogenetic stimulation to address a technical challenge facing the field. This software approach complements experimental approaches and can be readily combined with them. Similar computational approaches of thus far remained lacking, so the work is pioneering in this regard.\n\nClarity: The manuscript is clear with respect to how it poses the question being addressed, how it delineates its approach, and its statement of results.\n\nSignificance: The authors make a significant contribution to a timely question for systems neuroscience and will likely see significance use by labs performing holographic optogenetic stimulation, especially for applications in which minimizing off-target stimulation is crucial. The method could also see use in designing stimulation protocols to precisely pattern the activity of an ensemble of neurons, although that is a future direction.", "weaknesses": "Optimizing stimulation according to one parameter (in this case, minimization of off-target stimulation) often necessitates trade-off to other parameters. The authors bring up one such parameter, the time required to map the ORF and perform the optimization. While this consideration will vary from experimenter to experimenter, as the authors suggest, other considerations might not. For example, by moving stimulation locations off-center relative to the soma, stimulation reliability might become more sensitive to sample motion. Furthermore, the authors assume a ORF fixed during time; the manuscript would benefit from additional consideration for how the ORF changes with sustained stimulation, as desensitization of the channelrhodopsin could significant affect the ORF. In general, the authors should further explore the limitations and trade-offs of this method. This weakness is minor, however, since I assume that these other considerations could be incorporated as parameters to be optimized.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you again for spending the time to review our submission and for providing your comments!\n\n**Sensitivity to sample motion**\n\nThanks for raising this relevant point. It is true that sample motion could affect the optimality of the computationally identified stimuli. In practice, we would therefore recommend performing online motion correction [11] immediately before stimulation, so that the learned optogenetic receptive field (ORF) posterior is as closely aligned in space to the true \"real-time\" receptive field as possible.\n\n**Sensitivity to opsin nonstationarity**\n\nDesensitisation of opsins following prolonged illumination has indeed been observed [12]. However, after the ORFs have been mapped and holographic stimuli have been optimised, it is straightforward to periodically recompute the posterior as the experiment proceeds. This would provide an updated estimate of the laser power required to maintain optimal precision. Note that this does not require redoing the whole ORF mapping phase as we would not expect the shape of the ORF to change drastically beyond desensitisation.\n\nWe will update the manuscript to discuss these important points.\n\n**Are the real-world validations provided using soma-localized channelrhodopsin?**\n\nYes, the validation using experimental data is performed with ChroME2f, a state-of-the-art soma-targeted opsin [13].\n\n**Does the method provide a greater benefit in one situation versus the other?**\n\nWe have not explored applying our technique in the case of an opsin that is not localised to the soma, because in such experiments the precision of any optogenetic manipulation is extremely poor [7, 8]. We therefore recommend that any experimenter wanting to achieve single-cell precision first change to a soma-localised opsin before attempting any computational optimisation.\n\nThanks again for your comments.\n"}, {"review_id": "Ue5dwg69aI", "reviewer": "Reviewer_G2K6", "summary": "This paper proposes a new method for limiting off-target optogenetic stimulation based on Gaussian Process modeling. In holographic photostimulation, the goal is to excite specific neurons via targeted laser light, but widespread expression of opsins may result in additional neurons not in the desired population (\"off-target\" neurons) also firing action potentials. The proposed method uses approximate GP inference (MAP estimation) in combination with a novel gradient-based optimization method to refine target locations in order to minimize the $L_2$ distance between the desired and evoked patterns of stimulation.\n\nThis is a very nice paper that combines good modeling with a very thoughtful approach to experimental realities. While direct experiments will be needed to fully test its efficacy, it stands to make an important contribution to this particular neuroscience application.", "questions": "- Line 102 defines $\\sigma$ as the sigmoid function. Does that mean logistic? (E.g, $\\tanh$ is also sigmoid in shape but clearly not intended.)\n- It is not stated in the paragraph surrounding Eq. 2 that $I$ is the laser power. It might help to simply state $\\mathbf{x} = (c_1, c_2, I)$.\n- Line 105 states that one should have $g_n(x) \\ge 0$ so that stimulation should always have a positive effect, but in the supplement, this condition is waived in the single-target case. I assume this is because one can compensate in this case by taking $\\theta = \\inf_{x} g_n(x)$? If so, it might help to articulate this in the supplement.\n- what is the MAP version of the ORF $\\hat{g}_n$? This is just $\\hat{g}_n(\\mathbf{X})$, with $\\mathbf{X}$ the matrix of test data points?\n- Why use MSE and not cross-entropy, which would be the standard assumption for a desired set of observations under a Bernoulli model?", "limitations": "- Much of the inference algorithm is quite closely tied to the experimental setup. This is both a strength and a limitation, since the proposed work is somewhat unlikely to find broader applications.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- Careful modeling of features of real experiments.\n- Use of a scalable GP approach. \n- Clear exposition.", "weaknesses": "- Much of the inference algorithms seems specific to some of the particular modeling choices made (particularly the need for a convex loss in Line 1 of Algorithm 1).\n- Figs 4-5 should clearly be labeled as _simulation experiments_ based on real data; this is clearly explained in ll. 212-214 but elided in, e.g., ll. 254-56, which makes it seem as if the optimization was performed and validated as part of data collection.\n- I was suprised not to see references to Fletcher and Rangan (2014) and Draelos and Pearson (2020) in the related work section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks again for your helpful feedback.\n\n**Figs 4-5 should clearly be labeled as _simulation experiments_ based on real data**\n\nThanks for noting this ambiguity. We will update the figure captions and lines 254-256 accordingly.\n\n**Missing references**\n\nWe thank you for pointing out that we had not cited these papers, these were poor and unintended omissions on our part that we will fix. Notably, Draelos and Pearson (2020) is a recent inspiration for our work that we had absolutely meant to cite in the related work section.\n\n**Clarification of which sigmoid function is being used**\n\nYes, we did indeed mean the logistic sigmoid and we will update the main text to clarify this.\n\n**Clarification of laser power parameter**\n\nThanks for pointing out that this is unclear -- we will follow your suggestion!\n\n**Non-negativity is waived in the single-target case**\n\nWhen we waive the non-negativity constraint in the single-target case, we actually drop the threshold parameter entirely, much like in Gaussian process models of visual receptive fields [14]. This is because now if a point on the optogenetic receptive field (ORF) becomes inhibitory (by taking a negative value), it will not conflict with excitation from any other hologram, and therefore will still effectively model the response to optogenetic stimulation. We will clarify this in the supplement.\n\n**What is the MAP version of the ORF $\\hat g_n$?**\n\nThe MAP version of the ORF is $\\hat g_n(\\mathbf{X})$ (i.e. the ORF evaluated at $\\mathbf{X}$), but where $\\mathbf{X}$ is the set of \"training\" points probed during the ORF mapping phase (perhaps this is what you meant already). We will be explicit about this in section 3.2.\n\n**Why use MSE and not cross-entropy?**\n\nNo particular reason, other than that it appeared to be an effective loss function in practice!\n\nThanks again for all your comments."}]

Author rebuttal:
Thank you all for your excellent feedback and positive evaluation of our submission! We presented Bayesian target optimisation, a computational approach to overcoming off-target stimulation in two-photon optogenetics experiments. We are delighted that every reviewer clearly understood the motivation, methodological contributions, and experimental relevance of our work.

Below we consider the comments that are specific to each review, but first we would like to address the fact that multiple reviewers had concerns regarding experimental data. In most cases, the ideal data for calibrating the internal model components unfortunately do not yet exist or are not currently available to us. We are actively working with experimental collaborators to validate our technique _in vivo_. However, obtaining new experimental data to address the specific concerns of individual reviewers is (hopefully understandably) beyond the scope of this NeurIPS submission. We have therefore done our best to provide, where possible, additional references to what is known in the literature. Nevertheless, we believe that our chosen validation approach (combining data from separate real experiments) comes as close to a direct experimental validation as is currently feasible, and is more closely tied to real experimental data than any existing computational work at NeurIPS in this area [1-3].

Please note that we have attached some additional data in response to Reviewer 4tSm, and have collected all references in this general rebuttal in order to meet character limits. Thanks again all for your work in reviewing our paper.

**References**

1. Hu, T., & Chklovskii, D. (2009). Reconstruction of sparse circuits using multi-neuronal excitation (RESCUME). Advances in Neural Information Processing Systems, 22
2. Shababo, B., Paige, B., Pakman, A., & Paninski, L. (2013). Bayesian inference and online experimental design for mapping neural microcircuits. Advances in Neural Information Processing Systems, 26.
3. Draelos, A., & Pearson, J. (2020). Online neural connectivity estimation with noisy group testing. Advances in Neural Information Processing Systems, 33, 7437-7448.
4. Rickgauer, J. P., & Tank, D. W. (2009). Two-photon excitation of channelrhodopsin-2 at saturation. Proceedings of the National Academy of Sciences, 106(35), 15025-15030.
5. Losonczy, A., & Magee, J. C. (2006). Integrative properties of radial oblique dendrites in hippocampal CA1 pyramidal neurons. Neuron, 50(2), 291-307.
6. Russell, L. E., Dalgleish, H. W., Nutbrown, R., Gauld, O. M., Herrmann, D., Fişek, M., ... & Häusser, M. (2022). All-optical interrogation of neural circuits in behaving mice. Nature Protocols, 17(7), 1579-1620.
7. Baker, C. A., Elyada, Y. M., Parra, A., & Bolton, M. M. (2016). Cellular resolution circuit mapping with temporal-focused excitation of soma-targeted channelrhodopsin. Elife, 5, e14193.
8. Shemesh, O. A., Tanese, D., Zampini, V., Linghu, C., Piatkevich, K., Ronzitti, E., ... & Emiliani, V. (2017). Temporally precise single-cell-resolution optogenetics. Nature neuroscience, 20(12), 1796-1806.
9. Mardinly, A. R., Oldenburg, I. A., Pégard, N. C., Sridharan, S., Lyall, E. H., Chesnov, K., ... & Adesnik, H. (2018). Precise multimodal optical control of neural ensemble activity. Nature neuroscience, 21(6), 881-893.
10. Bounds, H. A., Sadahiro, M., Hendricks, W. D., Gajowa, M., Gopakumar, K., Quintana, D., ... & Adesnik, H. (2023). All-optical recreation of naturalistic neural activity with a multifunctional transgenic reporter mouse. Cell Reports, 42(8).
11. Pnevmatikakis, E. A., & Giovannucci, A. (2017). NoRMCorre: An online algorithm for piecewise rigid motion correction of calcium imaging data. Journal of neuroscience methods, 291, 83-94.
12. Marshel, J. H., Kim, Y. S., Machado, T. A., Quirin, S., Benson, B., Kadmon, J., ... & Deisseroth, K. (2019). Cortical layer–specific critical dynamics triggering perception. Science, 365(6453), eaaw5202.
13. Sridharan, S., Gajowa, M. A., Ogando, M. B., Jagadisan, U. K., Abdeladim, L., Sadahiro, M., ... & Adesnik, H. (2022). High-performance microbial opsins for spatially and temporally precise perturbations of large neuronal networks. Neuron, 110(7), 1139-1155.
14. Park, M., Horwitz, G., & Pillow, J. (2011). Active learning of neural response functions with Gaussian processes. Advances in neural information processing systems, 24.

Meta-review: This work aims to identify the use of Bayesian optimization to develop a targeting approach that minimizes the stimulation of non-target neurons that can occur during in-vivo neural stimulation using optical methods. The methods itself is sound, and the reviewers noted that the exploration is timely and complements well the development of optical methods in systems neuroscience.

The main weakness of this work is the difficulty in knowing just how well the method, with the approximations and modeling choices of the light-neuron interactions, will work in real datasets. The example data uses primarily simulated data (either toy data or sampling from in-vivo data). Granted that real data are difficult to obtain, there is more realistic simulation software available (Song et al. 2022 J. Neuo Methods), and simple characterizations like showing that the offsets optimized by the system are useful even under in vivo motion conditions, which can throw off targeting quite a bit. Regardless the reviewers considered the benefits of the work to outweight the shortcomings and so I recommend this work be accepted.
Decision: Accept (spotlight)