Research Article  Open Access
Identification of Functionally Interconnected Neurons Using Factor Analysis
Abstract
The advances in electrophysiological methods have allowed registering the joint activity of single neurons. Thus, studies on functional dynamics of complexvalued neural networks and its information processing mechanism have been conducted. Particularly, the methods for identifying neuronal interconnections are in increasing demand in the area of neurosciences. Here, we proposed a factor analysis to identify functional interconnections among neurons via spike trains. This method was evaluated using simulations of neural discharges from different interconnections schemes. The results have revealed that the proposed method not only allows detecting neural interconnections but will also allow detecting the presence of presynaptic neurons without the need of the recording of them.
1. Introduction
The microelectrodes array technology is a standard tool in the neuroscience field. This technology has allowed observing the activity of neuronal populations, establishing specific correlations between them, and revealing different strategies for processing of sensory information [1–3]. The temporal correlations among spike trains have been associated with the neural coding, stimuli discrimination, and process related to attention [4]. Furthermore, it has been suggested that such correlations are due to the anatomical and functional interconnection between neural networks of the underlying tissue [5]. Thus, the spatio/temporal knowledge of neural interconnections is currently of great interest in the study of the sensory neural code.
Processing techniques, such as crosscorrelation and Granger causality, are often used to establish connections between neurons [6, 7], while other multivariate methods, such as generalized linear model (GLM), directed transfer function (DTF), or partial directed coherence (PDC), allow identifying interconnections in neuronal populations [8–10]. All these methods require that the functional activities of neurons, or neuronal groups, are electrophysiologically registered. However, it is often necessary to identify functional interconnections between neurons whose activity is known with others whose activity is unknown, and, for this, very few processing techniques can be used. It is known that the neural facilitation/inhibition of retinal ganglion responses changes when these are stimulated outside of their receptive fields [11, 12] and that this is due to interconnections between ganglion cells with presynaptic neurons [13]. This situation cannot be studied by conventional methods, since the simultaneous recording of ganglion cells and presynaptic neurons cannot be performed. Thus, a method that reveals such interconnections would be of great interest in the neurosciences area.
Here, we propose a multivariate technique (factor analysis) to detect and to establish functional interconnections between neurons. The factor analysis is a wellknown statistical analysis technique; however, it was not used to detect interconnections in the sensory systems based on response of single neurons (spikes). The robustness of the proposed method was assessed through computational simulations. Our results show that factor analysis allows identifying neuronal interconnections among neurons with known response (i.e., ganglion cells), in addition to those with unknown response (i.e., presynaptic neurons).
2. Materials and Methods
2.1. Generation of Synapses Interconnection
Synaptic interactions and interconnections schemes were modeled in two ways: the first using a model of correlated currents and the second form is using neural networks.
The model of correlated currents consists in generating fluctuations of presynaptic currents using [14]where is the temporal average of the current. The second term represents fluctuations along time and it is composed of the weighted sum of two factors: is the current of th neuron, is the current of presynaptic neuron , () is the presynaptic correlation values of and currents, and is the variance of input current. These currents have a distribution of white noise, and in all cases these have a temporal length of 5 seconds.
The simulated neuronal interconnections are shown in Figure 1(a). The presynaptic current has connections with , , and neurons, and it is weighted by , , and , correlation values, respectively. Likewise, current has connections with and neurons.
(a)
(b)
(c)
(d)
(e)
The correlated current values, , were inputs of Leaky IntegrateandFire (LIF) Neuron Model which was used as spike generator. The parameters used in the model were as follows: membrane time constant (tau) tau = 10 ms, membrane potential threshold = −55 mV, resting potential = −70 mV, and reset value = −75 mV. Thus, the spike trains from , , , , and neurons were obtained. Then, the spike count within a window of 50 ms was determined.
Because the model of correlated current has very strong modeling assumptions about the structure of correlations and generation of postsynaptic currents, it was decided additionally to use a model based on neural networks.
Thus, the second method used to model the neural interconnections consists in using feedforward network [15]. Neural baseline activity is given by uncorrelated homogeneous Poisson processes. For different simulations, the rate parameter was varied (, 1/25, and 1/50) corresponding to level of spontaneous spiking activity. The synaptic weights () took different values according to different situations analyzed. In the simulations with the interconnections of Figures 1(b), 1(c), 1(d), and 1(e), the synaptic weights took deterministic values, while in the scheme of Figure 2(a) the weights took values with uniform distribution (parameters and ). The nonlinear combination of all inputs was deterministically performed through logistic sigmoid function (activation function): where and [16].
(a)
(b)
2.2. Factor Analysis
The factor analysis (FA) is based on the following model [17]:where is a vector of which contain the spike count values of all neurons, is a mean vector (), is a vector of which contain the unobservable factors—spike count values of presynaptic neurons— is a loading matrix of , and is a vector of with unobservable perturbations, which has a multivariate normal distribution of order , with the components of mean vector being equal to zero, and the variance matrix is equal to identity. Thus, vector has a multivariate normal distribution, .
The amount of factors is related to the amount of presynaptic neurons; that is, if two factors are considered in the model, then two presynaptic neurons are considered in the analysis. will be related to correlation values in the interconnection model proposed, and it will have five rows (one for each recorded neuron) and two columns (one for each presynaptic neuron). Finally, the maximum likelihood method is used to estimate the loading matrix . This method consists in finding and values which maximize the following function:where is the amount of time intervals where the spike count was done, is the sample covariance matrix, and indicates the trace of matrix. The parameters that maximize the convex function are obtained through their corresponding partial derivatives (Statistical toolbox, MATLAB).
2.3. Metrics
The matrix 2norm, , of the difference between calculated loading matrix and optimum loading matrix, was proposed to quantify the results obtained (see (4)). belongs to the interval [0, inf) and can be interpreted as the distance, or discrepancy, between the loading matrix and optimum loading matrix. = 0 when both matrices are identical. The optimum loading matrix, , reflects the preestablished interconnections; in the case of Figure 1(a) it is as follows:
2.4. Identification of Presynaptic Interconnections
The synaptic interconnections were determined by comparing the elements of loading matrix. These comparisons were performed in two steps: by rows and then by columns.
First, the position of the element with greatest value is identified through analysis of each of rows of the loading matrix. Thus, the presynaptic neuron which is connected to recorded neuron is determined. The proposed method will not be able establish any connection, if the value of any element of a row is not significantly higher compared to the others.
The previous procedure is applied to all recorded neurons. Then, each of the columns of the loading matrix is similarly analyzed. If there is one or more elements whose values are significantly higher compared to the remaining elements (in the same column), then, such positions (recorded neurons) are connected to the same presynaptic neuron.
Considering the above, a neural interconnection is determined if the following conditions are met:(i)The position of the element with the highest loading value of each row coincides with the elements position of “1” value of the corresponding rows of matrix.(ii)The loading values must be significantly higher compared to the remaining elements of the same column.Two values are considered significantly different if their difference exceeds a preestablished threshold. This threshold was established empirically in 0.03 for proposed case studies. For future implementations, it should be taken into account that a high threshold will allow establishing stronger connections, while a low threshold will allow establishing weaker connections. On the last case, the result could contain a greater number of false positives.
3. Results
3.1. Influence of Weightings
The influence of neuronal correlation values on the identification of interconnections was analyzed through simulations in which value varies while and remain constant ( and ). Then, the identification of interconnections was quantitatively assessed through .
Figure 3(a) shows values for interconnections scheme of Figure 1(a) and considering . It is possible to note that values decrease as , , , and values increase. This trend is most noticeable when , Figure 3(b).
(a)
(b)
3.2. Influence of the Samples Number
values for the neuronal interconnections of Figure 1(a), versus the samples number (), is shown in Figure 4. It is observed that values decrease and converge to a value when the samples number is increased. The convergence values are 1.1, 0.8, and 0.6 for values equaling 0.4, 0.6, and 0.8, respectively.
(a)
(b)
(c)
3.3. Identification of Interconnections from Different Loading Matrices
Different interconnections schemes were proposed for assessing the factor analysis procedure (Figures 1(b)–1(e)). Spike trains from five neurons were recorded in all cases. The average loading matrices for interconnection scheme of Figure 1(b) are shown in Tables 1(a) and 1(b). It is important to clarify that the weighting or correlation values should be varied independently to make a more intensive analysis of the proposed method. However, as a first approximation and to simplify the results, the case where all weights are equal to each other has been considered. For the two models of interconnections (Tables 1(a) and 1(b)) it is observed that the values of the first column are higher than those of second column for different weighting values indicating the existence of a presynaptic neuron. Likewise, the loading matrix reveals that is the presynaptic neuron (highest values) when comparing the elements of the first columns. This difference between the values of the first column is greater in the model correlated currents (Table 1(a)).
(a)  
 
