Identification of Functionally Interconnected Neurons Using Factor Analysis.

The advances in electrophysiological methods have allowed registering the joint activity of single neurons. Thus, studies on functional dynamics of complex-valued 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.


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][2][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 cross-correlation 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][9][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 well-known statistical analysis technique; however, it was not used to detect interconnections in the sensory systems based on response 2 Computational Intelligence and Neuroscience 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).

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 , (0 < < 1) 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 1 has connections with 1 , 2 , and 3 neurons, and it is weighted by 1 , 2 , and 3 , correlation values, respectively. Likewise, 2 current has connections with 4 and 5 neurons.
The correlated current values, , were inputs of Leaky Integrate-and-Fire (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 1 , 2 , 3 , 4 , and 5 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 feed-forward network [15]. Neural baseline activity is given by uncorrelated homogeneous Poisson processes. For different simulations, the rate parameter was varied ( = 1/10, 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 where = 1 and = 0.8 [16].

Factor
Analysis. The factor analysis (FA) is based on the following model [17]: where is a vector of × 1 which contain the spike count values of all neurons, is a mean vector ( × 1), is a vector of × 1 which contain the unobservable factors-spike count values of presynaptic neurons-Δ is a loading matrix of × , and is a vector of × 1 with unobservable perturbations, which has a multivariate normal distribution of order -(0, ), with the components of mean vector 0 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 tr( ) indicates the trace of matrix. The parameters that maximize the (Δ, ) convex function are obtained through their corresponding partial derivatives (Statistical toolbox, MATLAB).

Metrics.
The matrix 2-norm, , 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, Δ 0 , reflects the preestablished interconnections; in the case of Figure 1(a) it is as follows:   The black arrows indicate the direction of presynaptic current, while the red arrows indicate the neuron that generates the presynaptic current and whose activity is independent of the others. All other neurons to fire independently (without arrows).
Presynaptic neurons 1 2 3 4 5 Neurons records Neurons records 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 Table 1: (a) Average loading matrices obtained for the interconnection scheme of Figure 1(b) using the model of correlated current. These were calculated for different weighting values. Two hundred simulations were realized for each situation. 5 neuron activity is independent of the others (in bold font). (b) Idem (a) using neural network model. (c) Theoretical loading matrix for the interconnection scheme of Figure 1  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.

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 ( = 3 and = 1.5). Then, the identification of interconnections was quantitatively assessed through . Figure 3(a) shows values for interconnections scheme of Figure 1(a) and considering 3 = 0. It is possible to note that values decrease as 1 , 2 , 4 , and 5 values increase. This trend is most noticeable when 1 = 2 = 3 , Figure 3(b).

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.

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 Table 2: (a) Average loading matrices obtained for the interconnection scheme of Figure 1(c) using the model of correlated current. These were calculated for different weighting values. Two hundred simulations were realized for each situation. 2 neuron activity is independent of the others (in italic font) and 5 neuron turns out to be presynaptic to the others (in bold font). (b) Idem (a) using neural network model. (c) Theoretical loading matrix for the interconnection scheme of Figure 1(c).   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 5 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)). 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 5 neuron are the highest when comparing the elements of the first columns. This particularity indicates that 2 neuron would not be connected to other neurons and that it would fire independently (italic font values) and also that 5 is presynaptic to the other neurons (bold font values).
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 4 and 5 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 4 weighting. For this situation, it is not possible to determine the direction of the neural connection because of the similarity between the values of 4 and 5 neurons. Thus, the loading matrix only would indicate the existence of a common presynaptic neuron.
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).

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 = 0.8 value, while values for the other − 1 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.
In Tables 6(d), 6(b), and 6(f), it is observed that by increasing the number of presynaptic neurons values in the Table 3: (a) Average loading matrices obtained for the interconnection scheme of Figure 1(d) using the model of correlated current. These were calculated for different weighting values. Two hundred simulations were realized for each situation. 1 , 2 , and 3 fire independently (in italic font) and 5 neuron turns out to be presynaptic to the others (in bold font). (b) Idem (a) using neural network model. (c) Theoretical loading matrix for the interconnection scheme of Figure 1( Table 5: Theoretical loading matrix for the interconnection scheme of Figure 2(a).

Discussion and Conclusion
The FA is a multivariate statistical method that has traditionally been applied in the psychology areas and recently Computational Intelligence and Neuroscience  Table 7: Average loading matrices obtained for the interconnection scheme of Figure 2 Table 8: Theoretical loading matrix for the interconnection scheme of Figure 2  using an extension of FA (Gaussian-process 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][9][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.