BACKGROUND
connectomics, which seeks to identify the connectivity structure between all pairs of neurons not only in local circuits but over the entire brain, is a crucial research direction  <cit> , because the brain’s functions are believed to emerge within the connectivity structure of its constituent elements, the neurons.
toward such a direction, functional connectivity analyses provide a powerful bottom-up tool to investigate the neurophysiology of the relationships of multiple neurons .
moreover, a decoding study has quantified information amount that is encoded into functional connectivity of retinal ganglion cells  <cit> .
when we focus on information processing principles of local neuronal circuits, we need a sophisticated approach to unify top-down simulations and bottom-up experimental observations both in terms of functional connectivity and anatomical connectivity.
the reliable detection of functional connectivity is thus vital to identify the functions of local and global networks in the brain.
similar attempts have focused on the functional connectivity of the higher-order view of the brain: dynamic causal modeling of the brain regions  <cit> , for example.
such research is sometimes called macro-connectomics.
although our functional connectivity analysis method can be applied to such macro-connectome problems, in this study we focus on functional connectivity between neurons, one of the elementary levels of the brain’s hierarchy.
by definition, a pair of neurons is called functionally connected if the physiological activities of these neurons interact with each other.
granger causality provides a theoretical definition to a directed causal relationship that can be estimated from time course observations like spike trains from neuronal circuits  <cit> .
in contrast, structural connectivity implies anatomical observation of synaptic structures between the neurons found by electron microscopy or the imaging of synaptic markers  <cit> .
functional connectivity and structural connectivity have different theoretical implications.
structural connectivity does not necessarily imply functional connectivity; since most synapses are silent or inactive in vivo, few active synapses establish functional connectivity  <cit> .
functional connectivity does not necessarily require structural connectivity; since the observation cannot be complete, functional connectivity may be produced by indirect causality bypassed by unobservable neurons and glial cells and mediated by various transmission channels like electronic coupling and metabolic factors.
some studies have shown indirect evidence of relationships between functional and structural connections  <cit> .
because one major objective of the physiology of neuronal circuits is to determine their functions, functional connectivity analysis can be more important than its structural counterpart.
in order to estimate functional connectivity among many neurons, we require sufficient data of electrical activity of the neurons.
recent functional multi-neuron calcium imaging technique enabled us to observe hundreds or thousands of neurons’ activities simultaneously  <cit> .
however, there are following two trade-offs between information amounts that can affect the quality of functional connectivity estimation.
the first is a trade-off due to limited observation speed of a scanning microscopy; we can observe the larger number of neurons by the wider field of sight with the larger spatial resolution, leading to the lower temporal resolution.
the second is a trade off due to unavoidable photobleaching; higher signal to noise ratio requires stronger light emission, leading to stronger photobleaching of fluorescent dyes that shorten the observation time length.
thus, these two trade-offs limit information amount per neuron even when the number of simultaneously recorded neurons is large.
moreover, information amount per a pair of neurons is decreasing with growing number of neurons because the number of pairs of neurons is quadratic to the number of neurons.
we call this situation an information shortage for functional connectivity estimation.
in the current study, we focus on generalized linear models  of spike response functions  and their extensions, which are found powerful statistical tools to sketch the functional connectivity between multiple neurons .
an srf expresses the increase and decrease of the spike generation probability triggered by a pre-synaptic spike input as a function of the lag-time between pre-synaptic input and post-synaptic output spikes.
based on the glm framework, srfs for all possible pairs of neurons are estimated simultaneously, and the granger causality test substantially reduces false positives caused by indirect causality that cannot be distinguished if they are estimated for each local pair of neurons.
recently, several studies have presented sparse estimation methods that prefer zero connection weights for a large part of all neuron pairs to improve the estimation of glm-based srfs  <cit> .
on the other hand, there have been a few attempts of network structure inference without focusing on fitting glm, such as information theoretic analysis of granger causality  <cit>  and inference of dynamic bayesian networks  <cit> .
in case of information shortage, two major requirements to estimate srfs, model fitting and detection accuracy, are not compatible with each other.
estimation of glm, in case of information shortage, owes its performance much to regularization term, in which regularization coefficient is tuned for better model fitting.
recent sparse estimations of glm are not exceptions.
however, the best-tuned coefficient for model fitting does not generate the best results for detection accuracy; not only that there is no guarantee of compatibility, but also that there is a severe contradiction.
discrimination accuracy by the sparse estimation of srfs usually depends on the arbitrary tuning of a hyperparameter that controls the sparseness.
the hyperparameter tuning, however, involves severe problems that are not well understood.
namely, there is no guarantee that the best hyperparameter to minimize the fitting error in terms of cross-validation leads to the best discrimination accuracy to detect functional connectivity.
actually, the connectivity detection performance is not favorable when we tune the sparse model so as to maximize the likelihood.
moreover, no good method exists to control false positives with sparse estimation mainly because it is difficult to obtain closed-form null distributions of testing statistics.
there have been many studies focusing on model fitting, but few studies focusing on detection accuracy except for the likelihood-ratio testing of granger causality  <cit> .
however, they did not consider the case of information shortage and did not include any regularization factors.
the performance of their estimation might be sub-optimal, particularly with information shortage.
in this study, we propose a new functional connectivity analysis method focusing on improving the detection performance especially in the case of information shortage.
our method consists of two major procedures.
in the first, srfs are identified by estimating the parameters of glms with smooth bases to represent srfs and l <dig> regularization.
the second procedure is a granger causality analysis based on empirical bayesian testing.
because our situation for functional connectivity analysis is a typical multiple simultaneous testing, false positive control is definitely important.
empirical bayesian testing provides a reliable way to stably control the false positive proportion.
our research scope is summarized in table  <dig>  we focus on the case with information shortage in which the number of neurons is large but the observation length is relatively short.
because the amount of information available to estimate srf is proportional to the number of spikes per neuron pair, the statistical uncertainty is large in various imaging studies, like calcium imaging.
thus, we propose a methodology that is essentially different from best practice in electrophysiological studies.
namely, we apply multiple testing to determine the functional connectivity, rather than to estimate srf, for each pair of neurons.
the goal is to list a set of candidate connections while controlling false positives, rather than to build a single most likely model that maximizes the fitness or expected likelihood.
toward this goal, we apply a glm model with non-sparse regularization in contrast to one with sparse regularization.table  <dig> conceptual difference between two distinct typical cases of functional neuronal connectivity analyses
spike response estimation
generalized linear model of spike responses
in this section, we introduce a formulation of glm of spike responses according to  <cit> .
suppose that we have a spike train dataset, n={ni:i= <dig> …,c;t= <dig> …,t}, where ni denotes the number of spikes of the ith neuron  in the tth time-bin  with a common width over all the bins.
we assume the bin width is short enough so that the maximum spike number in a single bin is unity, and thus ni takes either one or zero.
we also assume that c neurons may receive a common input signal transmitted through l external channels: e={el∈r:l= <dig> …,l}, where el denotes a real-valued signal sent through the lth channel.
a stochastic spike response model represents the conditional probability to observe a spike of the ith neuron in the tth time bin pi≡pr=1|n,e,ri), where conditions n, e, and ri denote the past spike history up to time t- <dig> of all neurons c= <dig> …,c including the ith neuron itself, the external input at time t, and a set of parameters that define the response function of the ith neuron, respectively.
according to glm, the conditional probability is a function of total spike response λi,  <dig> pi=f), where the total spike response is given as a linear summation of all external inputs el and the past spikes from all neurons n,  <dig> λi≡λi,e,ri)=ri0+∑l=1lrileel+∑c=1c∑s=1mricnc.
ri <dig>  rile, and ric are respectively a background activity level of the ith neuron, a response coefficient to the lth external input, and the spike response with lag-time s to the cth neuron, of the ith neuron.
natural number m, called a window, is the maximum time-lag considered in the history.
f is a  link function from total response λi to probability pi, and we concentrate on a case of logistic link f=1/).
a reasonable alternative may be f=1-exp), called the complementary log–log link function.
the latter case makes this spike response model equivalent to a poisson spike model at a limit of infinitesimal bin-width  <cit> .
we may consider common external input e that effectively represents the total effect received by the observable neurons from unobserved ones and/or other external inputs.
if we can directly observe those external inputs, the estimation of the spike response model is reduced to an auto-regression  problem; by regarding glial activities as observable external inputs to a neuronal network, a glm-based spike response model for a neuron-glia system can be estimated  <cit> .
if we cannot observe the external inputs, the estimation is of an ar type with a moving average , in which we need to simultaneously estimate the external inputs and the response functions of individual neurons.
estimation and errors
we may obtain a small but non-zero estimation of response functions for some pairs of neurons that in fact have no functional connectivity because of the finite amount of available data.
such limited data causes statistical uncertainty preventing a clear-cut determination and two types of statistical errors, false negatives and false positives.
to deal with such statistical uncertainty in our application of functional connectivity analysis, we have several options in each of the following two steps that constitute functional connectivity analysis:  the estimation method of a glm-based spike response model including spike response ric between neurons i and c, and  the selection of statistics to test the null hypothesis, where the spike response ric is zero for any time-lag s if the null hypothesis is true, and ric is not zero for some s if the alternative hypothesis is true.
our ways of dealing with these two issues will be discussed in the following subsections.
application of smooth bases and regularization term
this subsection describes a way to estimate a glm-based spike response model.
we assume that each spike response function, ric, is represented as a linear combination of a small number of smooth bases:  <dig> ric=∑k=1kaickbk, where bk and aick denote the kth basis function that is shared by all neurons and the kth basis loading coefficient between neurons c and i, respectively.
to reflect the consistently positive  character of the facilitatory  post-synaptic current, epsc , and the consistent profile of epsc  induced by a single pre-synaptic spike, each basis function is given by the following gamma density function:  <dig> bk=ga.
here, ga denotes the probability density function of a gamma distribution whose mean and variance are set at mk=vk=12k <dig>  the gamma density function as a filter basis was first proposed in  <cit> .
a basis function with a small index k peaks at a small s value, which induces a large but short-delayed epsc or ipsc, and that with a large index is broad with a large mean value, which induces a small but long-lasting epsc or ipsc.
number of bases k is arbitrarily determined so that the glm model fits the real spike response well.
smoother response function is preferred in the estimation when a smaller number k is set, which can improve the estimation by preventing over-fitting to statistical uncertainty especially in cases with information shortage.
from eqs.
and , the total spike response of the ith neuron is given by  <dig> λi=ri0+∑l=1lreilel+∑c=1c∑k=1k∑s=1maickbknc.
when estimating the spike response model parameters, the following regularized log-likelihood function is independently maximized with respect to the unknown parameters ri <dig>  reil, and aick, l= <dig> …,l,c= <dig> …,c,k= <dig> …,k, for each neuron i:  <dig> li=∑t=1tlti-ηregi, where lti≡nilog)+)log)) is log-likelihood at time t, η is a hyperparameter that controls the strength of regularization, and regi is the regularization term.
we may use the following l <dig> regularization to encourage non-sparse estimation of the parameter aick:  <dig> l2:regi≡12∑c=1c∑k=1kaick <dig>  we alternatively consider the following l <dig> or group lasso regularization for facilitating sparse estimation:  <dig> l1:regi≡∑c=1c∑k=1kaick,   <dig> gl:regi≡∑c=1c∑k=1kaick <dig>
the common input term el in our glm-based spike response model  is fixed at an estimated value in a preprocess  prior to the model estimation.
the parameter estimation to maximize the regularized likelihood function  is implemented using a dual augmented lagrangian method  <cit> .
estimation of external input
when some common external inputs el,l= <dig> …,l were supposed to exist but not given, they were estimated using principal component analysis   <cit>  of the observed time course.
let us imagine a typical case that there are some synchronized spikes of many observed neurons.
in this case, we would detect facilitative connectivity between all neurons that emit the synchronized spikes.
however, we can reduce the effect of the synchronized spikes to the estimation of functional connectivity between the observed neurons if we regard the synchronized spikes as an effect of common external input.
we can extract the synchronized signal as the first principal component of a set of smoothed spike sequences of the neurons, or multiple signals as some components if there are some distinct sets of neurons that are synchronized.
note that we assume these common external inputs as a minimum device to reduce the harmful effect to functional connectivity estimation, rather than to infer any external reality.
we first applied a moving average filter of a certain window size to the spike time course ni of all the observed neurons i= <dig> …,c to obtain smoothed spike density profiles.
next, we calculated l principal components with eigenvalue decomposition of c×c covariance matrix for the smoothed spike density profiles of the c observed neurons and regarded them as the estimated l external inputs.
finally, we fit the glm model involving the common external inputs and assessed the fitness using akaike’s information criteria   <cit>  to determine the number of principal components l, where aic is defined as the log-likelihood value minus the number of free parameters of the glm.
if there is no external input, the smoothed spike density is based on the internal fluctuation stemming from the constant term  of the observed neurons.
then there is no significant principal component of the smoothed spike density profiles.
if there are principal components of the smoothed spike density profiles, they cannot be explained by the bias of the spontaneous activities of the observed neurons; therefore, they are based on the external inputs.
estimation of the external inputs by the above method corresponds to estimation of slow behaviors of the network and is useful for removing factors that are not able to be represented by our glm; it is beneficial to our main objectives of estimating the spike response functions and their significance measure in a stable manner.
statistical significance
here, we focus on statistical tests for functional connectivity analyses of glm-based spike response models.
our statistical tests are based on granger causality, but there are several options when applied to glm-based spike response models.
granger causality and a corresponding test statistic
granger causality is a criterion to determine the existence of a directed causal relationship between two nodes, from both of which we observed time sequences.
by definition, it is said that there is granger causality from nodes a to b if we can significantly better predict the future time sequence of node b using the information of the past and current time sequences of nodes a and b than only using that of node b. in the context of the functional connectivity analysis of our glm-based spike response models, the granger causality from the cth to the ith neurons can be examined by hypothesis testing with null hypothesis h0:ric= <dig> for all s= <dig> …,m against alternative hypothesis h1:ric≠ <dig> for some s= <dig> …,m.
kim et al.
<cit>  proposed a likelihood-ratio test of all h <dig> i,c= <dig> …,c by a log-likelihood-ratio statistic:  <dig> lric=- <dig>  where l0ic and li are the log-likelihood of the null and alternative hypotheses, respectively.
this statistic reflects the loss of model fitness by omitting the functional connectivity from neuron c to neuron i; the alternative likelihood li is first calculated for the post-synaptic neuron i and the null likelihood l0ic is calculated by omitting the contribution of each of the pre-synaptic neurons c= <dig> …,c.
in their likelihood-ratio test, the alternative log-likelihood is given by eq.
without the l <dig> regularization, that is, with η= <dig>  and the null log-likelihood is given by eq.
with ric= <dig> for all s= <dig> …,m .
note that the null log-likelihood depends on the pre-synaptic neuron index c, but the alternative log-likelihood does not, because the null likelihood has an additional constraint dependent on the pre-synaptic neuron index, which is not involved in the alternative likelihood.
we can calculate the p value of the likelihood-ratio test based on the fact that the test statistic lric obeys a chi-squared distribution with m degree of freedom in a large sample limit  <cit> .
considering the situation of multiple hypothesis testing, q values were also calculated based on the p values for all h <dig> i,c= <dig> …,c based on the standard procedure  <cit> .
kim et al.’s study  <cit>  followed the standard methodology of statistics; especially for a relatively small amount of data, however, parameter estimation without any regularization may suffer from failure or instability, producing an unreliable calculation of the likelihood-ratio.
one possible way to overcome this difficulty is to introduce regularization to the null and alternative likelihood functions, but then the test statistic no longer obeys the asymptotic chi-square distribution.
moreover, parameter estimation for the null likelihood must be performed for every pair of pre- and post-synaptic neurons, which is computationally expensive.
empirical bayesian testing
empirical bayesian testing uses empirical null samples in place of null distributions that are analytically defined like the one described in the previous sub-section.
in our particular application, we do not know the null distribution analytically because of the introduction of regularization and kernel functions to the glm-based spike response models; in this case, empirical bayesian testing is a practical choice.
in our empirical bayesian testing, a certain number of empirical null samples of test statistics are obtained by the following surrogation method.
we assume several surrogate neurons, each of which emits an artificial spike train, by time-shifting a real neuron’s spike train: nc∗=nc for some c∈ <dig> …,c.
here, we also assume that there is no prominent periodic activity in the network and set the time-shift ts to a number larger than the time-lag m. then, the spike train of a surrogate neuron c∗ becomes independent from that of any real neuron i= <dig> …,c, in the time lag.
after adding a certain number of surrogate neurons to the set of real neurons, we estimate the glm-based spike response model.
in this estimated spike response model, the test statistic between a real neuron a and a surrogate neuron c, which was generated by time-shifting the spike train of a real neuron other than a, can be regarded as a null test statistic, because a and c are independent, with no functional connectivity between them.
note that there is a non-zero probability to include some statistic values from the non-independent  pairs of neurons in the surrogating process, which can weaken detection power.
however, the risk is small because the strength of ordinary response is negligible for a time lag larger than a certain limit and the risk can be smaller by setting large enough time-shift for the surrogate.
more importantly, inclusion of false negatives does never violate conservativeness.
in multiple simultaneous hypothesis testing like functional connectivity analysis, q value, which is an estimation of false discovery rate, is often used as an alternative significance measure to a p value  <cit> .
when test statistics are available along with empirical null samples for calculating them, q values can be estimated directly from them without calculating the p values.
in our empirical bayesian multiple simultaneous hypothesis testing, we used this direct method to obtain q values for all the multiple hypotheses.
this procedure is described in the additional file  <dig>
shape-related statistics
according to empirical bayesian testing, we can use an arbitrary test statistic.
we also examined the following statistics since they are easier to compute compared to likelihood-ratio statistic after estimating the glm-based spike response model, but they still represent the non-zero character of the spike response function.
we prepared several statistics to characterize the spike response function  between neurons c and i.
surface sicsurface≡∑s=1m|ric| peak sicpeak≡maxs|ric| parameter vector norm with a special metric design  we evaluate the norm of parameter vector aic· with a special metric design Σkk′.
sicmd≡∑k,k′aickaick′Σkk′, where Σkk′ is the -th element of a regularized inverse of the following covariance matrix s <dig> of the parameter vectors of null links s0∝∑∈h0aic·taic·.
here, aic· is a k-dimensional parameter vector whose k-th element is aick and t denotes a transpose.
in the equation above, the summation was taken over every pair of pre- and post-synaptic neurons whose functional connectivity was determined not to exist by the statistical test.
thus, the md incorporates the anisotropic nature of the null distribution and evaluates the distance between each case and the typical null cases.
delay sicdelay≡argmaxs|ric|
the last term, delay, is a bit different, because it is used for describing a functional connection detected by a certain statistical test, rather than for defining a test statistic to detect the functional connectivity.
statistical tests for functional connectivity analysis
in this study, we compared the following statistical tests in a scenario of functional connectivity analysis based on a glm-based spike response model.
they examine a pairwise functional connectivity, that is, whether a direct relationship exists from a pre-synaptic neuron c to a post-synaptic neuron i.
chi <dig> if there is no regularization term, a chi-square test can be applied to likelihood-ratio statistic lric, because it asymptotically obeys a chi-square distribution with m degree of freedom.
then, we can approximately use a chi-square test by setting the regularization hyperparameter η to a small positive value to avoid a divergence in the optimization procedure.
a p value is simply obtained by integrating the chi-square distribution.
a q value is estimated using the p values of all null hypotheses h <dig>  based on a standard procedure.
dof a p value and hence a q value were obtained by the same procedure as in chi <dig> above, except that the degrees of freedom were determined to fit an empirical null distribution of the likelihood-ratio statistic lric.
here, the empirical null distribution was produced by the surrogation method described in “shape-related statistics” section.
eb empirical bayesian testing was applied to the likelihood-ratio statistic after the model parameters were estimated by the optimization with a larger regularization hyperparameter η. by comparing the results of eb and chi <dig> , we can see the contribution of the l <dig> regularization and employing smooth bases.
eb-arbitrary stat because empirical bayesian testing can be applied to an arbitrary test statistic, the statistics described in “shape-related statistics” section were also examined in its framework.
note that we did not apply empirical bayesian testing for estimation results with l <dig> and group lasso regularization, because the empirical null distributions of the statistics in “shape-related statistics” section were found to be quite different among neurons, which was not the case with the l <dig> regularization.
confidence interval  and maxz a 95 % wald confidence interval is calculated for each estimated value of glm weight wk, wk± <dig> σkk, where σkk is an asymptotic standard deviation that is derived using hessian of the objective function  <cit> .
according to  <cit> , directed connectivity is detected if any ci corresponding to the connectivity is significantly different from zero.
we also consider the z-score of a weight value zk=|wk|/σkk and that of the response function z=∑k=1kzkbk.
we integrate them into a maxz statistic, a maximum of z,s= <dig> …,m corresponding to the directed connection, so that a criteria, maxz> <dig> , is equivalent to the application of ci in  <cit> .
as the straight application of the 95 % ci is no longer a valid multiple testing, we utilized our eb-arbitrary stat framework in order to include the maxz statistic in the following comparison study.
RESULTS
experiment 1: poisson model simulation
here, we examined the performance of the functional connectivity analysis methods using a dataset generated by simulating a network of glm-based spike response models.
we designed an artificial neuronal network of  <dig> neurons  consisting of three groups of neurons, g <dig>  g <dig>  and g <dig>  each of which consisted of five neurons.
this is a recurrent network: g <dig> → g <dig> → g <dig> → g <dig>  from g <dig> to g <dig>  there are ten excitatory connections among  <dig> pairs of pre- and post-synaptic neurons.
similarly from g <dig> to g <dig>  there are ten excitatory connections among  <dig> pairs, and from g <dig> to g <dig>  there are ten inhibitory connections among  <dig> pairs.
there is no connection between neurons belonging to the same group.
the activity of each neuron was simply determined by a glm-based spike response model, eqs.
and , with the complementary log-log link function that simulates a poisson spike model.
each neuron also received an external white poisson input of around 50–150 hz independently from the external input to the other neurons.
this assumption that the external input was independent between neurons is slightly different from the one assumed by our glm model, in which the external input is shared by the observed neurons.
however, as this situation of independent poisson inputs is much simpler, our glm-based spike response model can deal with it by just removing the estimation procedure of the external input .
the activities of all  <dig> neurons were observed at a sampling frequency of 200 hz.
with these settings, we obtained from  <dig>  to  <dig> spikes on average per neuron per 1 s .fig.
<dig>
a simulated neural network that was used in experiment  <dig>  three groups of neurons g <dig>  g <dig>  and g <dig> construct a recurrent network: g <dig> → g <dig> → g <dig> → g <dig>  gray scale blocks in the matrix denote a direct functional connectivity between source  and destination  neurons.
b–g performance comparison between glms based on sparse estimation with group lasso  and on l <dig> regularization .
comparison between logistic regression  and poisson regression  is also performed.
the glm parameters were estimated with various settings of hyperparameters and the best hyperparameter value for each case was determined so as to maximize likelihood on a validation dataset.
the results were evaluated with sensitivity and false positive proportion calculated by using a true directed link set for the simulation.
experiments with short  and long  datasets are shown in the upper  and lower  sets of panels, respectively.
b, e roc curves are drawn for each case at the best hyperparameter that maximizes likelihood on the validation data.
c, f auc score for each setting of regularization hyperparameter.
markers denote the values at the best hyperparameter that maximizes likelihood on the validation data.
d, g accuracy of model fitting, measured by likelihood on the validation data, for each setting of regularization hyperparameter.
markers denote the values at the best hyperparameter that optimizes model fitting.
we found that the differences were small between logistic and poisson regression models.
roc curves or auc were not necessarily the largest at the best-tuned values of hyperparameters for each of the four cases.
sparse estimation with the group lasso  was more sensitive to the hyperparameter settings than the l <dig> regularization
first, we compared results of the sparse and non-sparse estimation methods for glm.
we selected spike train data of  <dig> and  <dig>  samples from the simulated data for training and  <dig>  samples for validation.
the sample numbers in both cases are considerably shorter than some typical cases; as a reference, kim et al.
<cit>  used  <dig>  samples to demonstrate their method.
when using  <dig>  samples, we found our method achieved 100 % accuracy, although the results are not shown here.
we compared two types of regularization, group lasso  and l <dig>  for glm, corresponding to sparse and non-sparse estimation, respectively.
we also used l <dig> regularization, but did not show its results because the performance with l <dig> was comparable to that of group lasso.
we also compared two link functions, logistic  and poisson ; here poisson stands for the complementary log-log link function that simulates poisson spikes.
the srfs were estimated using the training data with various hyperparameters that determine the strength of regularization.
the estimated srfs were evaluated with two criteria, detection accuracy of the true functional connectivity and fitness to unseen data; detection accuracy is shown as a receiver operating characteristic  curve and the area under the curve  score, whereas fitness is given by the log-likelihood for the  <dig>  samples for validation.
in sparse estimation like with group lasso regularization, some statistics like peak may be exactly zero for some connections.
although a naive threshold to detect functional connectivity based on such statistics would be zero, setting a non-zero threshold can in fact produce better functional connectivity estimations.
therefore, we evaluated roc and auc by changing the threshold value even for sparse estimation methods.
figure  <dig> shows the results.
in all aspects, logistic and poisson link functions did not show any significant differences whereas the difference between sparse and non-sparse estimation was fairly large.
when the hyperparameter was determined to increase the fitness to the validation data, the connectivity detection performance was poor as shown in the roc curves  and .
as shown in the panels  and , the best hyperparameters in terms of the detection accuracy or auc were smaller than those for the largest fitness .
the auc could be made larger by arbitrarily setting the hyperparameter to a smaller value.
with respect to the fitness, the best model was achieved by group lasso in both sizes of training data,  <dig> and  <dig>  as shown in  and , respectively.
with respect to the detection accuracy, the best auc with the l <dig> regularization  was better than that with group lasso  for t= <dig>  and the best aucs were  <dig>  with both l <dig> and group lasso for t= <dig>  as shown in  and , respectively.
more importantly, l <dig> regularization performed well for a wide range of hyperparameters, while group lasso with non-optimal hyperparameters performed noticeably worse than the best performance.
we then compared the results obtained by the three statistical tests, chi <dig>  dof, and eb, described in “statistical tests for functional connectivity analysis” section.
the glm-based spike response model was estimated to maximize the l <dig> regularized log-likelihood .
when performing the chi <dig> and dof, the regularization hyperparameter was set at η= <dig> .
when performing the eb, we compared two cases eb-lr and eb-regularized-lr, in which the regularization hyperparameter was set to η= <dig>  and η= <dig> , respectively.
the value η= <dig>  was arbitrarily determined to obtain similar results to those of η= <dig> but with computational stability.
the value η= <dig>  was also arbitrarily set to a small positive value, because there was no predominant way to set it for enlarging the detection power.
we also used smooth bases in the eb-regularized-lr.
figures  <dig> and  <dig> show the results.
the upper panels show the roc curves with auc scores in the title, in both of which the higher and the larger are the better ones, respectively.
in the lower panels, the false discovery proportion  is shown for each q value threshold on the horizontal axis; fdp is the proportion of false positives  in such links that were detected as functionally connected because the corresponding q values were smaller than the threshold on the horizontal axis.
the lower panels evaluate the quality of the q value estimation.
because a q value is defined as the fdp estimation, when the fdp is smaller than the q value threshold, the corresponding q value estimation is said to be conservative.
from the upper panels of fig.
<dig>  the roc curves and the auc scores of the proposed method  consistently outperformed those of its un-regularized version , and such a benefit of regularization was prominent especially when the observation length t was small.
in the lower panels of fig.
<dig>  the conservativeness of the q value estimation was compared among eb-lr , chi <dig> , dof , and the proposed eb-regularized-lr .
if a line is close to the diagonal x=y line, the corresponding q value estimation is faithful, and if a line is below the diagonal line, the q value estimation is conservative.
we prefer q value estimation that is faithful but conservative; that is, we do not prefer a line that goes over the diagonal line.
we found that the empirical bayesian q values  were stably good, while a classic chi-square test without dof adjustment  almost failed to perform false positive control when the observation length t was small.
when the observation length t was large enough, the results depicted by the blue lines were consistent with those reported by  <cit> .
in the lower panels, the fdp values sometimes went over the diagonal line, violating the conservativeness.
this violation comes from the estimation variance of q values; a large variance in the q value estimation may arise especially when there is significant correlation between hypotheses, which is the case in the functional connectivity analysis of neuronal networks.fig.
<dig> results of experiment  <dig> in which total error and fdp control were compared using several statistical testing methods.
upper panels show receiver operating characteristic  curves, where horizontal and vertical axes denote the specificity and sensitivity of the statistical test, respectively.
red and blue
lines show the eb-regularized-lr  and eb-lr  results, respectively.
in the title of each column, observation length t and the area under the roc curve  scores of the two methods are shown.
lower panels display false positive control, where horizontal and vertical axes denote q value threshold and false discovery proportion  when the q value threshold on the horizontal axis was used, respectively.
because the q value threshold is regarded as an estimation of the fdp value, when fdp is located under and over the thin black diagonal line of x=y, the corresponding q value estimations are considered conservative and aggressive, respectively.
a conservative line is preferable, because by definition, the q value should be a conservative estimation of the fdp value.
if the line is close to the diagonal line, the q value estimation is faithful.
blue broken, dotted, and solid lines represent the q values estimated by chi2-lr, dof-lr, and eb-lr for the likelihood-ratio statistics without regularization  is the likelihood-ratio statistic with larger regularization.
the six panels from the leftmost to the rightmost columns correspond to the experimental settings of observation lengths t= <dig> , <dig> , <dig> , <dig>  and  <dig>
fig.
<dig> results of experiment  <dig> in which total error and false positive control were compared among four test statistics: peak , maxz , surface , and md  .
see caption of fig.
<dig> for explanation of horizontal and vertical axes of upper and lower panels
fig.
<dig> results of experiment  <dig> in which spike trains ware observed from a recurrent network of  <dig> non-linear conductance-based  neurons.
notations are identical as in fig.
2
figure  <dig> shows the results of statistical tests employing shape-related statistics .
we did not find any substantial differences either in total detection accuracy  or false positive control  among the statistics of the peak , the maxz , the surface , and the md .
it is noteworthy that we numerically confirmed the stable false discovery rate control regardless of possibly different null distributions of various test statistics.
this result may allow us to select arbitrary test statistics in the proposed empirical bayesian testing framework.
we may expect good roc by selecting a relevant test statistic that fits particular shapes of true spike response functions although the rocs were similar between different statistics in the current case partly because we assumed simple spike response function shape such that caused no contradiction between test statistics that we compared here.
the roc curve of maxz statistic included 95 % ci criteria as a special case , which is specified as a marker ’o’ on the curve .
we found that the 95 % ci provided a moderate balance between sensitivity and specificity in these cases.
however, note that the 95 % ci had no guarantee of conservativeness that eb could provide as q value.
experiment 2: a network of non-linear conductance-based neurons
next, we applied functional connectivity analysis to a model mismatch case.
a spike train dataset was obtained by simulating a recurrent network of non-linear conductance-based neurons.
each neuron was modeled by single-compartment hodgkin–huxley  equations  <cit> , consisting of sodium, potassium, and leak channels.
thus, the neuronal activities did not follow glms.
we, however, expect that the proposed false positive control procedure works reasonably well because our framework mainly depends on the empirical null distribution of test statistics built upon a set of surrogate neuron pairs and hence does not depend on the glm’s fitness.
the simulated recurrent network resembled that in experiment 1: three layers, each of which had five hh-type neurons.
the simulation was conducted using the nest simulator  <cit> , and the parameters of the hh neurons were set by the default setting of the nest simulator.
the results are shown in fig.
<dig>  in which similar behaviors of the statistical tests can be seen to those in experiment  <dig>  the q value estimation by chi-square tests, chi2-lr  and dof-lr , was substantially poor even when observation length t was large.
when employing non-linear neurons like the hh-type ones, glm could no longer well represent their non-linear behaviors, and the model of the alternative hypotheses would contain bias, which was not resolved even after very long observation.
although the glm-based spike response model had bias, the empirical bayesian testing with regularized model estimation, eb-regularized-lr , showed reasonably good roc curves and auc values  and faithful estimation of the q values , based on an empirical null distribution constructed by the surrogation method.fig.
<dig> results of experiment  <dig> that considers partial observation with an external input.
a external input signals.
the artificially designed and an estimated signals are shown as blue and red lines, respectively.
<dig> time samples out of  <dig>  time samples are shown here.
b model selection via aic criterion.
vertical and horizontal axes denote aic criterion and number of external inputs that are involved in the glm.
aic values are normalized so as to that with no external input is zero.
c the estimated spike response functions.
five panels correspond to the following categories of neuron pairs;  not connected,  directly connected,  connected via one interneurons,  connected via two or more interneurons, and  connection from/to surrogate neurons.
red line denotes statistically significant directed neuron pairs of q< <dig>
fig.
<dig> results of experiment  <dig> that assessed functional connectivity detection performances at various spike frequency and connectivity weights with fixing time sample number at t= <dig>  spike frequencies and connectivity weights are controlled by parameters a and w, respectively.
rows and columns of the  <dig> panels denote variations of a and w values, respectively.
in each panel, we compared four roc curves; red and blue line color denote estimation methods of glm, l <dig> regularization and variational bayes estimation, respectively.
real and dotted line styles denote peak and maxz statistics, respectively.
the combination of glm with l <dig> regularization and peak statistic performed the best for most of the  <dig> panels
fig.
<dig> analysis of calcium imaging data.
a averaged calcium fluorescent image of a rat ca <dig> region.
b  <dig> neurons  and  <dig> probable directed connections  between the  <dig> neurons.
red and blue arrows denote putative facilitatory and suppressive connections, respectively.
statistical properties of the empirical bayesian testing
fig.
<dig> statistical properties of the empirical bayesian testing.
a the estimated distribution  and surrogate null distribution  of the test statistic, peak.
the blue curve denotes the estimated total distribution, which is a normalized histogram of the estimated statistic of all possible pairs of observed neurons.
the green curve denotes the estimated null distribution, which is a normalized histogram of the estimated statistics of the pairs of the observed and surrogate neurons.
b the estimated local false discovery rate.
vertical and horizontal axes denote local false discovery rate and the test statistic .
each blue point corresponds to a pair of source and destination neurons.
c a matrix of the estimated test statistic, peak.
rows and columns denote destination and source neurons, where the source neurons include ten surrogate neurons.
the color code denotes the value of the test statistic, peak
fig.
<dig> spike response functions of  <dig> top significant functional connections, arranged from the top-left to bottom-right panels in ascending q value order.
in each panel, a pair of source and destination neurons’ code numbers and the corresponding q value are shown.
horizontal and vertical axes denote the lag-time s  and the spike response function, respectively, where  <dig> sample corresponds to 10 ms
fig.
<dig> spike response functions of six pairs of significant bi-directional connections.
in each panel, both directed spike response functions are shown; for 0<s< <dig>  each spike response function represents ric, while for -80<s< <dig>  it represents rci, instead.
from the definition of spike response functions ric≈rci, which corresponds to the correlation of neurons i and c
fig.
<dig> spike response function and epsp possibly mediated by ampa and nmda receptors.
a a mixture of fast and slow epsps, represented by two gamma distributions peaking at  <dig> and 50 ms, respectively.
according to an electrophysiological study, rat ca <dig> pyramidal neurons showed sharp, strong and relatively slow epsps peaking at around 10 ms and lying over 10–100 ms, respectively  <cit> , which were possibly mediated by ampa and nmda receptors.
in the mixture, the amount of the slow epsp was set as small as one third of the blue broken line.
b the delayed response of the calcium indicator induced by the membrane voltage was represented by a convolution kernel, an alpha function of the delay time of 50 ms. c by convolving the convolution kernel b to the typical epsp , we have a typical response function obtained by our functional connectivity analysis.
contributions in the response function from ampa- and nmda-mediated epsps are displayed as blue lines
experiment 3: partial observation case
we applied the proposed method to a partial observation case involving external input.
we generated an artificial spike train by a simulation of  <dig>  time samples using the same artificial neural network as the ones in experiments  <dig> and  <dig> again.
in the simulation, we added an artificial external input signal  to all neurons, where connection weight from the external input to the neurons were set at random positive values.
in order to consider partial observation, we regarded  <dig> out of  <dig> neurons were considered as non-observed neurons, and spike response functions between  <dig> observed neurons, namely those of  <dig> times  <dig> pairs of pre- and post-synaptic neurons, were estimated by our procedure.
the external input was estimated by applying moving average filter of window size  <dig> and pca.
in fig 5a, we show the estimated external input denoted as red line.
we compared aic values of seven models with  <dig> ,..., and  <dig> external inputs, and found that the model with one external input was the best .
the estimated spike response functions are shown in fig 5c, where those function shapes corresponding to the  <dig> times  <dig> directed neuron pairs were divided into the following four cases:  not connected,  directly connected,  connected via one interneuron, and  connected via two or more interneurons.
in addition, we also show spike response functions corresponding to a pair of a surrogate and any observed neurons, that is called a category .
the spike response functions were estimated as significant  for a large part of direct connection, red lines in case  in fig 5c.
on the other hand, the estimated spike response functions of indirect connections, cases  and , were very weak and difficult to be distinguished from those of cases  and .
the results were identical when the number of principal component was set at zero, two, and more .
experiment 4: spike frequency, connection weight, and detection accuracy
how do different spike frequency and connection weight affect the detection accuracy in case of information shortage?
does the proposed framework keep its advantage against sparse estimation methods in these cases?
we generated an artificial spike train of  <dig> time samples by a simulation of a simple network involving three neurons,  <dig>   <dig>  and  <dig>  the network had two directed facilitative connections from  <dig> to  <dig> and from  <dig> to  <dig>  whose connection weights, or peaks of true spike response functions, were  <dig>  and w, respectively.
here we compared six values of w in { <dig> .
<dig> .
<dig> .
<dig> .
<dig> }.
all the three neurons had the same background activity level a, that was ri <dig> in equation , which controls total spike frequency.
we compared six values of a in { <dig> , <dig> , <dig> , <dig> , <dig> , <dig> }.
as the network size is small, we repeated the same simulation  <dig> times with different random seeds and averaged the performances in order to draw stable roc curves.
we compared two glm estimation methods with the l <dig> regularization with an arbitrary setting of λ= <dig>  and with the variational bayesian  estimation  <cit> .
vb approximates a hierarchical bayesian model with an automatic relevance determination  technique that determines sparse solution so as to optimize model fitting.
for test statistics, we compared peak, as a representative of the simplest cases, and maxz, as a case that is equivalent to the application of 95 % confidence interval  <cit> .
figure  <dig> shows the results of experiment  <dig>  the  <dig> panels correspond to the 6× <dig> variations of simulation setting parameters, w and a. four roc curves in each panel illustrate detection powers of the connection from  <dig> to  <dig> of weight w by the four methods, l <dig> with peak, l <dig> with maxz, vb with peak, and vb with maxz.
in total, we found that l <dig> prominently outperformed vb in most of the cases.
we also found that larger a tend to improve the detection accuracy as we expected.
interestingly, change of connection weight w did not affect much to the detection accuracy except for the case showing chance level w= <dig>  test statistic peak performed a little bit better than maxz indicating non-omnipotency of the confidence interval criteria although we could not conclude general superiority of peak statistic, neither.
note that the prominent performance of l <dig> against those of vb is consistent to the result of experiment  <dig>  it, however, does not stand for the general inferiority of vb.
further discussion with general comparison will be found in “alternative options in smoothing and regularization” section.
analysis of calcium imaging data
we applied functional connectivity analysis with glm-based spike response models to a real calcium imaging dataset.
a high-throughput calcium imaging experiment was performed on a cultured slice of the ca <dig> region of rat hippocampus  <cit> .
using the nipkow disc, the sampling frequency was as high as 100 hz, and the number of time samples was t= <dig> ; that is, the observation time was 10 min.
during preprocessing,  <dig> regions of interest , each of which showed its individual activity, and the corresponding time-courses of the fluorescence level were obtained.
a visual inspection found that each roi mostly corresponded to a neuron.
then, for each fluorescence time-course, we detected the peak timings by applying a fixed template that represented a calcium profile associated with a single spike.
a single detected peak may not correspond to a single real neuronal action potential; actually, there were some cases in which only one spike was detected for burst-like calcium activity.
we simply regarded such burst activity as a single spike because their number was small.
more detail on the above preprocessing is shown in the additional file  <dig>  we selected  <dig> neurons out of the  <dig> rois that showed frequent spike activities and found that the highest spike frequency was  <dig> times larger than the lowest one within them.
after this preprocessing, we evaluated the functional connectivity between the  <dig> neurons and their statistical significance using our method empirical bayesian testing after the regularized parameter estimation of the glm-based spike response model.
for each neuron pair, the spike response function of at most an  <dig> samples  time-lag was estimated.
we prepared ten surrogate neurons and assumed two channels of external inputs, which were estimated by applying pca analysis  to all of the  <dig> neurons’ activities.
we found that the false positive control based on shape-related statistics worked as well as the likelihood test statistic, so the q value threshold was determined using the peak statistic .
figure  <dig> shows the result of the functional connectivity analysis, in which  <dig> probable directed connections are shown by arrows.
the number of connections  was determined so that the estimated false discovery rate was  <dig> .
in this figure, some neuron pairs connected by red arrows  and being located close to one another may be single neurons segmented into two parts , such as soma and axon hillocks.
even if such pairs are removed, the detected functional connectivity still tends to connect two neurons near each other.
there were five neuron pairs connected by blue arrows  being located distant to each other.
statistical properties of the proposed empirical bayesian testing for this dataset are shown in fig.
<dig>  the probability distribution of a test statistic is a mixture  of unknown null and alternative distributions, which would correspond to the left mass and the right tail of the density function shape in panel , respectively.
the surrogate distribution of the test statistic  fits well to the left mass part of the empirical distribution , indicating that the surrogate distribution well simulated the unknown null distribution.
the local false discovery rate  shows that the estimated total proportion of null connections was π0= <dig>  and that the local false discovery rate was smaller than  <dig>  when the statistic was larger than  <dig> .
the significant pairs of neurons are shown as white dots in the source/destination matrix , so that their degrees of significance are seen by comparing to the background noise level in the surrogate area of the matrix.
the spike response functions of significant functional connectivity are shown in fig.
<dig> in ascending order of q values.
many response functions of the top significant functional connectivity had a positive peak at lag-times of around  <dig> samples , and some response functions included a tail around the 50th sample , which can be consistently considered as epsps possibly mediated by ampa and nmda receptors .
moreover, obvious continuity of response functions between pre-post connections, ric, and their reverse ones, rci, was found for many pairs of neurons, , in fig.
<dig>  this continuity  occurs partly because it reflects the correlation of neurons i and c. such correlation might have been overestimated due to jitter in spike timing detection; in some cases, neuron i emits a spike slightly earlier than that of neuron c, and in other cases, neuron c’s spike can be earlier, because of the jitter in the binning.
discussion
information shortage and detection performance
to detect functional connectivity, we need sufficient events of delayed spike pairs at the corresponding pair of neurons.
moreover, the number of events should exceed the background level of spontaneous activity of the post-synaptic neuron.
thus, various factors can reduce detection performance via information shortage; shorter observation time length, lower total spike frequency, weaker connection weight, an unknown external input to post-synaptic neuron, high background activity of post-synaptic neuron, and so on.
any of these factors unavoidably exists in real situations and consideration of information shortage is always vital for any case even when the observation time length is fairly long.
this point of view motivated the current study.
and, a part of them were demonstrated in experiments  <dig> and  <dig> in this paper.
alternative options in smoothing and regularization
recently, two distinct approaches have been proposed for functional connectivity analysis from the neuronal spike train data.
in the first approach, hierarchical bayesian methods were applied to glm-based models to obtain sparse solutions  <cit> .
a prior in the hierarchical bayesian setting allows the statistical estimation of glm-based models to prefer sparsely connected networks, such as reflecting the sparse structure of neuronal networks.
the sparseness also reduces the complexity of each spike response function.
the second approach was based on granger causality with a likelihood-ratio test  <cit> .
functional connectivity analysis based on granger causality can be performed for any regression model; the authors used the simplest ar model without any smoothing kernel or regularization to enable the null distribution of the likelihood-ratio test to be analytically calculated.
these two approaches are based on different statistical approaches but share the same objective of functional connectivity analysis; combining these approaches has options, such as, how to introduce regularization and how to define null and alternative hypotheses in the likelihood-ratio test.
in this study, we evaluated several reasonable options in terms of two criteria in statistical testing:  appropriate false positive control, which is based on good estimation of the false discovery rate in multiple simultaneous testing, and  high detection accuracy which is determined by increasing the number of true positives while fixing the level of false positives.
further comparison to chen et al.
<cit>  may better illustrate the difference between the two approaches.
their work proposed two distinct ideas, automatic-relevance-determination  with variational bayes  to estimate connection weights of glm and a use of 95 % confidence interval  criteria to determine the connection.
their ard with vb technique estimates sparse model as like the l <dig> regularization does with optimizing sparseness hyperparameter, and thus the result tends to be as good as the best among the l <dig> of various regularization parameters.
the use of the 95 % ci may resemble our work in the consideration of false positive risk control.
however, note that both of these ideas were designed based on the first approach for the purpose of fitting a generative model to the observed data rather than to detect functional connectivity.
the best model, measured by likelihood, is not necessarily the best model, measured by roc curve.
furthermore, the gap between these two criteria tends to be larger at information shortage.
we showed this contradiction in fig  <dig> of experiment  <dig> and more prominently in experiment  <dig>  our approach remarkably outperformed the vb if the performance is measured by roc.
consequently, we recommend choosing appropriate approach with determining purposes in advance to the analysis.
model mismatch bias
model mismatch bias is an important issue that needs to be addressed here.
any bayesian methodologies using prior probability to improve estimation performance always, at least implicitly, presume that the parametric model includes the truth.
although such a presumption can actually improve the estimation performance in ideal situations without any model mismatch, most models include a mismatch to the reality to a greater or lesser extent.
there is no guarantee that the glm represents actual neural behaviors well enough even with the best setting of the model parameters.
however, our framework is based on statistical testing that allows model mismatch to some extent.
namely, the alternative model does not necessarily include the truth because the functional connectivity is determined by rejection of the corresponding null hypothesis.
better and worse alternative models lead to decrease and increase in false negative rate, but do not affect false positive control.
the false positive is well controlled if the surrogate distribution well simulates the null process, that is, spike response between functionally unconnected neurons in our case.
experiment  <dig>  showed a reasonable performance in a typical model mismatch case where the data were generated by the hodgkin–huxley equations and analysis was done based on the glm model.
the analysis of calcium imaging data  suggested that the surrogate distribution well simulated the null part of the distribution of the peak statistic.
we thus conclude that our method can perform conservative statistical testing in many practical cases including that of model mismatch.
computational cost
computational cost for the glm parameter estimation is proportional to the observation length and to the square of the number of neurons.
when we add surrogate neurons for the empirical bayesian testing, it increases the effective number of possible connections and causes a corresponding increase in computational cost.
note that the number of surrogate neurons cs may not be so large when the real neuron number c is large because the variance in the statistical test depends on the number of null samples of their inter-connections, csc.
we set cs so that the number of null samples was around csc= <dig> in our experiments.
partial observation
partial observation is an important practical issue because it is difficult for many physiological techniques to cover all neurons recruited by the target system; functional connectivity analyzes should cope with this issue.
in experiment  <dig>  we demonstrated that spike response functions of indirectly connected pairs of neurons were far weaker than those of directly connected pairs.
we conducted similar experiments with various other settings and found similar results , where the variants included different observation length,  <dig> ,  <dig> , and  <dig> , different sets of observed and non-observed neurons, different moving average windows in pca analysis, different strengths of external inputs, and different network structures.
we consistently observed the following tendencies:when a strong external input exists, it is detected by the proposed procedure with pca and aic.
the proposed procedure with pca with aic hardly ever detected a total activity of non-observed neurons as a distinct external input.
the effect of non-observed neurons was indistinguishable to a background activity level of each observed neuron.
spike response was very weak and not significant for indirect connections.
spike response of direct connection was stronger than that of indirect connections but weaker than that in fully observed cases.
in analysis of partial observation case with the glm model, the total effect of non-observed neurons is roughly summed up into the scalar value of background activity level, which causes severe loss of information compared to the fully observed case.
this information loss has probably caused the higher estimation variance of spike response and the lower sensitivity to detect direct and indirect functional connectivities.
in total, experiment  <dig> demonstrated an additional reason to consider information shortage that might be present in case of long observation length and to need the proposed procedure with empirical bayesian test with surrogate neurons to control false positives.
further analysis would be needed as future studies to clarify the effect of non-observed neurons to the estimation variance although it is out of scope of the current study.
biological suggestions
when we applied the proposed method to the rat ca <dig> calcium imaging dataset, we extracted  <dig> directed connections .
many of their response functions exhibited sharp peaks at around  <dig> samples , while some  were also accompanied by dull and weak tails lying over 10– <dig> samples  .
the delayed activation corresponding to the weak tails may reflect sequential spike chains across multiple neurons, which prevail in hippocampal networks and often last for more than 100 ms  <cit> .
more physiologically, the sharp peaks and weak tails may reflect epsp of rat hippocampal pyramidal neurons revealed by an electrophysiological study; reportedly, epsp of rat ca <dig> pyramidal neurons showed sharp strong peaks at around 10 ms and relatively weak tails of 10–100 ms delay  <cit> , which were putatively mediated by ampa and nmda receptors.
in fact, when we convolved an alpha function with a delay time of 50 ms , which represents a chelating property of our calcium indicator , to the mixture of the fast ampa-mediated epsp and the late nmda-mediated epsp , we could reproduce a typical shape of our response functions .
changes in the mixture of ampa- and nmda-mediated epsps would produce various kinds of response functions.
therefore, the variety in the response functions found by our data-driven functional connectivity analysis  may have reflected to some extent the difference in the receptor distributions between the pyramidal neurons in the rat ca <dig> circuit.
there were two prominent negative responses, 16→ <dig> and 18→ <dig>  that also included positive peaks .
the positive peak of 18→ <dig> was weak and could not be statistically significant, while that of 16→ <dig> was still prominent.
seeing the response function of the counter direction 46→ <dig>  the positive peak in 16→ <dig> was considered to be an artifact coming from the significant positive peak in 46→ <dig> .
there were no other bi-polar responses that were significantly strong enough.
CONCLUSIONS
in this study, we presented a new combination of glm-based spike response modeling and empirical bayesian testing to perform functional connectivity analysis between neurons.
even when the observation period was relatively short, our method showed reasonably good detection accuracy while keeping good false positive control.
empirical bayesian testing effectively estimated the q values for multiple simultaneous hypotheses testing, leading to good false positive control.
a regularized but non-sparse estimation for the glm-based spike response model improved the detection accuracy.
conventional testing procedures have suffered from difficulties in approximating the false discovery rate, particularly when a likelihood-ratio statistic is biased, for example, when with a regularization; our approach based on empirical bayesian testing is a reasonable solution to the difficulty.
our new method’s contribution is prominent, especially when the sample size is relatively small and there are short observation periods like in many in vivo and ex vivo imaging experiments.
a typical example can be seen in functional multi-neuron calcium imaging, where high-temporal resolution restricts the observation length.
in addition, we found that empirical bayesian testing on arbitrary statistics that represent the shapes of spike response functions attained a similar performance to that using the well-established likelihood-ratio statistic.
this finding is important for increasing computational efficiency, because the likelihood-ratio statistics must be calculated after fitting glms for all pairs of null and alternative hypotheses, and hence are computationally heavy.
when applied to a functional multi-neuron calcium imaging dataset from a rat hippocampal ca <dig> region, we found significant functional connections that are possibly mediated by ampa and nmda receptors.
accordingly, our method exhibited reasonably good functional connectivity results even from relatively short observation times and could become a powerful statistical tool in studies of connectomics.
additional file
<dig> /s12868-016-0255-x supplementary document.
abbreviations
glmgeneralized linear model
srfspike response function
arauto-regression
epsc excitatory  postsynaptic current
pcaprincipal component analysis
aicakaike’s information criterion
ciconfidence interval
rocreceiver operating characteristic curve
aucarea under the roc curve
fdpfalse discovery proportion
fdrfalse discovery rate
hhhodgkin–huxley
vbvariational bayesian
ardautomatic relevance determination
roiregion of interest
authors' contributions
so, kn, and yi designed the pre-processing procedure.
so, kn, sa, si, and jy designed the main statistical procedure and wrote the source code.
jy conducted the hodgkin–huxley neuron simulation.
so, si, yi, kn, and jy wrote the manuscript.
all authors read and approved the final manuscript.
competing interests
the authors declare that they have no competing interests.
availability of data and material
the corresponding data and software are available on request from s. oba at oba@i.kyoto-u.ac.jp.
funding
this work was supported by kakenhi grant number  <dig> from mext-japan, crest from japan science and technology agency , platform project for supporting in drug discovery and life science research  from japan agency for medical research and development , program for brain mapping by integrated neurotechnologies for disease studies  also from amed, and project of next-generation core robot and ai technology developmentthe from new energy and industrial technology development organization .

 connectomics, which seeks to identify the connectivity structure between all pairs of neurons not only in local circuits but over the entire brain, is a crucial research direction, because the brain’s functions are believed to emerge within the connectivity structure of its constituent elements, the neurons. toward such a direction, functional connectivity analyses provide a powerful bottom-up tool to investigate the neurophysiology of the relationships of multiple neurons. moreover, a decoding study has quantified information amount that is encoded into functional connectivity of retinal ganglion cells. when we focus on information processing principles of local neuronal circuits, we need a sophisticated approach to unify top-down simulations and bottom-up experimental observations both in terms of functional connectivity and anatomical connectivity. the reliable detection of functional connectivity is thus vital to identify the functions of local and global networks in the brain. similar attempts have focused on the functional connectivity of the higher-order view of the brain: dynamic causal modeling of the brain regions, for example. such research is sometimes called macro-connectomics. although our functional connectivity analysis method can be applied to such macro-connectome problems, in this study we focus on functional connectivity between neurons, one of the elementary levels of the brain’s hierarchy. by definition, a pair of neurons is called functionally connected if the physiological activities of these neurons interact with each other. granger causality provides a theoretical definition to a directed causal relationship that can be estimated from time course observations like spike trains from neuronal circuits. in contrast, structural connectivity implies anatomical observation of synaptic structures between the neurons found by electron microscopy or the imaging of synaptic markers. functional connectivity and structural connectivity have different theoretical implications. structural connectivity does not necessarily imply functional connectivity; since most synapses are silent or inactive in vivo, few active synapses establish functional connectivity. functional connectivity does not necessarily require structural connectivity; since the observation cannot be complete, functional connectivity may be produced by indirect causality bypassed by unobservable neurons and glial cells and mediated by various transmission channels like electronic coupling and metabolic factors. some studies have shown indirect evidence of relationships between functional and structural connections. because one major objective of the physiology of neuronal circuits is to determine their functions, functional connectivity analysis can be more important than its structural counterpart. in order to estimate functional connectivity among many neurons, we require sufficient data of electrical activity of the neurons. recent functional multi-neuron calcium imaging technique enabled us to observe hundreds or thousands of neurons’ activities simultaneously. however, there are following two trade-offs between information amounts that can affect the quality of functional connectivity estimation. the first is a trade-off due to limited observation speed of a scanning microscopy; we can observe the larger number of neurons by the wider field of sight with the larger spatial resolution, leading to the lower temporal resolution. the second is a trade off due to unavoidable photobleaching; higher signal to noise ratio requires stronger light emission, leading to stronger photobleaching of fluorescent dyes that shorten the observation time length. thus, these two trade-offs limit information amount per neuron even when the number of simultaneously recorded neurons is large. moreover, information amount per a pair of neurons is decreasing with growing number of neurons because the number of pairs of neurons is quadratic to the number of neurons. we call this situation an information shortage for functional connectivity estimation. in the current study, we focus on generalized linear models of spike response functions and their extensions, which are found powerful statistical tools to sketch the functional connectivity between multiple neurons. an srf expresses the increase and decrease of the spike generation probability triggered by a pre-synaptic spike input as a function of the lag-time between pre-synaptic input and post-synaptic output spikes. based on the glm framework, srfs for all possible pairs of neurons are estimated simultaneously, and the granger causality test substantially reduces false positives caused by indirect causality that cannot be distinguished if they are estimated for each local pair of neurons. recently, several studies have presented sparse estimation methods that prefer zero connection weights for a large part of all neuron pairs to improve the estimation of glm-based srfs. on the other hand, there have been a few attempts of network structure inference without focusing on fitting glm, such as information theoretic analysis of granger causality and inference of dynamic bayesian networks. in case of information shortage, two major requirements to estimate srfs, model fitting and detection accuracy, are not compatible with each other. estimation of glm, in case of information shortage, owes its performance much to regularization term, in which regularization coefficient is tuned for better model fitting. recent sparse estimations of glm are not exceptions. however, the best-tuned coefficient for model fitting does not generate the best results for detection accuracy; not only that there is no guarantee of compatibility, but also that there is a severe contradiction. discrimination accuracy by the sparse estimation of srfs usually depends on the arbitrary tuning of a hyperparameter that controls the sparseness. the hyperparameter tuning, however, involves severe problems that are not well understood. namely, there is no guarantee that the best hyperparameter to minimize the fitting error in terms of cross-validation leads to the best discrimination accuracy to detect functional connectivity. actually, the connectivity detection performance is not favorable when we tune the sparse model so as to maximize the likelihood. moreover, no good method exists to control false positives with sparse estimation mainly because it is difficult to obtain closed-form null distributions of testing statistics. there have been many studies focusing on model fitting, but few studies focusing on detection accuracy except for the likelihood-ratio testing of granger causality. however, they did not consider the case of information shortage and did not include any regularization factors. the performance of their estimation might be sub-optimal, particularly with information shortage. in this study, we propose a new functional connectivity analysis method focusing on improving the detection performance especially in the case of information shortage. our method consists of two major procedures. in the first, srfs are identified by estimating the parameters of glms with smooth bases to represent srfs and l regularization. the second procedure is a granger causality analysis based on empirical bayesian testing. because our situation for functional connectivity analysis is a typical multiple simultaneous testing, false positive control is definitely important. empirical bayesian testing provides a reliable way to stably control the false positive proportion. our research scope is summarized in table we focus on the case with information shortage in which the number of neurons is large but the observation length is relatively short. because the amount of information available to estimate srf is proportional to the number of spikes per neuron pair, the statistical uncertainty is large in various imaging studies, like calcium imaging. thus, we propose a methodology that is essentially different from best practice in electrophysiological studies. namely, we apply multiple testing to determine the functional connectivity, rather than to estimate srf, for each pair of neurons. the goal is to list a set of candidate connections while controlling false positives, rather than to build a single most likely model that maximizes the fitness or expected likelihood. toward this goal, we apply a glm model with non-sparse regularization in contrast to one with sparse regularization.table conceptual difference between two distinct typical cases of functional neuronal connectivity analyses spike response estimation generalized linear model of spike responses in this section, we introduce a formulation of glm of spike responses according to. suppose that we have a spike train dataset, n={ni:i= …,c;t= …,t}, where ni denotes the number of spikes of the ith neuron in the tth time-bin with a common width over all the bins. we assume the bin width is short enough so that the maximum spike number in a single bin is unity, and thus ni takes either one or zero. we also assume that c neurons may receive a common input signal transmitted through l external channels: e={el∈r:l= …,l}, where el denotes a real-valued signal sent through the lth channel. a stochastic spike response model represents the conditional probability to observe a spike of the ith neuron in the tth time bin pi≡pr=1|n,e,ri), where conditions n, e, and ri denote the past spike history up to time t- of all neurons c= …,c including the ith neuron itself, the external input at time t, and a set of parameters that define the response function of the ith neuron, respectively. according to glm, the conditional probability is a function of total spike response λi, pi=f), where the total spike response is given as a linear summation of all external inputs el and the past spikes from all neurons n, λi≡λi,e,ri)=ri0+∑l=1lrileel+∑c=1c∑s=1mricnc. ri rile, and ric are respectively a background activity level of the ith neuron, a response coefficient to the lth external input, and the spike response with lag-time s to the cth neuron, of the ith neuron. natural number m, called a window, is the maximum time-lag considered in the history. f is a link function from total response λi to probability pi, and we concentrate on a case of logistic link f=1/). a reasonable alternative may be f=1-exp), called the complementary log–log link function. the latter case makes this spike response model equivalent to a poisson spike model at a limit of infinitesimal bin-width. we may consider common external input e that effectively represents the total effect received by the observable neurons from unobserved ones and/or other external inputs. if we can directly observe those external inputs, the estimation of the spike response model is reduced to an auto-regression problem; by regarding glial activities as observable external inputs to a neuronal network, a glm-based spike response model for a neuron-glia system can be estimated. if we cannot observe the external inputs, the estimation is of an ar type with a moving average, in which we need to simultaneously estimate the external inputs and the response functions of individual neurons. estimation and errors we may obtain a small but non-zero estimation of response functions for some pairs of neurons that in fact have no functional connectivity because of the finite amount of available data. such limited data causes statistical uncertainty preventing a clear-cut determination and two types of statistical errors, false negatives and false positives. to deal with such statistical uncertainty in our application of functional connectivity analysis, we have several options in each of the following two steps that constitute functional connectivity analysis: the estimation method of a glm-based spike response model including spike response ric between neurons i and c, and the selection of statistics to test the null hypothesis, where the spike response ric is zero for any time-lag s if the null hypothesis is true, and ric is not zero for some s if the alternative hypothesis is true. our ways of dealing with these two issues will be discussed in the following subsections. application of smooth bases and regularization term this subsection describes a way to estimate a glm-based spike response model. we assume that each spike response function, ric, is represented as a linear combination of a small number of smooth bases: ric=∑k=1kaickbk, where bk and aick denote the kth basis function that is shared by all neurons and the kth basis loading coefficient between neurons c and i, respectively. to reflect the consistently positive character of the facilitatory post-synaptic current, epsc, and the consistent profile of epsc induced by a single pre-synaptic spike, each basis function is given by the following gamma density function: bk=ga. here, ga denotes the probability density function of a gamma distribution whose mean and variance are set at mk=vk=12k the gamma density function as a filter basis was first proposed in. a basis function with a small index k peaks at a small s value, which induces a large but short-delayed epsc or ipsc, and that with a large index is broad with a large mean value, which induces a small but long-lasting epsc or ipsc. number of bases k is arbitrarily determined so that the glm model fits the real spike response well. smoother response function is preferred in the estimation when a smaller number k is set, which can improve the estimation by preventing over-fitting to statistical uncertainty especially in cases with information shortage. from eqs. and, the total spike response of the ith neuron is given by λi=ri0+∑l=1lreilel+∑c=1c∑k=1k∑s=1maickbknc. when estimating the spike response model parameters, the following regularized log-likelihood function is independently maximized with respect to the unknown parameters ri reil, and aick, l= …,l,c= …,c,k= …,k, for each neuron i: li=∑t=1tlti-ηregi, where lti≡nilog)+)log)) is log-likelihood at time t, η is a hyperparameter that controls the strength of regularization, and regi is the regularization term. we may use the following l regularization to encourage non-sparse estimation of the parameter aick: l2:regi≡12∑c=1c∑k=1kaick we alternatively consider the following l or group lasso regularization for facilitating sparse estimation: l1:regi≡∑c=1c∑k=1kaick, gl:regi≡∑c=1c∑k=1kaick the common input term el in our glm-based spike response model is fixed at an estimated value in a preprocess prior to the model estimation. the parameter estimation to maximize the regularized likelihood function is implemented using a dual augmented lagrangian method. estimation of external input when some common external inputs el,l= …,l were supposed to exist but not given, they were estimated using principal component analysis of the observed time course. let us imagine a typical case that there are some synchronized spikes of many observed neurons. in this case, we would detect facilitative connectivity between all neurons that emit the synchronized spikes. however, we can reduce the effect of the synchronized spikes to the estimation of functional connectivity between the observed neurons if we regard the synchronized spikes as an effect of common external input. we can extract the synchronized signal as the first principal component of a set of smoothed spike sequences of the neurons, or multiple signals as some components if there are some distinct sets of neurons that are synchronized. note that we assume these common external inputs as a minimum device to reduce the harmful effect to functional connectivity estimation, rather than to infer any external reality. we first applied a moving average filter of a certain window size to the spike time course ni of all the observed neurons i= …,c to obtain smoothed spike density profiles. next, we calculated l principal components with eigenvalue decomposition of c×c covariance matrix for the smoothed spike density profiles of the c observed neurons and regarded them as the estimated l external inputs. finally, we fit the glm model involving the common external inputs and assessed the fitness using akaike’s information criteria to determine the number of principal components l, where aic is defined as the log-likelihood value minus the number of free parameters of the glm. if there is no external input, the smoothed spike density is based on the internal fluctuation stemming from the constant term of the observed neurons. then there is no significant principal component of the smoothed spike density profiles. if there are principal components of the smoothed spike density profiles, they cannot be explained by the bias of the spontaneous activities of the observed neurons; therefore, they are based on the external inputs. estimation of the external inputs by the above method corresponds to estimation of slow behaviors of the network and is useful for removing factors that are not able to be represented by our glm; it is beneficial to our main objectives of estimating the spike response functions and their significance measure in a stable manner. statistical significance here, we focus on statistical tests for functional connectivity analyses of glm-based spike response models. our statistical tests are based on granger causality, but there are several options when applied to glm-based spike response models. granger causality and a corresponding test statistic granger causality is a criterion to determine the existence of a directed causal relationship between two nodes, from both of which we observed time sequences. by definition, it is said that there is granger causality from nodes a to b if we can significantly better predict the future time sequence of node b using the information of the past and current time sequences of nodes a and b than only using that of node b. in the context of the functional connectivity analysis of our glm-based spike response models, the granger causality from the cth to the ith neurons can be examined by hypothesis testing with null hypothesis h0:ric= for all s= …,m against alternative hypothesis h1:ric≠ for some s= …,m. kim et al. proposed a likelihood-ratio test of all h i,c= …,c by a log-likelihood-ratio statistic: lric=- where l0ic and li are the log-likelihood of the null and alternative hypotheses, respectively. this statistic reflects the loss of model fitness by omitting the functional connectivity from neuron c to neuron i; the alternative likelihood li is first calculated for the post-synaptic neuron i and the null likelihood l0ic is calculated by omitting the contribution of each of the pre-synaptic neurons c= …,c. in their likelihood-ratio test, the alternative log-likelihood is given by eq. without the l regularization, that is, with η= and the null log-likelihood is given by eq. with ric= for all s= …,m. note that the null log-likelihood depends on the pre-synaptic neuron index c, but the alternative log-likelihood does not, because the null likelihood has an additional constraint dependent on the pre-synaptic neuron index, which is not involved in the alternative likelihood. we can calculate the p value of the likelihood-ratio test based on the fact that the test statistic lric obeys a chi-squared distribution with m degree of freedom in a large sample limit. considering the situation of multiple hypothesis testing, q values were also calculated based on the p values for all h i,c= …,c based on the standard procedure. kim et al.’s study followed the standard methodology of statistics; especially for a relatively small amount of data, however, parameter estimation without any regularization may suffer from failure or instability, producing an unreliable calculation of the likelihood-ratio. one possible way to overcome this difficulty is to introduce regularization to the null and alternative likelihood functions, but then the test statistic no longer obeys the asymptotic chi-square distribution. moreover, parameter estimation for the null likelihood must be performed for every pair of pre- and post-synaptic neurons, which is computationally expensive. empirical bayesian testing empirical bayesian testing uses empirical null samples in place of null distributions that are analytically defined like the one described in the previous sub-section. in our particular application, we do not know the null distribution analytically because of the introduction of regularization and kernel functions to the glm-based spike response models; in this case, empirical bayesian testing is a practical choice. in our empirical bayesian testing, a certain number of empirical null samples of test statistics are obtained by the following surrogation method. we assume several surrogate neurons, each of which emits an artificial spike train, by time-shifting a real neuron’s spike train: nc∗=nc for some c∈ …,c. here, we also assume that there is no prominent periodic activity in the network and set the time-shift ts to a number larger than the time-lag m. then, the spike train of a surrogate neuron c∗ becomes independent from that of any real neuron i= …,c, in the time lag. after adding a certain number of surrogate neurons to the set of real neurons, we estimate the glm-based spike response model. in this estimated spike response model, the test statistic between a real neuron a and a surrogate neuron c, which was generated by time-shifting the spike train of a real neuron other than a, can be regarded as a null test statistic, because a and c are independent, with no functional connectivity between them. note that there is a non-zero probability to include some statistic values from the non-independent pairs of neurons in the surrogating process, which can weaken detection power. however, the risk is small because the strength of ordinary response is negligible for a time lag larger than a certain limit and the risk can be smaller by setting large enough time-shift for the surrogate. more importantly, inclusion of false negatives does never violate conservativeness. in multiple simultaneous hypothesis testing like functional connectivity analysis, q value, which is an estimation of false discovery rate, is often used as an alternative significance measure to a p value. when test statistics are available along with empirical null samples for calculating them, q values can be estimated directly from them without calculating the p values. in our empirical bayesian multiple simultaneous hypothesis testing, we used this direct method to obtain q values for all the multiple hypotheses. this procedure is described in the additional file shape-related statistics according to empirical bayesian testing, we can use an arbitrary test statistic. we also examined the following statistics since they are easier to compute compared to likelihood-ratio statistic after estimating the glm-based spike response model, but they still represent the non-zero character of the spike response function. we prepared several statistics to characterize the spike response function between neurons c and i. surface sicsurface≡∑s=1m|ric| peak sicpeak≡maxs|ric| parameter vector norm with a special metric design we evaluate the norm of parameter vector aic· with a special metric design Σkk′. sicmd≡∑k,k′aickaick′Σkk′, where Σkk′ is the -th element of a regularized inverse of the following covariance matrix s of the parameter vectors of null links s0∝∑∈h0aic·taic·. here, aic· is a k-dimensional parameter vector whose k-th element is aick and t denotes a transpose. in the equation above, the summation was taken over every pair of pre- and post-synaptic neurons whose functional connectivity was determined not to exist by the statistical test. thus, the md incorporates the anisotropic nature of the null distribution and evaluates the distance between each case and the typical null cases. delay sicdelay≡argmaxs|ric| the last term, delay, is a bit different, because it is used for describing a functional connection detected by a certain statistical test, rather than for defining a test statistic to detect the functional connectivity. statistical tests for functional connectivity analysis in this study, we compared the following statistical tests in a scenario of functional connectivity analysis based on a glm-based spike response model. they examine a pairwise functional connectivity, that is, whether a direct relationship exists from a pre-synaptic neuron c to a post-synaptic neuron i. chi if there is no regularization term, a chi-square test can be applied to likelihood-ratio statistic lric, because it asymptotically obeys a chi-square distribution with m degree of freedom. then, we can approximately use a chi-square test by setting the regularization hyperparameter η to a small positive value to avoid a divergence in the optimization procedure. a p value is simply obtained by integrating the chi-square distribution. a q value is estimated using the p values of all null hypotheses h based on a standard procedure. dof a p value and hence a q value were obtained by the same procedure as in chi above, except that the degrees of freedom were determined to fit an empirical null distribution of the likelihood-ratio statistic lric. here, the empirical null distribution was produced by the surrogation method described in “shape-related statistics” section. eb empirical bayesian testing was applied to the likelihood-ratio statistic after the model parameters were estimated by the optimization with a larger regularization hyperparameter η. by comparing the results of eb and chi, we can see the contribution of the l regularization and employing smooth bases. eb-arbitrary stat because empirical bayesian testing can be applied to an arbitrary test statistic, the statistics described in “shape-related statistics” section were also examined in its framework. note that we did not apply empirical bayesian testing for estimation results with l and group lasso regularization, because the empirical null distributions of the statistics in “shape-related statistics” section were found to be quite different among neurons, which was not the case with the l regularization. confidence interval and maxz a 95 % wald confidence interval is calculated for each estimated value of glm weight wk, wk± σkk, where σkk is an asymptotic standard deviation that is derived using hessian of the objective function. according to, directed connectivity is detected if any ci corresponding to the connectivity is significantly different from zero. we also consider the z-score of a weight value zk=|wk|/σkk and that of the response function z=∑k=1kzkbk. we integrate them into a maxz statistic, a maximum of z,s= …,m corresponding to the directed connection, so that a criteria, maxz>, is equivalent to the application of ci in. as the straight application of the 95 % ci is no longer a valid multiple testing, we utilized our eb-arbitrary stat framework in order to include the maxz statistic in the following comparison study.  experiment 1: poisson model simulation here, we examined the performance of the functional connectivity analysis methods using a dataset generated by simulating a network of glm-based spike response models. we designed an artificial neuronal network of neurons consisting of three groups of neurons, g g and g each of which consisted of five neurons. this is a recurrent network: g → g → g → g from g to g there are ten excitatory connections among pairs of pre- and post-synaptic neurons. similarly from g to g there are ten excitatory connections among pairs, and from g to g there are ten inhibitory connections among pairs. there is no connection between neurons belonging to the same group. the activity of each neuron was simply determined by a glm-based spike response model, eqs. and, with the complementary log-log link function that simulates a poisson spike model. each neuron also received an external white poisson input of around 50–150 hz independently from the external input to the other neurons. this assumption that the external input was independent between neurons is slightly different from the one assumed by our glm model, in which the external input is shared by the observed neurons. however, as this situation of independent poisson inputs is much simpler, our glm-based spike response model can deal with it by just removing the estimation procedure of the external input. the activities of all neurons were observed at a sampling frequency of 200 hz. with these settings, we obtained from to spikes on average per neuron per 1 s.fig. a simulated neural network that was used in experiment three groups of neurons g g and g construct a recurrent network: g → g → g → g gray scale blocks in the matrix denote a direct functional connectivity between source and destination neurons. b–g performance comparison between glms based on sparse estimation with group lasso and on l regularization. comparison between logistic regression and poisson regression is also performed. the glm parameters were estimated with various settings of hyperparameters and the best hyperparameter value for each case was determined so as to maximize likelihood on a validation dataset. the results were evaluated with sensitivity and false positive proportion calculated by using a true directed link set for the simulation. experiments with short and long datasets are shown in the upper and lower sets of panels, respectively. b, e roc curves are drawn for each case at the best hyperparameter that maximizes likelihood on the validation data. c, f auc score for each setting of regularization hyperparameter. markers denote the values at the best hyperparameter that maximizes likelihood on the validation data. d, g accuracy of model fitting, measured by likelihood on the validation data, for each setting of regularization hyperparameter. markers denote the values at the best hyperparameter that optimizes model fitting. we found that the differences were small between logistic and poisson regression models. roc curves or auc were not necessarily the largest at the best-tuned values of hyperparameters for each of the four cases. sparse estimation with the group lasso was more sensitive to the hyperparameter settings than the l regularization first, we compared results of the sparse and non-sparse estimation methods for glm. we selected spike train data of and samples from the simulated data for training and samples for validation. the sample numbers in both cases are considerably shorter than some typical cases; as a reference, kim et al. used samples to demonstrate their method. when using samples, we found our method achieved 100 % accuracy, although the results are not shown here. we compared two types of regularization, group lasso and l for glm, corresponding to sparse and non-sparse estimation, respectively. we also used l regularization, but did not show its results because the performance with l was comparable to that of group lasso. we also compared two link functions, logistic and poisson; here poisson stands for the complementary log-log link function that simulates poisson spikes. the srfs were estimated using the training data with various hyperparameters that determine the strength of regularization. the estimated srfs were evaluated with two criteria, detection accuracy of the true functional connectivity and fitness to unseen data; detection accuracy is shown as a receiver operating characteristic curve and the area under the curve score, whereas fitness is given by the log-likelihood for the samples for validation. in sparse estimation like with group lasso regularization, some statistics like peak may be exactly zero for some connections. although a naive threshold to detect functional connectivity based on such statistics would be zero, setting a non-zero threshold can in fact produce better functional connectivity estimations. therefore, we evaluated roc and auc by changing the threshold value even for sparse estimation methods. figure  shows the results. in all aspects, logistic and poisson link functions did not show any significant differences whereas the difference between sparse and non-sparse estimation was fairly large. when the hyperparameter was determined to increase the fitness to the validation data, the connectivity detection performance was poor as shown in the roc curves and. as shown in the panels and, the best hyperparameters in terms of the detection accuracy or auc were smaller than those for the largest fitness. the auc could be made larger by arbitrarily setting the hyperparameter to a smaller value. with respect to the fitness, the best model was achieved by group lasso in both sizes of training data, and as shown in and, respectively. with respect to the detection accuracy, the best auc with the l regularization was better than that with group lasso for t= and the best aucs were with both l and group lasso for t= as shown in and, respectively. more importantly, l regularization performed well for a wide range of hyperparameters, while group lasso with non-optimal hyperparameters performed noticeably worse than the best performance. we then compared the results obtained by the three statistical tests, chi dof, and eb, described in “statistical tests for functional connectivity analysis” section. the glm-based spike response model was estimated to maximize the l regularized log-likelihood. when performing the chi and dof, the regularization hyperparameter was set at η=. when performing the eb, we compared two cases eb-lr and eb-regularized-lr, in which the regularization hyperparameter was set to η= and η=, respectively. the value η= was arbitrarily determined to obtain similar results to those of η= but with computational stability. the value η= was also arbitrarily set to a small positive value, because there was no predominant way to set it for enlarging the detection power. we also used smooth bases in the eb-regularized-lr. figures  and show the results. the upper panels show the roc curves with auc scores in the title, in both of which the higher and the larger are the better ones, respectively. in the lower panels, the false discovery proportion is shown for each q value threshold on the horizontal axis; fdp is the proportion of false positives in such links that were detected as functionally connected because the corresponding q values were smaller than the threshold on the horizontal axis. the lower panels evaluate the quality of the q value estimation. because a q value is defined as the fdp estimation, when the fdp is smaller than the q value threshold, the corresponding q value estimation is said to be conservative. from the upper panels of fig.  the roc curves and the auc scores of the proposed method consistently outperformed those of its un-regularized version, and such a benefit of regularization was prominent especially when the observation length t was small. in the lower panels of fig.  the conservativeness of the q value estimation was compared among eb-lr, chi, dof, and the proposed eb-regularized-lr. if a line is close to the diagonal x=y line, the corresponding q value estimation is faithful, and if a line is below the diagonal line, the q value estimation is conservative. we prefer q value estimation that is faithful but conservative; that is, we do not prefer a line that goes over the diagonal line. we found that the empirical bayesian q values were stably good, while a classic chi-square test without dof adjustment almost failed to perform false positive control when the observation length t was small. when the observation length t was large enough, the results depicted by the blue lines were consistent with those reported by. in the lower panels, the fdp values sometimes went over the diagonal line, violating the conservativeness. this violation comes from the estimation variance of q values; a large variance in the q value estimation may arise especially when there is significant correlation between hypotheses, which is the case in the functional connectivity analysis of neuronal networks.fig. results of experiment in which total error and fdp control were compared using several statistical testing methods. upper panels show receiver operating characteristic curves, where horizontal and vertical axes denote the specificity and sensitivity of the statistical test, respectively. red and blue lines show the eb-regularized-lr and eb-lr results, respectively. in the title of each column, observation length t and the area under the roc curve scores of the two methods are shown. lower panels display false positive control, where horizontal and vertical axes denote q value threshold and false discovery proportion when the q value threshold on the horizontal axis was used, respectively. because the q value threshold is regarded as an estimation of the fdp value, when fdp is located under and over the thin black diagonal line of x=y, the corresponding q value estimations are considered conservative and aggressive, respectively. a conservative line is preferable, because by definition, the q value should be a conservative estimation of the fdp value. if the line is close to the diagonal line, the q value estimation is faithful. blue broken, dotted, and solid lines represent the q values estimated by chi2-lr, dof-lr, and eb-lr for the likelihood-ratio statistics without regularization is the likelihood-ratio statistic with larger regularization. the six panels from the leftmost to the rightmost columns correspond to the experimental settings of observation lengths t=,,, and fig. results of experiment in which total error and false positive control were compared among four test statistics: peak, maxz, surface, and md. see caption of fig.  for explanation of horizontal and vertical axes of upper and lower panels fig. results of experiment in which spike trains ware observed from a recurrent network of non-linear conductance-based neurons. notations are identical as in fig. 2 figure  shows the results of statistical tests employing shape-related statistics. we did not find any substantial differences either in total detection accuracy or false positive control among the statistics of the peak, the maxz, the surface, and the md. it is noteworthy that we numerically confirmed the stable false discovery rate control regardless of possibly different null distributions of various test statistics. this result may allow us to select arbitrary test statistics in the proposed empirical bayesian testing framework. we may expect good roc by selecting a relevant test statistic that fits particular shapes of true spike response functions although the rocs were similar between different statistics in the current case partly because we assumed simple spike response function shape such that caused no contradiction between test statistics that we compared here. the roc curve of maxz statistic included 95 % ci criteria as a special case, which is specified as a marker ’o’ on the curve. we found that the 95 % ci provided a moderate balance between sensitivity and specificity in these cases. however, note that the 95 % ci had no guarantee of conservativeness that eb could provide as q value. experiment 2: a network of non-linear conductance-based neurons next, we applied functional connectivity analysis to a model mismatch case. a spike train dataset was obtained by simulating a recurrent network of non-linear conductance-based neurons. each neuron was modeled by single-compartment hodgkin–huxley equations, consisting of sodium, potassium, and leak channels. thus, the neuronal activities did not follow glms. we, however, expect that the proposed false positive control procedure works reasonably well because our framework mainly depends on the empirical null distribution of test statistics built upon a set of surrogate neuron pairs and hence does not depend on the glm’s fitness. the simulated recurrent network resembled that in experiment 1: three layers, each of which had five hh-type neurons. the simulation was conducted using the nest simulator, and the parameters of the hh neurons were set by the default setting of the nest simulator. the results are shown in fig.  in which similar behaviors of the statistical tests can be seen to those in experiment the q value estimation by chi-square tests, chi2-lr and dof-lr, was substantially poor even when observation length t was large. when employing non-linear neurons like the hh-type ones, glm could no longer well represent their non-linear behaviors, and the model of the alternative hypotheses would contain bias, which was not resolved even after very long observation. although the glm-based spike response model had bias, the empirical bayesian testing with regularized model estimation, eb-regularized-lr, showed reasonably good roc curves and auc values and faithful estimation of the q values, based on an empirical null distribution constructed by the surrogation method.fig. results of experiment that considers partial observation with an external input. a external input signals. the artificially designed and an estimated signals are shown as blue and red lines, respectively. time samples out of time samples are shown here. b model selection via aic criterion. vertical and horizontal axes denote aic criterion and number of external inputs that are involved in the glm. aic values are normalized so as to that with no external input is zero. c the estimated spike response functions. five panels correspond to the following categories of neuron pairs; not connected, directly connected, connected via one interneurons, connected via two or more interneurons, and connection from/to surrogate neurons. red line denotes statistically significant directed neuron pairs of q< fig. results of experiment that assessed functional connectivity detection performances at various spike frequency and connectivity weights with fixing time sample number at t= spike frequencies and connectivity weights are controlled by parameters a and w, respectively. rows and columns of the panels denote variations of a and w values, respectively. in each panel, we compared four roc curves; red and blue line color denote estimation methods of glm, l regularization and variational bayes estimation, respectively. real and dotted line styles denote peak and maxz statistics, respectively. the combination of glm with l regularization and peak statistic performed the best for most of the panels fig. analysis of calcium imaging data. a averaged calcium fluorescent image of a rat ca region. b neurons and probable directed connections between the neurons. red and blue arrows denote putative facilitatory and suppressive connections, respectively. statistical properties of the empirical bayesian testing fig. statistical properties of the empirical bayesian testing. a the estimated distribution and surrogate null distribution of the test statistic, peak. the blue curve denotes the estimated total distribution, which is a normalized histogram of the estimated statistic of all possible pairs of observed neurons. the green curve denotes the estimated null distribution, which is a normalized histogram of the estimated statistics of the pairs of the observed and surrogate neurons. b the estimated local false discovery rate. vertical and horizontal axes denote local false discovery rate and the test statistic. each blue point corresponds to a pair of source and destination neurons. c a matrix of the estimated test statistic, peak. rows and columns denote destination and source neurons, where the source neurons include ten surrogate neurons. the color code denotes the value of the test statistic, peak fig. spike response functions of top significant functional connections, arranged from the top-left to bottom-right panels in ascending q value order. in each panel, a pair of source and destination neurons’ code numbers and the corresponding q value are shown. horizontal and vertical axes denote the lag-time s and the spike response function, respectively, where sample corresponds to 10 ms fig. spike response functions of six pairs of significant bi-directional connections. in each panel, both directed spike response functions are shown; for 0<s< each spike response function represents ric, while for -80<s< it represents rci, instead. from the definition of spike response functions ric≈rci, which corresponds to the correlation of neurons i and c fig. spike response function and epsp possibly mediated by ampa and nmda receptors. a a mixture of fast and slow epsps, represented by two gamma distributions peaking at and 50 ms, respectively. according to an electrophysiological study, rat ca pyramidal neurons showed sharp, strong and relatively slow epsps peaking at around 10 ms and lying over 10–100 ms, respectively, which were possibly mediated by ampa and nmda receptors. in the mixture, the amount of the slow epsp was set as small as one third of the blue broken line. b the delayed response of the calcium indicator induced by the membrane voltage was represented by a convolution kernel, an alpha function of the delay time of 50 ms. c by convolving the convolution kernel b to the typical epsp, we have a typical response function obtained by our functional connectivity analysis. contributions in the response function from ampa- and nmda-mediated epsps are displayed as blue lines experiment 3: partial observation case we applied the proposed method to a partial observation case involving external input. we generated an artificial spike train by a simulation of time samples using the same artificial neural network as the ones in experiments and again. in the simulation, we added an artificial external input signal to all neurons, where connection weight from the external input to the neurons were set at random positive values. in order to consider partial observation, we regarded out of neurons were considered as non-observed neurons, and spike response functions between observed neurons, namely those of times pairs of pre- and post-synaptic neurons, were estimated by our procedure. the external input was estimated by applying moving average filter of window size and pca. in fig 5a, we show the estimated external input denoted as red line. we compared aic values of seven models with,..., and external inputs, and found that the model with one external input was the best. the estimated spike response functions are shown in fig 5c, where those function shapes corresponding to the times directed neuron pairs were divided into the following four cases: not connected, directly connected, connected via one interneuron, and connected via two or more interneurons. in addition, we also show spike response functions corresponding to a pair of a surrogate and any observed neurons, that is called a category. the spike response functions were estimated as significant for a large part of direct connection, red lines in case in fig 5c. on the other hand, the estimated spike response functions of indirect connections, cases and, were very weak and difficult to be distinguished from those of cases and. the results were identical when the number of principal component was set at zero, two, and more. experiment 4: spike frequency, connection weight, and detection accuracy how do different spike frequency and connection weight affect the detection accuracy in case of information shortage? does the proposed framework keep its advantage against sparse estimation methods in these cases? we generated an artificial spike train of time samples by a simulation of a simple network involving three neurons, and the network had two directed facilitative connections from to and from to whose connection weights, or peaks of true spike response functions, were and w, respectively. here we compared six values of w in {.... }. all the three neurons had the same background activity level a, that was ri in equation, which controls total spike frequency. we compared six values of a in {,,,,, }. as the network size is small, we repeated the same simulation times with different random seeds and averaged the performances in order to draw stable roc curves. we compared two glm estimation methods with the l regularization with an arbitrary setting of λ= and with the variational bayesian estimation. vb approximates a hierarchical bayesian model with an automatic relevance determination technique that determines sparse solution so as to optimize model fitting. for test statistics, we compared peak, as a representative of the simplest cases, and maxz, as a case that is equivalent to the application of 95 % confidence interval. figure  shows the results of experiment the panels correspond to the 6× variations of simulation setting parameters, w and a. four roc curves in each panel illustrate detection powers of the connection from to of weight w by the four methods, l with peak, l with maxz, vb with peak, and vb with maxz. in total, we found that l prominently outperformed vb in most of the cases. we also found that larger a tend to improve the detection accuracy as we expected. interestingly, change of connection weight w did not affect much to the detection accuracy except for the case showing chance level w= test statistic peak performed a little bit better than maxz indicating non-omnipotency of the confidence interval criteria although we could not conclude general superiority of peak statistic, neither. note that the prominent performance of l against those of vb is consistent to the result of experiment it, however, does not stand for the general inferiority of vb. further discussion with general comparison will be found in “alternative options in smoothing and regularization” section. analysis of calcium imaging data we applied functional connectivity analysis with glm-based spike response models to a real calcium imaging dataset. a high-throughput calcium imaging experiment was performed on a cultured slice of the ca region of rat hippocampus. using the nipkow disc, the sampling frequency was as high as 100 hz, and the number of time samples was t=; that is, the observation time was 10 min. during preprocessing, regions of interest, each of which showed its individual activity, and the corresponding time-courses of the fluorescence level were obtained. a visual inspection found that each roi mostly corresponded to a neuron. then, for each fluorescence time-course, we detected the peak timings by applying a fixed template that represented a calcium profile associated with a single spike. a single detected peak may not correspond to a single real neuronal action potential; actually, there were some cases in which only one spike was detected for burst-like calcium activity. we simply regarded such burst activity as a single spike because their number was small. more detail on the above preprocessing is shown in the additional file we selected neurons out of the rois that showed frequent spike activities and found that the highest spike frequency was times larger than the lowest one within them. after this preprocessing, we evaluated the functional connectivity between the neurons and their statistical significance using our method empirical bayesian testing after the regularized parameter estimation of the glm-based spike response model. for each neuron pair, the spike response function of at most an samples time-lag was estimated. we prepared ten surrogate neurons and assumed two channels of external inputs, which were estimated by applying pca analysis to all of the neurons’ activities. we found that the false positive control based on shape-related statistics worked as well as the likelihood test statistic, so the q value threshold was determined using the peak statistic. figure  shows the result of the functional connectivity analysis, in which probable directed connections are shown by arrows. the number of connections was determined so that the estimated false discovery rate was. in this figure, some neuron pairs connected by red arrows and being located close to one another may be single neurons segmented into two parts, such as soma and axon hillocks. even if such pairs are removed, the detected functional connectivity still tends to connect two neurons near each other. there were five neuron pairs connected by blue arrows being located distant to each other. statistical properties of the proposed empirical bayesian testing for this dataset are shown in fig.  the probability distribution of a test statistic is a mixture of unknown null and alternative distributions, which would correspond to the left mass and the right tail of the density function shape in panel, respectively. the surrogate distribution of the test statistic fits well to the left mass part of the empirical distribution, indicating that the surrogate distribution well simulated the unknown null distribution. the local false discovery rate shows that the estimated total proportion of null connections was π0= and that the local false discovery rate was smaller than when the statistic was larger than. the significant pairs of neurons are shown as white dots in the source/destination matrix, so that their degrees of significance are seen by comparing to the background noise level in the surrogate area of the matrix. the spike response functions of significant functional connectivity are shown in fig.  in ascending order of q values. many response functions of the top significant functional connectivity had a positive peak at lag-times of around samples, and some response functions included a tail around the 50th sample, which can be consistently considered as epsps possibly mediated by ampa and nmda receptors. moreover, obvious continuity of response functions between pre-post connections, ric, and their reverse ones, rci, was found for many pairs of neurons,, in fig.  this continuity occurs partly because it reflects the correlation of neurons i and c. such correlation might have been overestimated due to jitter in spike timing detection; in some cases, neuron i emits a spike slightly earlier than that of neuron c, and in other cases, neuron c’s spike can be earlier, because of the jitter in the binning. discussion information shortage and detection performance to detect functional connectivity, we need sufficient events of delayed spike pairs at the corresponding pair of neurons. moreover, the number of events should exceed the background level of spontaneous activity of the post-synaptic neuron. thus, various factors can reduce detection performance via information shortage; shorter observation time length, lower total spike frequency, weaker connection weight, an unknown external input to post-synaptic neuron, high background activity of post-synaptic neuron, and so on. any of these factors unavoidably exists in real situations and consideration of information shortage is always vital for any case even when the observation time length is fairly long. this point of view motivated the current study. and, a part of them were demonstrated in experiments and in this paper. alternative options in smoothing and regularization recently, two distinct approaches have been proposed for functional connectivity analysis from the neuronal spike train data. in the first approach, hierarchical bayesian methods were applied to glm-based models to obtain sparse solutions. a prior in the hierarchical bayesian setting allows the statistical estimation of glm-based models to prefer sparsely connected networks, such as reflecting the sparse structure of neuronal networks. the sparseness also reduces the complexity of each spike response function. the second approach was based on granger causality with a likelihood-ratio test. functional connectivity analysis based on granger causality can be performed for any regression model; the authors used the simplest ar model without any smoothing kernel or regularization to enable the null distribution of the likelihood-ratio test to be analytically calculated. these two approaches are based on different statistical approaches but share the same objective of functional connectivity analysis; combining these approaches has options, such as, how to introduce regularization and how to define null and alternative hypotheses in the likelihood-ratio test. in this study, we evaluated several reasonable options in terms of two criteria in statistical testing: appropriate false positive control, which is based on good estimation of the false discovery rate in multiple simultaneous testing, and high detection accuracy which is determined by increasing the number of true positives while fixing the level of false positives. further comparison to chen et al. may better illustrate the difference between the two approaches. their work proposed two distinct ideas, automatic-relevance-determination with variational bayes to estimate connection weights of glm and a use of 95 % confidence interval criteria to determine the connection. their ard with vb technique estimates sparse model as like the l regularization does with optimizing sparseness hyperparameter, and thus the result tends to be as good as the best among the l of various regularization parameters. the use of the 95 % ci may resemble our work in the consideration of false positive risk control. however, note that both of these ideas were designed based on the first approach for the purpose of fitting a generative model to the observed data rather than to detect functional connectivity. the best model, measured by likelihood, is not necessarily the best model, measured by roc curve. furthermore, the gap between these two criteria tends to be larger at information shortage. we showed this contradiction in fig  of experiment and more prominently in experiment our approach remarkably outperformed the vb if the performance is measured by roc. consequently, we recommend choosing appropriate approach with determining purposes in advance to the analysis. model mismatch bias model mismatch bias is an important issue that needs to be addressed here. any bayesian methodologies using prior probability to improve estimation performance always, at least implicitly, presume that the parametric model includes the truth. although such a presumption can actually improve the estimation performance in ideal situations without any model mismatch, most models include a mismatch to the reality to a greater or lesser extent. there is no guarantee that the glm represents actual neural behaviors well enough even with the best setting of the model parameters. however, our framework is based on statistical testing that allows model mismatch to some extent. namely, the alternative model does not necessarily include the truth because the functional connectivity is determined by rejection of the corresponding null hypothesis. better and worse alternative models lead to decrease and increase in false negative rate, but do not affect false positive control. the false positive is well controlled if the surrogate distribution well simulates the null process, that is, spike response between functionally unconnected neurons in our case. experiment showed a reasonable performance in a typical model mismatch case where the data were generated by the hodgkin–huxley equations and analysis was done based on the glm model. the analysis of calcium imaging data suggested that the surrogate distribution well simulated the null part of the distribution of the peak statistic. we thus conclude that our method can perform conservative statistical testing in many practical cases including that of model mismatch. computational cost computational cost for the glm parameter estimation is proportional to the observation length and to the square of the number of neurons. when we add surrogate neurons for the empirical bayesian testing, it increases the effective number of possible connections and causes a corresponding increase in computational cost. note that the number of surrogate neurons cs may not be so large when the real neuron number c is large because the variance in the statistical test depends on the number of null samples of their inter-connections, csc. we set cs so that the number of null samples was around csc= in our experiments. partial observation partial observation is an important practical issue because it is difficult for many physiological techniques to cover all neurons recruited by the target system; functional connectivity analyzes should cope with this issue. in experiment we demonstrated that spike response functions of indirectly connected pairs of neurons were far weaker than those of directly connected pairs. we conducted similar experiments with various other settings and found similar results, where the variants included different observation length,,, and, different sets of observed and non-observed neurons, different moving average windows in pca analysis, different strengths of external inputs, and different network structures. we consistently observed the following tendencies:when a strong external input exists, it is detected by the proposed procedure with pca and aic. the proposed procedure with pca with aic hardly ever detected a total activity of non-observed neurons as a distinct external input. the effect of non-observed neurons was indistinguishable to a background activity level of each observed neuron. spike response was very weak and not significant for indirect connections. spike response of direct connection was stronger than that of indirect connections but weaker than that in fully observed cases. in analysis of partial observation case with the glm model, the total effect of non-observed neurons is roughly summed up into the scalar value of background activity level, which causes severe loss of information compared to the fully observed case. this information loss has probably caused the higher estimation variance of spike response and the lower sensitivity to detect direct and indirect functional connectivities. in total, experiment demonstrated an additional reason to consider information shortage that might be present in case of long observation length and to need the proposed procedure with empirical bayesian test with surrogate neurons to control false positives. further analysis would be needed as future studies to clarify the effect of non-observed neurons to the estimation variance although it is out of scope of the current study. biological suggestions when we applied the proposed method to the rat ca calcium imaging dataset, we extracted directed connections. many of their response functions exhibited sharp peaks at around samples, while some were also accompanied by dull and weak tails lying over 10– samples. the delayed activation corresponding to the weak tails may reflect sequential spike chains across multiple neurons, which prevail in hippocampal networks and often last for more than 100 ms. more physiologically, the sharp peaks and weak tails may reflect epsp of rat hippocampal pyramidal neurons revealed by an electrophysiological study; reportedly, epsp of rat ca pyramidal neurons showed sharp strong peaks at around 10 ms and relatively weak tails of 10–100 ms delay, which were putatively mediated by ampa and nmda receptors. in fact, when we convolved an alpha function with a delay time of 50 ms, which represents a chelating property of our calcium indicator, to the mixture of the fast ampa-mediated epsp and the late nmda-mediated epsp, we could reproduce a typical shape of our response functions. changes in the mixture of ampa- and nmda-mediated epsps would produce various kinds of response functions. therefore, the variety in the response functions found by our data-driven functional connectivity analysis may have reflected to some extent the difference in the receptor distributions between the pyramidal neurons in the rat ca circuit. there were two prominent negative responses, 16→ and 18→ that also included positive peaks. the positive peak of 18→ was weak and could not be statistically significant, while that of 16→ was still prominent. seeing the response function of the counter direction 46→ the positive peak in 16→ was considered to be an artifact coming from the significant positive peak in 46→. there were no other bi-polar responses that were significantly strong enough.  in this study, we presented a new combination of glm-based spike response modeling and empirical bayesian testing to perform functional connectivity analysis between neurons. even when the observation period was relatively short, our method showed reasonably good detection accuracy while keeping good false positive control. empirical bayesian testing effectively estimated the q values for multiple simultaneous hypotheses testing, leading to good false positive control. a regularized but non-sparse estimation for the glm-based spike response model improved the detection accuracy. conventional testing procedures have suffered from difficulties in approximating the false discovery rate, particularly when a likelihood-ratio statistic is biased, for example, when with a regularization; our approach based on empirical bayesian testing is a reasonable solution to the difficulty. our new method’s contribution is prominent, especially when the sample size is relatively small and there are short observation periods like in many in vivo and ex vivo imaging experiments. a typical example can be seen in functional multi-neuron calcium imaging, where high-temporal resolution restricts the observation length. in addition, we found that empirical bayesian testing on arbitrary statistics that represent the shapes of spike response functions attained a similar performance to that using the well-established likelihood-ratio statistic. this finding is important for increasing computational efficiency, because the likelihood-ratio statistics must be calculated after fitting glms for all pairs of null and alternative hypotheses, and hence are computationally heavy. when applied to a functional multi-neuron calcium imaging dataset from a rat hippocampal ca region, we found significant functional connections that are possibly mediated by ampa and nmda receptors. accordingly, our method exhibited reasonably good functional connectivity results even from relatively short observation times and could become a powerful statistical tool in studies of connectomics. additional file /s12868-016-0255-x supplementary document. abbreviations glmgeneralized linear model srfspike response function arauto-regression epsc excitatory postsynaptic current pcaprincipal component analysis aicakaike’s information criterion ciconfidence interval rocreceiver operating characteristic curve aucarea under the roc curve fdpfalse discovery proportion fdrfalse discovery rate hhhodgkin–huxley vbvariational bayesian ardautomatic relevance determination roiregion of interest authors' contributions so, kn, and yi designed the pre-processing procedure. so, kn, sa, si, and jy designed the main statistical procedure and wrote the source code. jy conducted the hodgkin–huxley neuron simulation. so, si, yi, kn, and jy wrote the manuscript. all authors read and approved the final manuscript. competing interests the authors declare that they have no competing interests. availability of data and material the corresponding data and software are available on request from s. oba at oba@i.kyoto-u.ac.jp. funding this work was supported by kakenhi grant number from mext-japan, crest from japan science and technology agency, platform project for supporting in drug discovery and life science research from japan agency for medical research and development, program for brain mapping by integrated neurotechnologies for disease studies also from amed, and project of next-generation core robot and ai technology developmentthe from new energy and industrial technology development organization.