(b)  
 
(c)  

Similarly, Table 2 shows the average loading matrices obtained for the interconnection scheme of Figure 1(c). Spike trains from five neurons were recorded for all the experiments, while the values of neuron are the highest when comparing the elements of the first columns. This particularity indicates that neuron would not be connected to other neurons and that it would fire independently (italic font values) and also that is presynaptic to the other neurons (bold font values).
(a)  
 
(b)  
 
(c)  

The average loading matrices obtained for the interconnection scheme of Figure 1(d) are shown in Table 3. It is observed that the values of and neurons (first column of matrices) are similar to each other and higher values than those belonging to other neurons. The difference of values between the first and second column increases with weighting. For this situation, it is not possible to determine the direction of the neural connection because of the similarity between the values of and neurons. Thus, the loading matrix only would indicate the existence of a common presynaptic neuron.
(a)  
 
(b)  
 
(c)  

In Tables 4(a) and 4(b) (loading matrix of Figure 1(e)) it is observed that the values of both columns are similar for all neurons. This indicates that the method cannot identify a pattern of neuronal interconnection which agrees with the scheme of Figure 1(e).
(a)  
 
(b)  

3.4. Influence of the Number of Presynaptic Neurons
The influence of presynaptic neurons number in the calculation of the “loading matrix” was analyzed. For this, neural network model was used, and scheme used interconnections shown in Figure 2(a). In this case only a presynaptic neuron connected to neurons 1 and 2 with value, while values for the other connections between presynaptic neurons and neurons registration set at random and with different values of . In Table 5 the theoretical loading matrix for the interconnection scheme of Figure 2(a) is observed. In Tables 6(a), 6(c), and 6(e), it is observed that by increasing the number of presynaptic neurons the difference between the maximum and minimum values tends to decrease but still identifying interconnections is correct.

