Inferring connectivity in neuronal networks remains a key challenge in statistical neuroscience. The “common input” problem presents a major roadblock: it is difficult to reliably distinguish causal connections between pairs of observed neurons versus correlations induced by common input from unobserved neurons. Available techniques allow us to simultaneously record, with sufficient temporal resolution, only a small fraction of the network. Consequently, naive connectivity estimators that neglect these common input effects are highly biased. This work proposes a “shotgun” experimental design, in which we observe multiple sub-networks briefly, in a serial manner. Thus, while the full network cannot be observed simultaneously at any given time, we may be able to observe much larger subsets of the network over the course of the entire experiment, thus ameliorating the common input problem. Using a generalized linear model for a spiking recurrent neural network, we develop a scalable approximate expected loglikelihood-based Bayesian method to perform network inference given this type of data, in which only a small fraction of the network is observed in each time bin. We demonstrate in simulation that the shotgun experimental design can eliminate the biases induced by common input effects. Networks with thousands of neurons, in which only a small fraction of the neurons is observed in each time bin, can be quickly and accurately estimated, achieving orders of magnitude speed up over previous approaches.
It is now possible to image hundreds of neurons simultaneously at high spatiotemporal resolution [1] or tens of thousands of neurons at low spatiotemporal resolution [2]. The number of recorded neurons is expected to continue to grow exponentially [3]. This, in principle, provides the opportunity to infer the “functional” (or “effective”) connectivity of neuronal networks, i. e. a statistical estimate of how neurons are affected by each other, and by a stimulus. The ability to accurately estimate large, possibly time-varying, neural connectivity diagrams would open up an exciting new range of fundamental research questions in systems and computational neuroscience [4]. Therefore, the task of estimating connectivity from neural activity can be considered one of the central problems in statistical neuroscience. Naturally, such a central problem has attracted much attention in recent years (see section 8). Perhaps the biggest challenge here involves the proper accounting for the activity of unobserved neurons. Despite rapid progress in simultaneously recording activity in massive populations of neurons, it is still beyond the reach of current technology to simultaneously monitor a complete large network of spiking neurons at high temporal resolution. Since connectivity estimation relies on the analysis of the the activity of neurons in relation to their inputs, the inability to monitor all of these inputs can result in persistent errors in the connectivity estimation due to model miss-specification. More specifically, “common input” errors, in which correlations due to shared inputs from unobserved neurons are mistaken for direct, causal connections, plague most naive approaches to connectivity estimation. Developing a robust approach for incorporating the latent effects of such unobserved neurons remains an area of active research in connectivity analysis (see section 8). In this paper we propose an experimental design which can greatly ameliorate these common-input problems. The idea is simple: if we cannot observe all neurons in a network simultaneously, perhaps we can instead observe many overlapping sub-networks in a serial manner over the course of a long experiment. Then we can use statistical techniques to patch the full estimated network back together, analogous to “shotgun” genetic sequencing [5]. Obviously, it is not feasible to purposefully sample from many distinct sub-networks at many different overlapping locations using multi-electrode recording arrays, since multiple re-insertions of the array would lead to tissue damage. However, fluorescence-based imaging of neuronal calcium [6,7] (or, perhaps in the not-too-distant future, voltage [8]) makes this approach experimentally feasible. For example, such a shotgun approach could be highly beneficial and relatively straightforward to implement using a 3D acousto-optical deflector microscope [1]. Using such a microscope, one can scan a volume of 400 × 400 × 500 μm, which contains approximately 8000 cells. In normal use, the microscope’s 50kHz sampling rate allows for a frame rate of about 6Hz when scanning the entire volume. Unfortunately, this frame rate is too low for obtaining reliable connectivity estimates, which requires a frame rate of at least 30Hz [9]. However, we can increase the effective frame rate to 30Hz by using a shotgun approach. We simply divide the experimental duration into segments, where in each segment we scan only 20% of the network. As a side benefit of this shotgun approach, photobleaching and phototoxicity (two of the most important limitations on the duration of these experiments [10]) are reduced, since only a subset of the network is illuminated and imaged at any given time. Connectivity estimation with missing observations is particularly challenging (section 9). Fortunately, as we show here, given the shotgun sampling scheme, we do not have to infer the unobserved spikes. We considerably simplify the network model loglikelihood using the expected loglikelihood approximation [11–13], and a generalized Central Limit Theorem (CLT) [14] argument to approximate the neuronal input as a Gaussian variable when the size of the network is large. This approximate loglikelihood and its gradients depend only on the empiric second order statistics of the spiking process (mean spike rate and spike correlations). Importantly, these approximate sufficient statistics can be calculated, even with partial observations, by simply “ignoring” any unobserved activity (section 3. 6). In order to obtain an accurate estimation of the connectivity, posterior distributions involving this simplified loglikelihood (along with various types of prior information about network connectivity) can be efficiently maximized. Using a sparsity inducing prior on the weights, we demonstrate numerically the effectiveness of our approach on simulated recurrent networks of spiking neurons. First, we demonstrate that the shotgun experimental design can largely eliminate the biases induced by common input effects (section 4). Then, we show that we can quickly infer connectivity for large networks, with a low fraction of neurons observed in each time bin (section 5). For example, our algorithm can be used to infer the connectivity of a sparse network with O (103) neurons and O (105) connections, given O (106) time bins of spike data in which only 10% − 20% of the neurons are observed in each time bin. On a standard laptop, simulating such a network takes about half an hour, while inference takes a few minutes. This is faster than previous approaches by orders of magnitude, even when all spikes are observed (section 6. 2). Our parameter scans suggest that our method is robust, and could be used for arbitrarily low observation ratios and an arbitrarily large number of neurons, given long enough experiments. We will discuss the outlook for experimental realizations of the proposed approach below, after presenting the basic methodology and simulated results. The supplementary material S1 Text contains the full details of the mathematical derivations and the numerical simulations.
We use a Bayesian approach to infer the unknown weights. Suppose initially, for simplicity, that all spikes are observed and that there is no external input (G = 0). In this case, the log-posterior of the weights, given the spiking activity, is ln P (W | S, b) = ln P (S | W, b) + ln P 0 (W) + C, (4) where ln P (S∣W, b) is the loglikelihood, P0 (W) is some prior on the weights (we do not assume a prior on the biases b), and C is some unimportant constant which does not depend on W or b. Our aim is to find the Maximum A Posteriori (MAP) estimator for W, together with the Maximum Likelihood (ML) estimator for b, by solving max W, b ln P (W | S, b). (5) If S is fully observed, this problem can be straightforwardly optimized without requiring an approximation (though the optimization procedure can be slow). However, our goal is to provide an estimate when only a subset of S is observed. This cannot be easily done using standard method. To see this, we examine the likelihood of a GLM (recalling Eqs (1) and (2) ), ln P (S | W, b) = ∑ i = 1 N ∑ t = 1 T ln [ e S i, t U i, t 1 + e U i, t ] (6) = ∑ i = 1 N ∑ t = 1 T [ S i, t U i, t - ln (1 + e U i, t) ]. (7) This likelihood (Eq 7), and its gradients, both contain a sum over weighted spikes in Ui, t (the WS⋅, t−1 term in Eq 2), that cannot be evaluated if some spikes are missing, unless the missing spikes are accurately inferred (section E in S1 Text). However, methods for inferring these missing spikes are typically slow, and do not scale well. To circumvent these issues, we will show the loglikelihood can be approximated with a simple form, under a few reasonable assumptions. Importantly, this simple form can be easily calculated even if there are missing observations (the full derivation is in section 2. 1). Using an extension of the techniques in [11–13], we develop an approximation to the likelihood based on the law of large numbers (the “expected loglikelihood” approximation) together with a generalized Central Limit Theorem (CLT) argument [14], in which we approximate the neuronal input to be Gaussian near the limit N → ∞; then we calculate the “profile likelihood” maxbln P (S∣W, b), in which the bias term has been substituted for its maximizing value. The end result is max b ln P (S | W, b) ≈ T ∑ i = 1 N [ ∑ j = 1 N [ W i, j Σ i, j (1) ] - h (m i) 1 + π 8 ∑ k, j W i, j Σ k, j (0) W i, k ], (8) where we defined the mean spike probability, spike covariance, and the entropy function, respectively: m i≜⟨ S i, t ⟩ T (9) Σ i, j (k) ≜⟨ S i, t S j, t - k ⟩ T - m i m j (10) h (m i) ≜- m i ln m i - (1 - m i) ln (1 - m i). (11) A few comments: Importantly, the profile loglikelihood (Eq 8) depends only on the first and second order moments of the spikes m and Σ (k) for k ∈ {0,1}. When all of the neurons in the network are observed, these moments can be computed directly, and therefore the empirical moments are approximate sufficient statistics, whose value contains all the information needed to compute any estimate of W. As we explain in section 3, these empirical moments can be estimated even if only a subset of the spikes is observed. As we show in section A in S1 Text, the profile loglikelihood (Eq 8) is concave, so it is easy to maximize the log-posterior and obtain the MAP estimate of W. This can be done orders of magnitude faster than in the standard MAP estimate (section 6. 2), since Eq 8 does not contain a sum over time, as the original loglikelihood (Eq 7). Moreover, the optimization problem of finding the MAP estimate can be parallelized over the rows of W. max b ln P (W | S, b) = ∑ i max b ln P (W i, · | S, b), (12) because the profile loglikelihood (Eq 8) decomposes over the rows of W, as does the L1 prior we will use here (Eq 46). As we show in section A in S1 Text, we can straightforwardly differentiate Eq 8 to analytically obtain the gradient, Hessian, and even the maximizer of this profile loglikelihood, which is the maximum likelihood estimate of W. However, due to the nature of the integral approximation we make in Eq 14, more accurate results are obtained if we first differentiate the original loglikelihood (Eq 7), and then use the expectation approximation (together with the generalized CLT argument). This results in an adjustment of the loglikelihood gradient (section D in S1 Text). A novel aspect of this work is that we apply the Expected LogLikelihood (ELL) approximation to a GLM with a bounded logistic rate function (Eq 1), which allows us to infer connectivity in recurrent neural networks. In contrast, previous works that used the ELL approximation [11–13] focused on single neuron responses, with an emphasis on either a Poisson neuron model with an exponential rate function, or simpler linear Gaussian models. Such models are less suitable for recurrent neural networks. Exponential rate functions cause instability, as the activity tends to to diverge, unless both the weights and the time bins are small. Linear networks are not a very realistic model for a neural network, and do not perform well in inferring synaptic connectivity [19]. Though we assumed a logistic neuron model (Eq 1), similar results can be derived for any spiking neuron model for which 1−f (x) = f (−x). This is explained in section A. 3 S1 Text. Though we assumed the network does not have a stimulus (G = 0), one can be incorporated into the inference procedure. To do so, we treat the stimulus X⋅, t simply as the activity of additional, fully observed, neurons (albeit Xi, t ∈ ℝ while Si, t ∈ {0,1}). Specifically, we define a new “spikes” matrix Snew ≜ (S⊤, X⊤) ⊤, a new connectivity matrix W new≜ (W G 0 D × N 0 D × D), and a new observation matrix Onew ≜ (O⊤, 1T×D) ⊤. Repeating the derivations for Snew, Wnew and Onew, we obtain the same profile loglikelihood. Once it is used to infer Wnew, we extract the estimates of W and G from their corresponding blocks in Wnew. For simplicity and efficiency, we chose to focus on MAP estimates. However, other types of estimators and Bayesian approaches (e. g. , MCMC, variational Bayes) might be used with this approximate loglikelihood, and should be explored in future work. As we showed in section 2, in order to infer network connectivity, we just need to estimate the first and second empiric spike statistics, defined in Eqs (9) – (10). These statistics cannot be calculated exactly if some observations are missing; in this case they must be estimated, as we discuss in section 3. 6 below. First, though, it is useful to discuss a few concrete examples of the partial network observation schemes we are considering (Fig 1). We discuss the pros and cons of each scheme in terms of both inferential and experimental constraints.
In this section we use a toy network with N = 50 neurons to visualize the common input problem, and its suggested solution—the “shotgun” approach. Errors caused by common inputs are particularly troublesome for connectivity estimation, since they can persist even as T → ∞. Therefore, for simplicity, we work in a regime where the experiment is long and data is abundant (T = 5⋅108 timebins). In this regime, any prior information we have on the connectivity becomes unimportant so we simply use the Maximum Likelihood (ML) estimator. We chose the weight matrix W to illustrate a “worst-case” common input condition (Fig 2A). Note that the upper-left third of W is diagonal (Fig 2B): i. e. , neurons i = 1, …, 16 share no connections to each other, other than the self-connection terms Wi, i. However, we have seeded this W with many common-input motifs, in which neurons i and j (with i, j ≤ 16) both receive common input from neurons k with k ≥ 17. If we use a “shotgun” approach and observe the whole network with pobs = 16/50 with a fully random observation scheme, we obtain a good ML estimate of the network connectivity, including the 16 × 16 upper-left submatrix (Fig 2C). Now, suppose instead we concentrate all our observations on these 16 neurons, so that pobs = 1 within that sub-network, but the other neurons are unobserved. If common input was not a problem, our estimation quality should improve on that submatrix (since we have more measurements per neuron). However, if common noise is problematic, then we will “hallucinate” many nonexistent connections (i. e. , off-diagonal terms) in this submatrix. Fig 2D illustrates this phenomenon. In contrast to the shotgun case, the resulting estimates are significantly corrupted by the common input effects. Next, we quantitatively test the performance of the Maximum A Posteriori (MAP) estimate of the network connectivity matrix W using a detailed network model with biologically plausible parameters from the mouse visual cortex. Details on the network parameters, simulation details and definitions of the quality measures are given in section B in S1 Text. We use the inference method described in section 2, with a sparsity inducing prior (section C in S1 Text) on a simulated network with GLM neurons (Eqs (1) – (2) ). First, in Fig 3, we examine a small GLM network with N = 50 observed neurons, with an experiment length of 5. 5 hours. As can be seen, the weight matrix can be very accurately estimated for high values of observation probability pobs, and reasonably well even for low value of pobs. For example, even if pobs = 0. 04, and only two neurons are observed in each timestep, we get a correlation of C ≈ 0. 84 between inferred weights and the true weights, and the signs of the non-zero weights are only wrong only for 4 weights (out of 448 non-zero weights). When pobs is decreased, the variance of the estimation increases, more weak weights are inferred as zero weights (and vice versa), and we also see more “shrinkage” of the non-diagonal weights (a decreased magnitude of the non-zero weights) due to the L1 penalty imposed on them (Eq 47 in S1 Text). In Fig 4 we demonstrate that our method works well even if the neuron model is not a GLM, as we assume, but a Leaky Integrate and Fire (LIF) neuron model (Fig 4). The model mismatch results in a weight mismatch by a global multiplicative constant, and in a worse estimate of the diagonal weights, due to the hard reset in the LIF model. Besides these issues, inference results are both qualitatively and quantitatively similar to results in the GLM network In Fig 5, we examine another GLM network with N = 1000 observed neurons, which is closer to the scale of the number of recorded neurons in current calcium imaging experiments (see activity simulation in Fig S1 S1 Text). The experiment duration is again 5. 5 hours. Results are qualitatively the same as the case of N = 50, except performance is somewhat decreased (as we have more parameters to estimate). Additional information is available in Fig 6. On the left (A, D, G), we see that the algorithm converges properly to a single solution. In the middle panels (B, E, H), we see that for pobs = 1 we have very good performance (in terms of area under the ROC), but this performance declines for the excitatory weights as pobs decreases. The inhibitory weights are correctly detected much better than the excitatory weights. This is because most excitatory weights are much weaker, as can be seen on the right column (C, F, I). In that column, we observe that strong weights are more easily detected than weak weights. Specifically, around the median of the excitatory weight distribution (0. 178), we detected 99. 9%, 34. 9% and 16. 1% of all the weights, when pobs = 1,0. 2 and 0. 1, respectively. Next, in Fig 7 we quantify how inference performance changes with parameters. We vary the number of neurons, N, observation probability pobs, mean firing rate m and connection sparsity pconn. For the given parameters N, pobs and m, performance monotonically improves when T increases. These scans suggest we can maintain a good quality of connectivity estimation for arbitrarily large or small values of N or pobs, respectively—as long as we sufficiently increase T. Note there is a lower bound on T, below which estimation does not work. Looking at Fig 7, we find that approximately, this lower bound scales as T ∝ N p obs 2. (18) Above this lower bound, estimation quality gradually improves with T. Moreover, in order to maintain good estimation quality (up to some saturation level) above this bound, T should be scaled as T ∝ N p obs 2 m 2. (19) This scaling can be explained intuitively. Our main sufficient statistic is the partially observed spike covariance Σ ˜ (k) (Eq 17). Each component (i, j) of Σ ˜ (k) contains a sum of all the observed spike pairs (T⟨Oi, t Oj, t−k Si, t Sj, t−k⟩T) divided by the number of observed neurons (T⟨Oi, t Oj, t−k⟩T). The total number of observed neuron pairs is approximately N T p obs 2 (ignoring observation correlations), and the total number of observed spike pairs is approximately N T p obs 2 m 2 (ignoring spike correlations, and assuming the firing rate is not very high), where T is measured in time bins. The total number of components in Σ ˜ (k) is N2. Therefore, in each component of Σ ˜ (k), the average number of observed neuron pairs is T p obs 2 / N, while the average number of observed spike pairs is approximately T p obs 2 m 2 / N (except on the diagonal of Σ (0), where we have Tpobs/N neuron pairs and Tpobs m/N spikes). We conclude that the number of both observed neuron pairs and spike pairs must be above a certain threshold so that inference will be able to work properly. Above these thresholds, performance improves further when pconn is decreased (Fig 7), as this reduces the effective number of parameters we are required to estimate. For analytic results on this issue see [28]. If our goal is to infer all the input connections of a single neuron, then performance can be significantly improved if we always observe the output of that neuron. This is demonstrated in Fig 8. In this figure, we examine a single neuron with O (104) observed inputs, O (103) of which are non-zero (implementation details in S1 Text, section B. 2). The inputs are partially observed (with pobs = 1,0. 1,0. 01), but we always observe the output neuron. Therefore, the average number of observed neuron pairs and spike pairs in Σ (1) is increased to NTpobs and NmTpobs, respectively. This can improve the scaling relations in Eqs (18) and (19) to T/pobs and T/ (pobs m2), respectively, if the off-diagonal terms of Σ (0) are not too strong (since the number of observations for these components still scales with ∝ p obs 2). Thus, in Fig 8, we see that even when pobs = 0. 01 it is still possible to estimate strong weights with some accuracy, despite the large number of connections. In a real imaging experiment, we would not have direct access to spikes, as we have assumed for simplicity so far. Next, we test the estimation quality when we only have direct access to the fluorescence traces of activity (Fig 10A). The fluorescence traces were generated using a model of GCaMP6f calcium fluorescence indicator. Implementation details are described in section B. 3 in S1 Text. As can be seen in Fig 10A, 10B, our spike inference algorithm works reasonably well, both in high and low noise regimes. We then infer network connectivity both from the inferred spikes and the true spikes. As can be seen in Fig 10C-10H, using the inferred spikes usually somewhat reduces estimation performance. This is due to the temporal inaccuracy in the spike estimation. For example, in the inhibitory neurons, the higher firing rates result in more missing spikes in the inference. This causes shrinkage in the magnitude of the inferred weights, since the cross-correlation is weakened by these missing spikes. Combining this information into the inference algorithm (as in [9]), it may be possible to correct for this; we have not pursued this question further here. However, even at low observation probabilities (pobs = 0. 1), strong weights are inferred reasonably well, and the sign of synapse is usually inferred correctly for almost all nonzero weights. Therefore, weight inference is still possible at low firing rates, using current generation fluorescence imaging methods.
Neural connectivity inference has attracted much attention in recent years. One approach to this problem is direct anatomical tracing [31,32]. However, this method is computationally challenging [33]; moreover, the magnitudes of the synaptic connections (which also vary over time [34]) currently cannot be inferred this way. Another approach, on which we focus here, aims to infer the synaptic connectivity from neural activity. This activity can be either action potentials (“spikes”) [16,35–37] or calcium fluorescence traces [9,18,38–40] which are approximately a noisy and filtered version of the spikes. Various inference procedures have been suggested for this purpose. Some works use model-free empirical scores [38,40,41]. Others assume an explicit generative model for the network activity [16,30,35–37], and then infer connectivity by estimating model parameters. So far, only few works have validated the connectivity estimate with some form of “ground truth”. Gerhard et al. [19] inferred small scale anatomical connectivity, comparing different methods. A Generalized Linear Model (GLM) approach was successful, while linear models and model-free approaches failed. Volgushev et al. [42] estimated the weights of fictitious synapses (injected current). Again, a GLM-based approach outperformed simple correlation-based approaches. Lastly, Latimer et al. [43] was able to infer the magnitude of intracellular synaptic conductances, using a modified GLM. These results indicate that a GLM-based approach should be the method of choice for estimating synaptic connectivity. The task of inferring synaptic connectivity is severely hindered by technical limitations on the number of neurons that can be simultaneously observed with sufficient quality. Typically, the scanning speed of the imaging device is limited, so we cannot cover the entire network with a high enough frame rate and signal-to-noise ratio to infer spikes from the observed fluorescence traces. Previous studies indicate that at low frame rates (below 30Hz [9]), synaptic connectivity cannot be inferred. In such low frame rate regimes, one may use spike correlations or simple dynamical systems as a coarse measure of effective connectivity (e. g. , [44]), but such measures are not claimed to predict synaptic connectivity, only provide a statistical description of the network dynamics. Therefore, common approaches to infer connectivity of a neural network focus all the observations in one experiment on a small part of the network, in which all neurons are fully observed at a high frame rate. However, unobserved input into this sub-network can generate significant error in the estimation, and this error does not vanish with longer experiments. Various works aimed to deal with this persistent error: [17] inferred connectivity in a simulated two-neuron network in which one neuron was never observed; [45] inferred connectivity in a simulated network with two observed neurons and an unobserved common input; [46] inferred unobserved common input in an experimentally recorded network of 250 neurons using a GLM network with latent variables; [47] inferred connectivity in a simulated network with 100 neurons where 20−50 were never observed, with a varying degree of success. To help deal with the “common input” problem, we propose a “shotgun” approach, in which we reconstruct network connectivity by serially observing small parts of the network—where each part is observed at a high frame rate for a limited duration. Thus, despite the limited scanning speed of the imaging device, by using this method, we can extend the number of the neurons covered by the scanning device and effectively decrease the number (and therefore the effect) of unobserved common inputs. Additionally, as only a small part of the network is illuminated together, this method can potentially reduce phototoxicity and photobleaching, and allow long, possibly chronic [48], imaging experiments. We showed here that the proposed method is capable of incorporating prior information about the sparsity of synaptic connections. More specific information could be included. An abundance of such prior information is available for both connection probabilities and synaptic weight distributions as a function of cell location and identity [52]. Cutting edge labeling and tissue preparation methods such as Brainbow [53] and CLARITY [54] are beginning to provide rich anatomical data about “potential connectivity” (e. g. , the degree of coarse spatial overlap between a given set of dendrites and axons) that can be incorporated into these priors. Exploiting such prior information can significantly improve inference quality, as demonstrated in previous network inference papers [9,17,55]. For example, by adjusting the L1 regularization parameters, we can reflect such additional priors: that the probability of having a connection between two neurons typically decreases with the distance between two neurons, and that it is affected by the neuronal type. Another way to improve connectivity estimates is to use stimulus information. For example, increasing the firing rate can improve quality (Eq 19 and Fig 7), up to a limit. If the firing rate is too high, it becomes harder to infer spikes from fluorescence. A more sophisticated spatio-temporal stimulus scheme can potentially lead to significant improvements in estimation quality [56]. The type of stimulus used can also affect performance. Sensory stimulation usually affects the measured network indirectly, potentially through many layers of neuronal processing. This may result in undesirable common input (“noise correlations”). Optogenetic stimulation does not have this problem, since it stimulates neurons directly by using light sensitive ion channels. However, this type of optical stimulation can potentially interfere with optical recording. Such cross-talk can be minimized by using persistent ion channels [57] (which require only a brief optical stimulus to be activated), or more sophisticated types of stimulation schemes [58,59]. Such optogenetic approaches, coupled with the inference and experimental design methods described here, have the potential to lead to significantly improved connectivity estimates. Even if all the neuronal inputs are eventually observed, if the observation probability pobs is low then the variance due to the unobserved inputs may still be high, since, at any given time, most of the inputs to each neuron will be unobserved (see also [28]). As a result, the duration of the experiment required for accurate inference increases quadratically with the inverse of the observation probability (Eqs (18) – (19) and Fig 7), and weak weights become much harder to infer (Fig 6). Note this variance may be significantly reduced if we only aim to infer the input connections to only a few neurons (Fig 8). However, in many cases we wish to infer the entire network. In those cases the variance issue will persist, for any fixed observation strategy that does not take into account any prior information on the network connectivity. However, there might be a significant improvement in performance if we can focus the observations on synaptic connections which are more probable. This way, we can effectively reduce input noise from unobserved neurons, and improve the signal to noise ratio. As a simple example, suppose we know the network is divided into several disconnected components. In this case, we should scan each sub-network separately, i. e. , there is no point in interleaving spike observations from two disconnected sub-networks. How should one focus observations in the more general case, making use of past observations in an online manner? Again, we leave this “active learning” problem as an important direction for future research. In this work we suggest a “shotgun” experimental design, in which we infer the connectivity of a neural network from highly sub-sampled spike data. This is done in order to overcome experimental limitations stemming from the bounded scanning speed of any imaging device. To do this, we develop a statistical expected loglikelihood-based Bayesian method. This method formally captures the intuitive notion that empiric spike correlations and mean spike rates are approximately the sufficient statistics for connectivity inference. Exploiting these sufficient statistics, our method has two major advantages over previous related approaches: (1) it is orders of magnitude faster (2) it can be used even when the spike data is massively sub-sampled. We show that by using a double serial scanning scheme, all spike rates and correlations can be eventually inferred (and therefore neural connectivity). We demonstrate numerically that our method works efficiently in a simulated model with highly sub-sampled data and thousands of neurons. We conclude that the limited scanning speed of an imaging device recording neuronal activity is not a fundamental barrier which prevents consistent estimation of network connectivity.
Optical imaging of the activity in a neuronal network is limited by the scanning speed of the imaging device. Therefore, typically, only a small fixed part of the network is observed during the entire experiment. However, in such an experiment, it can be hard to infer from the observed activity patterns whether (1) a neuron A directly affects neuron B, or (2) another, unobserved neuron C affects both A and B. To deal with this issue, we propose a “shotgun” observation scheme, in which, at each time point, we observe a small changing subset of the neurons from the network. Consequently, many fewer neurons remain completely unobserved during the entire experiment, enabling us to eventually distinguish between cases (1) and (2) given sufficiently long experiments. Since previous inference algorithms cannot efficiently handle so many missing observations, we develop a scalable algorithm for data acquired using the shotgun observation scheme, in which only a small fraction of the neurons are observed in each time bin. Using this kind of simulated data, we show the algorithm is able to quickly infer connectivity in spiking recurrent networks with thousands of neurons.
Efficient "Shotgun" Inference of Neural Connectivity from Highly Sub-sampled Activity Data