Reconstructing neuronal circuitry from parallel spike trains

State-of-the-art techniques allow researchers to record large numbers of spike trains in parallel for many hours. With enough such data, we should be able to infer the connectivity among neurons. Here we develop a method for reconstructing neuronal circuitry by applying a generalized linear model (GLM) to spike cross-correlations. Our method estimates connections between neurons in units of postsynaptic potentials and the amount of spike recordings needed to verify connections. The performance of inference is optimized by counting the estimation errors using synthetic data. This method is superior to other established methods in correctly estimating connectivity. By applying our method to rat hippocampal data, we show that the types of estimated connections match the results inferred from other physiological cues. Thus our method provides the means to build a circuit diagram from recorded spike trains, thereby providing a basis for elucidating the differences in information processing in different brain regions.

neurons whose product of firing rates λ pre λ post is close to a given value (deviation less than 10%).
The jittering method rarely produces FPs but instead it produces a large number of FNs. Our GLMCC consistently produces a smaller number of errors in total than the two other estimation methods. (A) Excitatory connectivity: FPs represent directed links that were mistakenly assigned as excitatory, whereas FNs represent excitatory connections that were assigned as disconnected or inhibitory. (B) Inhibitory connectivity: FPs represent directed links that were mistakenly assigned as inhibitory, whereas FNs represent inhibitory connections that were assigned as disconnected or excitatory.  Data obtained from Leaky Integrate-and-Fire neurons 4,5 . We confirmed that the method estimates the connectivity accurately for these data as well.

Supplementary Note 1: Large-scale simulation using the NEST simulator
We conducted a large-scale simulation of a cortical network consisting of 10,000 neurons (See Supplementary Tables 1, 2, and 3

Supplementary Note 2: Derivation of the GLMCC from a two-body GLM
The GLMCC (equation 6 in the main text) can be derived from a two-body GLM describing individual neurons interacting with each other. This is given by the firing rates of two neurons, λ 1 (t) and λ 2 (t): where u i (t) represents extrinsic fluctuations for each neuron mediated through many marginal neurons. J ij is the neuronal connection from the jth neuron to the ith neuron.
s j (t) represents a temporal profile of monosynaptic impact: where t j k is the time of the kth spike of the jth neuron. The monosynaptic interaction of the timescale τ is modelled here with a fixed time profile: f (t) = exp(− t−d τ ) for t > d and f (t) = 0 otherwise, where d is the synaptic transmission delay.
If the precise time profiles of extrinsic fluctuations u 1 (t) and u 2 (t) are known, the mutual interactions J 12 and J 21 may be estimated by fitting equations (S1) and (S2) to a pair of spike trains. Detailed information for u 1 (t) and u 2 (t) may be available in some initial sensory systems. In an analysis of retinal neural networks, information on optical stimuli was efficiently utilized to estimate neuronal connectivity 8 . However, in cortical networks with denser connectivity, the temporal fluctuations may be induced not only extrinsically but also intrinsically and spontaneously because of mutual excitation between neurons 9,10 . Thus, fluctuating inputs to individual neurons u 1 (t) and u 2 (t) are generally unavailable. Though it is in principle possible to estimate external stimuli from the spike trains themselves, the temporal resolution cannot be finer than the typical inter-spike interval. This means that the influence of brain waves in subsecond time period cannot be estimated from spike trains of a few Hz.
Nevertheless, systematic fluctuations in subsecond periods are actually observed in CC, as demonstrated in Fig. 5B. This occurs because the CC is made by piling up spikes of one neuron that occurred in the neighbourhood of spikes of another neuron; accordingly, the firing rate of one neuron is in practice multiplied by the total number of spikes of another neuron. To make the most of this information, we derive another GLM that may adapt to the CC data.
A CC of spike trains generated from the underlying rates (S1) and (S2) is given as where the over-line · · · represents the time-average over s. By separating the time dependent fluctuations as u i (t) = u 0 i + δu i (t) and ignoring the higher order terms in the cumulant of CC, we obtain the GLMCC (equation(S5)): where a(t) represents the fluctuation induced by the external stimuli: Thus a(t) represents a large-scale modulation of the CC caused extrinsically and carried out by marginal neurons. The monosynaptic impacts J ij and J ji originally defined in the two-body GLM (equations (S1) and (S2)) appear in GLMCC (equation (S5)).

Supplementary Note 3: The Levenberg-Marquardt method
The Levenberg-Marquardt (LM) method, which interpolates between the Newton method and the gradient descent method works efficiently to maximize a function. Here, we apply this method to maximizing the log posterior distribution function log p(θ) ≡ log p(θ|{t k }) defined by the equation (10). The LM algorithm is given as where θ k is the kth iteration parameters, ∇ log p and H(log p) are the gradient and Hessian, The typical excitatory PSP (EPSP) and inhibitory PSP (IPSP) were calculated using a neuron model: where g eff = 0.27 mS/cm 2 and E eff = −67 mV are the effective conductance and resting voltage, respectively, and I AMPA/GABA (t) is a synaptic current mediated by either the AMPA or GABA receptor, respectively representing excitatory or inhibitory synapses. This model is known to reproduce the sub-threshold voltage trace of HH-type neurons and cortical neurons recorded in vitro 11,12 . Typical EPSP (mV) and IPSP (mV) values were well-fitted by the following equations ( Supplementary Fig. 8): On the basis of these empirical relations, the model conductances g AMPA and g GABA are translated into EPSP and IPSP, respectively.
If V (t−) < V th and V (t+) ≥ V th , 1. Set t * = t and V (t) = V res , and 2. Emit spike with time stamp t * .