(a)  
 
(b)  
 
(c)  
 
(d)  
 
(e)  
 
(f)  

In Tables 6(d), 6(b), and 6(f), it is observed that by increasing the number of presynaptic neurons values in the first column of the loading matrix take similar values to each other. This causes the correct identification of the presynaptic neuron interconnections with greater weight () being limited by two factors, first by the number of presynaptic neurons and second by the value of lambda. The correct identifications can be made until the number of presynaptic neurons is equal to , 65, and 80 for values 1/10, 1/25, and 1/50, respectively.
3.5. Influence of Correlation between Trains Spikes Presynaptic Neurons
The simulations above meet the condition that the discharges of presynaptic neurons are independent. Now we will study the robustness of FA method when trains spikes of presynaptic neurons are correlated.
The neuronal interconnections for scheme of Figure 2(b) were established by using the correlated current model. In this scheme, the discharge of and presynaptic neurons is correlated by the discharges of neuron. The loading matrices for this situation were obtained applying the FA method to 1, 2, 3, 4, and 5 recorded neurons (Table 7); in Table 8 the theoretical loading matrix for the interconnection scheme is observed. It is observed that, for low correlation values (), FA can correctly identify the structure of neuronal interconnections. For greater correlation values () FA presents problems to correctly detect the structure of interconnections.