we also showed that empirical bayesian testing on arbitrary statistics in place of likelihood-ratio statistics reduce the computational cost without decreasing the detection performance. functional connectivity analyses of multiple neurons provide a powerful bottom-up approach to reveal functions of local neuronal circuits by using simultaneous recording of neuronal activity. in this study, we propose a new combination of glm and empirical bayesian testing for the estimation of spike response functions that enables both conservative false discovery control and powerful functional connectivity detection. our method achieved high detection performance of functional connectivity with conservative estimation of false discovery rate and q values in case of information shortage due to short observation time. functional connectivity of neuronsspike statisticsgeneralized linear modelempirical bayesian testingissue-copyright-statement© the author 2016 when our proposed method was applied to a functional multi-neuron calcium imaging dataset from the rat hippocampal region, we found significant functional connections that are possibly mediated by ampa and nmda receptors. we compared our proposed method’s performance with those of sparse estimation of glm and classical granger causality testing. a statistical methodology, generalized linear modeling of the spike response function, is one of the most promising methodologies to reduce false link discoveries arising from pseudo-correlation based on common inputs. the proposed empirical bayesian testing framework with glm is promising especially when the amount of information per a neuron pair is small because of growing size of observed network. although recent advancement of fluorescent imaging techniques has increased the number of simultaneously recoded neurons up to the hundreds or thousands, the amount of information per pair of neurons has not correspondingly increased, partly because of the instruments’ limitations, and partly because the number of neuron pairs increase in a quadratic manner.
 functional connectivity analyses of multiple neurons provide a powerful bottom-up approach to reveal functions of local neuronal circuits by using simultaneous recording of neuronal activity. a statistical methodology, generalized linear modeling of the spike response function, is one of the most promising methodologies to reduce false link discoveries arising from pseudo-correlation based on common inputs. although recent advancement of fluorescent imaging techniques has increased the number of simultaneously recoded neurons up to the hundreds or thousands, the amount of information per pair of neurons has not correspondingly increased, partly because of the instruments’ limitations, and partly because the number of neuron pairs increase in a quadratic manner. consequently, the estimation of glm suffers from large statistical uncertainty caused by the shortage in effective information.  in this study, we propose a new combination of glm and empirical bayesian testing for the estimation of spike response functions that enables both conservative false discovery control and powerful functional connectivity detection. we compared our proposed method’s performance with those of sparse estimation of glm and classical granger causality testing. our method achieved high detection performance of functional connectivity with conservative estimation of false discovery rate and q values in case of information shortage due to short observation time. we also showed that empirical bayesian testing on arbitrary statistics in place of likelihood-ratio statistics reduce the computational cost without decreasing the detection performance. when our proposed method was applied to a functional multi-neuron calcium imaging dataset from the rat hippocampal region, we found significant functional connections that are possibly mediated by ampa and nmda receptors.  the proposed empirical bayesian testing framework with glm is promising especially when the amount of information per a neuron pair is small because of growing size of observed network. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords functional connectivity of neuronsspike statisticsgeneralized linear modelempirical bayesian testingissue-copyright-statement© the author 2016 