4. Discussion and Conclusion
The FA is a multivariate statistical method that has traditionally been applied in the psychology areas and recently in neuroscience areas. Thus, for example, specific functional interconnections between neuronal groups have been established through the factor analysis [18]. Yu et al. propose using an extension of FA (Gaussianprocess factor analysis, GPFA) to detect and model modulatory or underlying neural processes that modify the response of neuronal populations over time [19, 20].
The correct identification of connections between neurons and/or neuronal groups is a problem of interest in neuroscience. In this aspect many techniques to identify functional interconnections have been proposed, but most are limited to the interconnection of neuronal groups with a large number of neurons [8–10, 21].
New advances in electrophysiological methods have allowed registering the joint activity of single neurons, so that a more specific functional analysis could be conducted. Thus, the methods for identifying neuronal interconnections via spike trains are a growing demand in the area of neurosciences. Thus, for example, the influence of the activity of interneurons, or presynaptic neurons (no recorded), on the activity recorded from other neurons, could be of great interest. Echtermeyer et al. have proposed a technique capable of detecting the presence of interneurons whose activity was unknown but interconnected with other neurons whose activity was known. However, it was not capable of detecting interconnections between neurons whose activities are known with presynaptic neurons whose activities are unknown [22].
In this study we have proposed a factor analysis to identify functional interconnections among neurons by using spike trains. Factor analysis is a statistical technique, and to be applied it is necessary to verify the fulfillment of its hypotheses. In this aspect the robustness of the technique applied to neural responses was validated by comparing the results obtained with ideal results, particularly for this study comparing interconnections obtained with FA and imposed by the model.
In order that the FA can be used to identify neural interconnections the following assumptions must be met:(i)Trains presynaptic neurons spikes should be independent or have a low correlation between them.(ii)Recorded spikes trains should not have temporal correlation.It is not possible to verify a priori the independence between spike trains from presynaptic neurons by using real recordings. Therefore the assumptions necessary to validate the results obtained with the proposed method must be corroborated through the anatomical/functional knowledge of the analyzed system. For example, it is known that the amacrine cells of the retina, which are presynaptic of the ganglion cells, fire independently. Thus, the proposed method could be used in retina.
Our results, based on computational simulations, have revealed that the proposed method not only allows detecting neural interconnections but will also allow detecting the presence of presynaptic neurons without the need of the recording of them.
The neuronal interconnection schemes proposed in this study were chosen because of their similarity to those commonly found in the nervous system. Thus, for example, the neuronal interconnections of Figure 1(a) are biologically plausible in the human retina and are given by ganglion and amacrine cells [13], whereas those proposed in Figures 1(b), 1(c) and 1(d) could be found in cortical areas [23].
Although the results found are based on computational simulations, FA could be applied to real neural recordings. For this, neuronal recordings must meet specific conditions (listed above) such as those related to temporal correlation.
In the simulations we have used two methods to generate trains of spikes and interconnections. The results have shown that the efficiency of FA in identifying interconnections does not depend on the method used in generating trains of spikes but rather depends on the probabilistic structure (whether or not there is correlation). In addition the results obtained revealed that there are two factors that influence the efficiency of the method, first the number of presynaptic neurons and secondly the weight of synapses.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work has been supported by grants from Agencia Nacional de Promoción Científica y Tecnológica (ANPCYT), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), and Consejo de Investigaciones de la Universidad Nacional de Tucumán (CIUNT), as well as by Institutional funds from Instituto Superior de Investigaciones Biológicas (INSIBIO).
References
 M. Stopfer, S. Bhagavan, B. H. Smith, and G. Laurent, “Impaired odour discrimination on desynchronization of odourencoding neural assemblies,” Nature, vol. 390, no. 6655, pp. 70–74, 1997. View at: Publisher Site  Google Scholar
 J.M. Alonso, W. M. Usrey, and R. C. Reid, “Precisely correlated firing in cells of the lateral geniculate nucleus,” Nature, vol. 383, no. 6603, pp. 815–819, 1996. View at: Publisher Site  Google Scholar
 J. Shlens, G. D. Field, J. L. Gauthier et al., “The structure of multineuron firing patterns in primate retina,” Journal of Neuroscience, vol. 26, no. 32, pp. 8254–8266, 2006. View at: Publisher Site  Google Scholar
 R. Christopher Decharms and M. M. Merzenich, “Primary cortical representation of sounds by the coordination of actionpotential timing,” Nature, vol. 381, no. 6583, pp. 610–613, 1996. View at: Publisher Site  Google Scholar
 C. Diekman, K. Dasgupta, V. Nair, and K. P. Unnikrishnan, “Discovering functional neuronal connectivity from serial patterns in spike train data,” Neural Computation, vol. 26, no. 7, pp. 1263–1297, 2014. View at: Publisher Site  Google Scholar
 W. Hesse, E. Möller, M. Arnold, and B. Schack, “The use of timevariant EEG Granger causality for inspecting directed interdependencies of neural assemblies,” Journal of Neuroscience Methods, vol. 124, no. 1, pp. 27–44, 2003. View at: Publisher Site  Google Scholar
 D. H. Perkel, G. L. Gerstein, and G. P. Moore, “Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains,” Biophysical Journal, vol. 7, no. 4, pp. 419–440, 1967. View at: Publisher Site  Google Scholar
 M. J. Kaminski and K. J. Blinowska, “A new method of the description of the information flow in the brain structures,” Biological Cybernetics, vol. 65, no. 3, pp. 203–210, 1991. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. W. Pillow, J. Shlens, L. Paninski et al., “Spatiotemporal correlations and visual signalling in a complete neuronal population,” Nature, vol. 454, no. 7207, pp. 995–999, 2008. View at: Publisher Site  Google Scholar
 D. Y. Takahashi, L. A. Baccal, and K. Sameshima, “Connectivity inference between neural structures via partial directed coherence,” Journal of Applied Statistics, vol. 34, no. 910, pp. 1259–1273, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 B. P. Ölveczky, S. A. Baccus, and M. Meister, “Segregation of object and background motion in the retina,” Nature, vol. 423, no. 6938, pp. 401–408, 2003. View at: Publisher Site  Google Scholar
 J. T. McIlwain, “Receptive fields of optic tract axons and lateral geniculate cells: peripheral extent and barbiturate sensitivity,” Journal of Neurophysiology, vol. 27, pp. 1154–1173, 1964. View at: Google Scholar
 M. N. Geffen, S. E. J. de Vries, and M. Meister, “Retinal ganglion cells can rapidly change polarity from off to on,” PLoS Biology, vol. 5, no. 3, article e65, 2007. View at: Publisher Site  Google Scholar
 J. De La Rocha, B. Doiron, E. SheaBrown, K. Josić, and A. Reyes, “Correlation between neural spike trains increases with firing rate,” Nature, vol. 448, pp. 802–806, 2007. View at: Publisher Site  Google Scholar
 P. Dayan and L. F. Abbott, Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems, MIT, Cambridge, Mass, USA, 1st edition, 2005.
 D. K. Warland, P. Reinagel, and M. Meister, “Decoding visual information from a population of retinal ganglion cells,” Journal of Neurophysiology, vol. 78, no. 5, pp. 2336–2350, 1997. View at: Google Scholar
 D. Peña, Análisis de Datos Multivariantes, McGrawHill Interamericana de España S.L., 2002.
 J. R. Manning, R. Ranganath, K. A. Norman, and D. M. Blei, “Topographic factor analysis: a bayesian model for inferring brain networks from neural data,” PLoS ONE, vol. 9, no. 5, Article ID e94914, 2014. View at: Publisher Site  Google Scholar
 B. M. Yu, J. P. Cunningham, G. Santhanam, S. I. Ryu, K. V. Shenoy, and M. Sahani, “Gaussianprocess factor analysis for lowdimensional singletrial analysis of neural population activity,” Journal of Neurophysiology, vol. 102, no. 1, pp. 614–635, 2009. View at: Publisher Site  Google Scholar
 J. H. Macke, L. Buesing, J. P. Cunningham, B. M. Yu, K. V. Shenoy, and M. Sahani, “Empirical models of spiking in neural populations,” in Advances in Neural Information Processing Systems 24, pp. 1350–1358, MIT Press, 2011. View at: Google Scholar
 X.J. Wang and G. Buzsáki, “Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model,” The Journal of Neuroscience, vol. 16, no. 20, pp. 6402–6413, 1996. View at: Google Scholar
 C. Echtermeyer, T. V. Smulders, and V. A. Smith, “Causal pattern recovery from neural spike train data using the Snap Shot Score,” Journal of Computational Neuroscience, vol. 29, no. 12, pp. 231–252, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 V. A. Klyachko and C. F. Stevens, “Connectivity optimization and the positioning of cortical areas,” Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 13, pp. 7937–7941, 2003. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2017 Jorge H. Soletta et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.