[1] 12w Dimensional reduction of a V1 ring model with simple and complex cells

[1] 190w In this paper, we extend a framework for constructing low-dimensional dynamical systems models of mammalian primary visual cortex to a cortical network model that incorporates the full nonlinear effects of complex cells. The procedure consists of capturing the essential dynamics in a low-dimensional subspace using empirical methods, then recasting the equations in the reduced vector space. Previously, we considered visual cortical network models consisting of only simple cells with nearly linear responses to external stimuli. Here we show that fully nonlinear effects can be incorporated by examining the dimensional reduction of an idealized ring model of V1 with both simple and complex cells. We found it expedient to divide the subspace into four separate neuronal populations: excitatory simple, excitatory complex, inhibitory simple and inhibitory complex. In order to reproduce the fluctuation-driven dynamics in this reduced space, we incorporated (1) white noises with different intensities into individual neuronal populations, and (2) firing rate estimates to capture the probability of firing due to subthreshold fluctuations. With a more accurate, fitted connectivity, our modified dimensional reduced models can reproduce the firing rates, circular variances and modulation ratios observed in the original ring model.

[1] 112w Systems-level neuronal computations are often found to be emergent from the coherent activity of the underlying neural circuits. A major theoretical challenge in neuroscience modeling is to develop reduced descriptions of such dynamical activity and to capture the essential, network or circuit-wide computations. In recent papers, we have undertaken an empirical, data-driven approach to the dimensional-reduction of a large-scale, recurrently connected, neuronal network model of macaque primary visual cortex (V1). We have shown that in the cases where the dynamics is indeed coherent and lowdimensional, our approach provided a principled method to identify a reduced-dimensional system that can accurately model the dynamics of a large-scale network of neurons (Tao and Sornborger 2010).

[2] 103w Dimensional reduction techniques have found many applications in computational neuroscience. At the single neuron level, the FitzHugh-Nagumo equations represent a reduceddimensional version of the Hodgkin-Huxley (HH) neuron (FitzHugh 1961;Nagumo, Arimoto et al. 1962). Various versions (e.g., linear, quadratic, exponential) of the integrateand-fire neuron (I&F) have also been put forth as approximations of the HH neuron (Knight 1972;Hansel, Mato et al. 1998;Fourcaud-Trocmé, Hansel et al. 2003). At the network level, kinetic theoretical considerations have led to population density methods (Gerstner 2000;Knight, Omurtag et al. 2000;Nykamp and Tranchina 2000) and Fokker-Planck equations (Cai, Tao et al. 2006) describing coarse-grained versions of large-scale neuronal network dynamics.

[3] 230w From a data analysis point of view, the most common method for identifying suitable subspaces in multivariate data is the singular value decomposition (SVD). The SVD is a decomposition of a dataset into a set of orthogonal eigenvectors that captures covarying activity in the data. The usual application identifies a set of eigenvectors that represents most of the variance in the dataset, forming the basis for a reduced dimensional representation of the data. Dimensional reduction through identifying the appropriate mathematical subspaces of large dynamical systems has been extensively used in fluid turbulence (Sirovich 1987), weather modeling (Lorenz 1956) and molecular dynamics (Antoulas 2005). In the era of big data, and with the accelerating development of high-resolution neural imaging techniques, there arises an urgent need to develop systematic procedures that extract functional information directly from neural data. Multivariate analysis methods (e.g., SVD, independent component analysis, etc.) have been used effectively in the analysis of fMRI (Friston, Frith et al. 1995;McKeown, Jung et al. 1998), EEG (Subasi and Ismail Gursoy 2010), intrinsic optical imaging (Everson, Prashanth et al. 1998;Sornborger et al. 2003a, Sornborger et al. 2003b), calcium (Sornborger, Broder et al. 2008;Xu, Sornborger et al. 2008) and voltage-sensitive dye imaging (Sornborger et al. 2003a, Sornborger et al. 2003b). Here, we make use of these techniques to look for a mathematical description to bridge detailed cellular-level variables (unresolved or coarse-grained) and system-level variables.

[4] 230w Following the method of empirical eigenfunctions introduced by Sirovich and Rodriguez (Rodriguez and Sirovich 1990), we performed an SVD on the output of a large-scale numerical simulation, identifying the dominant eigenvectors and the subspace of relevance, and then performed a linear change of variables to project the original equations into the reduced-dimensional space. In the results presented in this paper, we used an idealized version of a model of an input layer of macaque V1 to extend our empirical, data-driven approach to a fully nonlinear network. Previously, we showed that our dimensionally-reduced models extracted from timeseries data were capable of predicting new input, that the computed couplings between systems-level variables could be used to reconstruct functional connectivities (Tao and Sornborger 2010), and that, by modeling the residuals with stochastic noise, we could reproduce the variance of firing rates in the original simulations (Tao, Praissman et al. 2012). These low-dimensional models were computationally efficient and could correctly reproduce a large-scale numerical simulation of simple cells in V1. Here we show that, by suitably accounting for the firing induced by subthreshold fluctuations, even a highly nonlinear network containing both simple and complex cells can be reproduced with our dimension-reduction framework. By empirically estimating an input-output relation (i.e., a firing rate function), we were able to accurately reproduce the distributions of circular variance, modulation ratios, and firing rates of the original largescale simulation.

[5] 72w By itself, empirical projection alone (2.2.1 and 2.2.2) cannot fully capture firing rates as observed in the original simulations. The reason was discussed in Tao, Praissman and Sornborger (Tao, Praissman et al. 2012): the projection and inverse-projection operators E T E, I T I, G T G in the DRM, having been truncated, are no longer identity matrices and fail to capture all synaptic fluctuations (and, most importantly, the spike-inducing synaptic fluctuations).

[6] 47w One way to compensate for the loss of fluctuation is to add noise to the conductances in the DRM. To do this, we assume the noise can be modeled by a spatially homogeneous, temporally uncorrelated, Gaussian process. First, residuals are calculated from the ring model simulation via

[7] 24w Then, standard deviations (SD) are calculated. Due to the different excitatory coupling strength of four populations, excitatory residuals were separated into four populations: σ-

[8] 22w Finally, the white noises ε E ; ε I ; ε G with different SDs were added to the conductances in DRM

[9] 73w DRMs added with noise are called stochastic dimensionally-reduced model (SDRM). In particular, Eq. ( 5) to ( 8) and (10) define our modified stochastic dimensionallyreduced model (Improved SDRM). Note that the white noises are orientation-averaged, that is to say, the SDs of the residuals were obtained from all neurons regardless of their preferred orientation. The orientation-averaged noise made neurons in the SDRM less orientation selective than the neurons in the original ring model.

[10] 171w the number of reduced dimensions Following Tao and Sornborger (Tao and Sornborger 2010), we use the Singular Value spectrum, the Lilliefors test (results not shown here), and the Akaike Information Crierion (AIC) to determine a suitable number of reduced dimensions. AIC is used to assess the most likely number of parameters in the DRM which can reproduce the original ring model output. As four populations of neurons have firing rate with different means and variances, we therefore obtain the residuals of firing rate separately for each population. The AIC is thus AIC M ð Þ ¼ ∑ P N P In RSS P =N P ð Þ ½ þ 2M . M represents the number of dimensionally-reduced variables, N P is the total number of time points of one population, P∈{ES,EC,IS,IC}, and the residual of one neuronal population is RSS P ¼ ∑ p;t ε 2 p;t , where ε p,t is the difference between the firing rate of the ring model and the DRM for neuron p at time t.

[1] 97w In the original model (Eq. ( 1) to ( 4)), individual spike trains were computed from the threshold passing times of the voltages. However, fine details of spike timing may be unimportant to modeling the average activity of the network and the dynamics can still be well-resolved on a time-scale of O 10 ð Þ msecs. Since the network contains complex cells, which can be fluctuation-driven, we coarse-grained Eq. ( 1) to (2) in time using an improved estimation of firing rate relative to Tao & Sornborger (Tao and Sornborger 2010) by taking subthreshold fluctuations into consideration:

[2] 137w are the effective reversal potential and the total conductance, respectively, and f P is a sixth degree polynomial function of a single neuronal population We introduce a sixth degree of polynomial to fit the V s -Rate curve when V s <1 in Eq. ( 5) (compared with no firing in our previous work, (Tao and Sornborger 2010;Tao, Praissman et al. 2012)). The reason for this modification comes from the fact that the neurons in our large-scale simulation may still cross threshold during the sampling time interval when the temporal average of V s <1. To see this point clearly, we plot the average firing rate versus the effective reversal potential in the original simulation (see Fig. 3 in Results). The mean firing rate is clearly non-negligible when V s is in the range of (0.5, 1).

[3] 53w In general, we can use Eq. ( 5) to provide an improved estimation of firing rate in our DRM. Even though the model we consider has time-varying sinusoidal LGN input, this estimation is reasonable on a timescale longer than the synaptic time-scales and shorter than the stimulus time-scale (Tao, Praissman et al. 2012).

[4] 63w In our study, we performed extensive comparison between the original ring model and its DRMs as constructed by the procedure we detailed above in Methods. The original network was stimulated by a sinusoidally-modulated grating drifting at 8 Hz with each orientation presented for 1 s. A typical trial had 32 orientations in a clockwise direction, i.e., 4 repeats for each of 8 orientations.

[5] 127w As the first step of dimensional reduction, we determined the optimal dimensions of reduced variables γ G , γ E , γ I in our DRM. We first examined the singular value spectrum of conductance. Figure 1 shows the semi-logarithmic singular value spectra of the LGN, intra-cortical excitatory and inhibitory conductances. (Semi-logarithmic plots of singular values for single neuronal population, e.g., excitatory simple, excitatory complex, inhibitory simple and inhibitory complex, have similar patterns.) As the square of each singular value represents the amount of variance in the dataset of the corresponding eigenvector, we note that the first 7 eigenvectors of G hold 94 % of the total variance, while the first 3 eigenvectors of E and I account for 91.1 and 93.7 % of the variances, respectively.

[6] 128w We make further confirmation of the dimension of intracortical conductances by using the AIC. The AIC estimates the tradeoff between the bias and variance of each model and determines the goodness of fit of our DRM in capturing the dynamics of the ring model. In Fig. 2, AIC is plotted as a function of M, the number of total variables of intra-cortical excitatory and inhibitory conductances trained with 16 grating directions. Here we sampled the original data at 80 Hz, fixed the number of LGN eigenvectors at 7, and thus M=2(n+1), where n represents the number of excitatory or inhibitory conductance eigenvectors. The minimum of the AIC curve is found at M=12, and we therefore choose n=5 as the final reduced dimension of intra-cortical excitatory and inhibitory conductances.

[7] 147w The integration time step in our DRM is of great importance for capturing the ring model dynamics. First, the time step should be shorter than the shortest decay time constants of post-synaptic conductances (here 5 ms for AMPA). Second, this time step is also the sampling interval for estimating firing rates. Thus this sampling interval should be long enough to avoid the costly step of repeatedly estimating instantaneous firing rates, and yet short enough to capture the response to time-varying stimuli. Here we choose 2.5 ms as the time step. In Fig. 3, we plot the relation between the effective reversal potential V s and the firing rate sampled at 400Hz. The blue lines are direct averages of firing rates from 64 s of ring model data; red dot dash lines are the expectation of firing rates using the mean-driven firing rate estimation (Shelley and McLaughlin 2002).

[8] 130w In Fig. 4, a sixth degree of polynomial is used to fit the V s -Rate curve of EC or IC cells when V s is less than 1, as discussed in Methods. In principle, the spiking activity in this regime is dependent on the mean synaptic conductances, but we find that the effect is weak. Therefore, we focus on modeling the complex cells because they are driven by large synaptic fluctuations and their responses are strongly nonlinear. We model the EC and IC cells separately. The polynomial functions for V s less than 1, together with a log function for V s larger than 1, constitute the modified estimation of coarsegrained firing rates (Eq. ( 5)). The DRMs which use this kind of modified estimation are called Improved DRMs.

[9] 85w Our Improved SDRM can reproduce firing rates, circular variance (CV) and modulation ratios F1/F0 in the ring model better than DRM or SDRM. In Fig. 5, we show scatter diagrams of neuronal firing rates between the original model and the DRMs. It is clear to see the improvement by adding noise and modifying the firing rate function when we compared the results between the original DRM (Tao and Sornborger 2010), the SDRM (Tao, Praissman et al. 2012), and our new DRM shown in Fig. 5.

[10] 190w Let us first examine the scatter diagrams between individuals cells in the original ring model and in the DRM for the ES and EC populations (Fig. 5a and b), respectively). Without appropriately accounting for synaptic fluctuations, the EC cells in DRM fire much less than in the original model. On the other hand, the ES cells in the DRM are less influenced by the loss of fluctuation and fire slightly more than those in the ring model. To understand this, we first add white noise only to the simple cells and find that their firing rates remain higher than in the original simulations. Then we add white noise to the complex cells and find firing rate of simple cells in the DRM match that of the original simulations. These two observations indicate that the truncation in the DRM leads to low firing rates of the complex cells. The main effect is that the low firing rates of the inhibitory complex cells lead in turn to low inhibitory conductances in the simple cells, causing the simple cells to have a higher firing rate in the DRM than in the original simulations.

[11] 103w Next, let us examine the scatter diagrams between individual cells in the original ring model and in the SDRM for the ES and EC populations (Fig. 5c and d), respectively). After adding noise, the firing rates in the SDRM are closer to the original simulations. Finally, we plot scatter diagrams of the firing rates of individual cells in the ring model and in the Improved SDRM for the ES and EC cells (Fig. 5e and f), respectively). The Improved SDRM much more closely follows the firing rate of the original ring model. The modified firing rate function does lead to an improved model.

[12] 111w Orientation selectivity is one of the major fundamental functional characteristics of V1 neurons. V1 is the first cortical area along the mammalian visual pathway where individual neurons show selectivity for stimulus orientation. Here we use the circular variance 1 (CV) (Ringach, Shapley et al. 2002) of mean firing rates to quantify the degree of selectivity. In Fig. 6, we plot the scatter diagrams of CV between neurons in the ring model and the DRMs. Again, we notice the improvement in the models as we incorporate the effects of synaptic fluctuations (either by modeling the residuals as noise or by using Eq. ( 5) to (6) to better capture the spike-inducing fluctuations).

[13] 85w 1 Orientation selectivity for drifting grating stimuli is measured by CV. Let m k denote the time-averaged firing rate with respect to stimulus angle θ k . The angles θ k spanned the range from 0 to 180°with equally spaced intervals. CV is defined as Notice that the main improvement comes from the selectivity of EC cells. The reason is that, in the original model, the complex cells are more fluctuation-driven and thus are more vulnerable to the neglect of synaptic fluctuations in the DRM.

[14] 179w Once we partially compensate for the fluctuations in the system, the firing rates, as well as the orientation selectivity, of complex cells in the DRM can be restored. Spatial summation is another fundamental property of V1 neurons. Hubel and Wiesel (Hubel and Wiesel 1962) described two major classes of cells (Simple and Complex) based on their receptive field properties. Later studies (Skottun, De Valois et al. 1991) (applying linear systems theory) showed that the key attribute that Hubel and Wiesel used to distinguish between the two classes was the linearity of summation: those neurons that showed nearly linear dependences were called Simple and the rest Complex. As many researchers used the drifting grating stimulus to probe linearity, it was proposed that a suitable explicit measure of linearity is the modulation ratio 2 (F1/F0) at optimal orientation, namely, the ratio between the amplitude of the first harmonic (F1) of the response and the mean response (F0) when the neuron was driven by a sinusoidal drifting grating at its optimal orientation. (The larger the first harmonic, the more linear the summation.)

[15] 60w Our DRM also reproduces this key property well. Figure 7 the scatter diagrams of F1/F0 for firing rates between neurons in the ring model and in the DRMs for the three different models. Again, we see the improvement, as we add fluctuations, either through modeling the residuals or through modeling the spike-inducing synaptic fluctuations using Eq. ( 5) to (6).

[16] 82w A key ingredient to the success of our DRM procedure lies in the effective connectivities. In our data-driven procedure, the fitted effective connectivities represent an empirically determined relationship between neuronal firing rates and the DRM system variables. Here we fit the connectivity matrices using a least-squares method, typically using 32 s of data from the original simulations. In Fig. 8, we show the relative error between original connectivity and effective connectivity as a function of sampling frequency. Relative error is defined as

[17] 25w where C and C are the real and fitted connectivity matrices, respectively. The matrix norm we use is the Frobenius norm, which is defined as

[18] 96w . We can see that the error in the excitatory connectivity matrix estimate drops to nearly zero when sampling frequency climbs above 40 Hz, whereas for the inhibitory connectivity matrix, the error drops when the sampling frequency is at or above 10 Hz. The difference in the location of the "knee" of the fitted excitatory and inhibitory connectivity matrices is due mainly to the different timescales in the respective synaptic conductances. Note that by using fitted connectivity matrices, we can build DRMs without the original cellular-level connectivities. Thus, our dimensional reduction procedure can be entirely data-driven.

[1] 73w In this paper, we showed that our dimensional reduction procedure is effective even when applied to a fully nonlinear neuronal network model of V1 containing both simple and complex cells. While we have come to expect that mean (i.e., coarse-grained) conductances and firing rates can be accurately represented in our DRM procedure, the results of Figs. 5 to 8 demonstrated that the orientation selectivity and modulation 2 Modulation ratio F1/F0 is defined as

[2] 143w , where R n is the cycle-averaged response to a sinusoidal drifting grating. Usually, F1/F0 is in the range from 0 to 1. and 2). Here, since our ring model contains four populations of neurons, the DRM models explicitly the corresponding populations. As may be expected, the original coarse-graining procedure indeed fails to capture all synaptic fluctuations, especially the important fluctuations that lead to complex cell-like behavior in the original simulation (Fig. 3, Fig. 5b, Fig. 6b, Fig. 7b). We thus estimated coarse-grained firing rates directly from the data. The improved SDRM reproduced firing rate, the CV and the modulation ratio F1/F0 in the original simulations (Fig. 5e and f, Fig. 6e and f, Fig. 7e and f). Finally, we showed that fitted connectivity matrices converge to original connectivity matrices when the sampling interval is smaller than the synaptic time constants (Fig. 8).

[3] 102w By modeling the sub-and near threshold voltage fluctuations that contribute to neuronal firing, we were able to extend the range of applicability of our original dimensional reduction procedure (Tao and Sornborger 2010). V1 complex cells are known to be highly nonlinear neurophysiologically (De Valois, Albrecht et al. 1982), and it was not clear that temporal coarse-graining would be successful at capturing the fluctuation sensitive nature of their dynamics. The coarse-graining (here, empirically determined as a function that explicitly transforms the voltage to neuronal spike rates) allowed us to neglect the details of spike timings that would have been costly to compute numerically.

[4] 166w In this and previous work, we have sought to use an empirical, data-driven method to the dimensional reduction of a neuronal network. In dynamics, they fitted the diffusion term in their Fokker-Planck equation. In this case and others, the detailed microscopic dynamics may be mapped to projected onto a known macroscopic equation (see, for instance, (Laing, Frewen et al. 2007)). Many efforts, equation-free or not, have been applied towards finding suitable lower-dimensional spaces and modeling the nonlinear dynamics within. Our procedure of modeling the subthreshold firing rates is along the same vein, and, as such, points to the generality of our dimensional reduction framework. While the SVD is a linear dimensional reduction method, the role that the SVD plays in our framework is primarily to finding a suitable reduced-dimensional subspace of the original dynamics; however, the dynamical equations of our DRMs within this lower-dimensional space are nonlinear. Thus our approach presents the possibility of dimensional reduction when a simple macroscopic description is not a priori available.

[5] 150w Our work has been aimed at a semi-automated data analysis framework for dimensional reduction of large-scale neuronal networks. The success of our procedure is based on the effective representation of the contribution of microscopic information to the dynamics of systems-level variables. In the DRM, the essential information is contained in the transformations that relate the membrane potential to the firing rate and in the inferred connectivities, here empirically computed from the data itself. (In Fig. 8, we show that the detailed, neuronal-level connectivity information is accurately reproduced.) In our DRM framework, the potential-firing rate transfer function and the inferred connectivities represent the input-output transformations at the synaptic level and provide a functional description between the signal integration at the cellular level and the dynamics at the network or systems level. Thus the dimensionally-reduced description of the neuronal connectivities allowed us to extract the functional couplings between the large-scale coherent modes.

[6] 215w The importance of modeling the subthreshold regime of the coarse-grained firing rates can be seen from Figs. 5 to 7. The success of our DRM in reproducing the CV and F1/F0 of the EC cells shows that an empirical estimate of the activity caused by subthreshold synaptic fluctuations allows our DRM to capture the essential dynamical features of the original network dynamics and, thereby, reproduce the functional aspects (specifically, the orientation selectivity and the degree of linearity in the spatial summation). However, work remains to be done to extend this procedure to networks containing other neuronal The mathematical framework that we have developed here and in previous work describes method by which a fully nonlinear neuronal network model may be constructed from firing rate data (and some information about the types of neurons and their synaptic conductances that make up the network) (see also, (Tao, Lauderdale et al. 2011)). Recent experiments using confocal or light-sheet microscopy in larval zebrafish transgenic for calcium indicators indicate that, soon, we will have (virtually) complete information concerning firing rates in an entire vertebrate brain. Such data in combination with our dimensional reduction methods should enable the construction of functionally accurate, large-scale neuronal network models, allowing us to make predictions and test hypotheses about network dynamics and their neural functions.

[1] 53w We follow Tao & Sornborger (Tao and Sornborger 2010) and Tao, et al. (Tao, Praissman et al. 2012) to perform dimensional reduction of a numerical model of V1. Here we briefly summarize the procedure and describe in detail the estimation of the subthreshold firing rate function and a new estimation of network connectivities.

[1] 50w Our visual cortical model consists of a system of N=1,024 conductance-based, integrate-and-fire (I&F) point neurons, forming a one-dimensional ring to model the activities of a single orientation hyper-column in mammalian primary visual cortex (V1). Under I&F dynamics, the temporal evolution of the membrane potential of each neuron is governed by

[2] 115w where P=E,I for excitatory and inhibitory neurons, j=1,…,N, V R , V E and V I are the rest, excitatory and inhibitory reversal potentials, respectively. We use dimensionless units except for time. The leakage conductance, g L =50/sec. We set V R =0 and use the difference between rest and threshold potentials to scale the reversal potentials, thus V E = 4.6667 and V I = -0.6667 (McLaughlin, Shapley et al. 2000). Once the j th neuron reaches the firing threshold and releases its l th spike, the exact time t l j is recorded and the membrane potential is reset to V R and held there for an absolute refractory period, τ ref .

[3] 19w We set τ ref =3 (1) msec for an excitatory (inhibitory) cell. The intra-cortical excitatory and inhibitory conductances are

[4] 45w We use drifting gratings as our visual stimulus. The j th simple cell in our ring model receives Poisson spike trains summed over all of its presynaptic LGN cells. Following (Shelley and McLaughlin 2002), we model the summed LGN excitation by its total firing rate

[5] 44w where ε represents the stimulus contrast and θ is the stimulus orientation. The preferred orientation and spatial phase of the j th neuron are θ j and ϕ j . The angular frequency ω=2πf, with drifting grating frequency f=8. The LGN conductance is then

[6] 10w where f 0 is the strength of each LGN spike.

[7] 82w The orientation preference of each neuron is given by θ j ¼ j N π . That is, the network has a one-dimensional ringlike architecture. The spatial phase preference of each neuron is randomly set to mimic the distribution observed in V1 (DeAngelis, Ghose et al. 1999). We assume that the LGN input is mediated only by AMPA receptors, and choose C= 800, ε =100 %, f 0 = 1. The term f INH j is a stimulus-nonspecific background inhibitory noise term.

[8] 237w We model the simple-complex property of individual neurons with the parameter λ j . That is, we assume simple cells receive strong LGN input (and set λ j to 1 in Eq. (2)), while complex cells, on the other hand, receive little or no LGN input (λ j = 0) (Hubel and Wiesel 1962) and are driven mainly by intra-cortical excitation (Chance, Nelson et al. 1999;Tao, Shelley et al. 2004). In Eq. ( 2), S PE 0 is the intra-cortical excitatory coupling strength for simple cells; S PE is the intra-cortical excitatory coupling strength for complex cells; S I is the intra-cortical inhibitory coupling strength. The values of the coupling strengths are S EE 0 =2, S EE =7.4, S IE 0 =4.5, S IE =9.5, S I =8. K j,k E (K j,k I ) denotes the excitatory (inhibitory) connectivity kernel between neurons j and k. The excitatory connectivity kernel is a normalized, two-dimensional Gaussian distribution with a standard deviation of 65°in orientation to model the local cortical excitation, while the inhibitory connectivity kernel is a normalized, two-dimensional uniform distribution to model untuned suppression (Xing, Ringach et al. 2011). The network is also sparsely connected: neuron k has probability p=0.15 to connect with neuron j. If neuron j and k are connected, p j,k =1; or else p j,k = 0. Post-synaptic conductances (PSCs) can be described by two coupled first order ordinary differential equations (ODEs):

[9] 32w where τ r (τ d ) are the rise (decay) time constants. We assume both AMPA and NMDA receptors contribute to the excitatory post-synaptic conductance (EPSC), thus G E take the form:

[10] 81w , where α stands for the fraction of AMPA receptors. We choose: α=0.5, τ r AMPA = 1 ms, τ d AMPA =5 ms, τ r NMDA =2 ms, τ d NMDA =80 ms. Inhibitory post-synaptic conductances (IPSC) are controlled by GABA A receptors and with time constants τ r GABA =1 ms, τ d GABA =10 ms. In addition, each neuron in our ring model also received an inhibitory noise conductance f INH j (t), which can be written as

[11] 19w with spike times drawn from a homogenous Poisson spike train with constant rate of 160/sec with S INH =0.14.

[12] 5w The improved stochastic dimensionally-reduced model

[13] 15w Following (Tao and Sornborger 2010), we perform a linear change of variables on the conductances:

[14] 51w where E, I and G are the projection matrices obtained from SVDs on time series data from the original ring model simulation. Dimensional reduction is achieved by truncating the number of neural eigenvectors of the full transformation. Note that g E =αg AMPA +(1-α)g NMDA , and we can further define

[15] 13w In the dimensionally reduced subspace, the system of ODEs can be written as

[16] 11w where τ AMPA,NMDA,I are the decay time constants of PSC and

[17] 62w Here k=1,…,n, where n is the number of eigenvectors kept in the truncation. We solve the above system of ODEs using Euler time-stepping. After one time step, we project back to the full, neuronal space to calculate the firing rate. The updated firing rate is again substituted into Eq. ( 7) and (8) for the next time step (Tao and Sornborger 2010).

[18] 27w Following Tao and Sornborger (Tao and Sornborger 2010), we fit the connectivity matrices using a least squares method. To begin, we write the intra-cortical excitatory conductance as

[19] 18w Here, we introduce the operator Φ to present the solution of Eq. ( 4). For a single spike,

[20] 34w where Θ(t) is the Heaviside step function. We write the excitatory connectivity matrix as C j,k E ≡S E j (p j,k /p)K j,k and the excitatory spike trains that neuron k receives as

[21] 14w . Then the second order excitatory input rate is defined in matrix form as

[22] 12w and the Eq. ( 11) can be rewritten in matrix form as

[23] 52w Note that we use m E because its integration over time is the same as time integration of the excitatory firing rate, m E . Similarly, the inhibitory connectivity matrix C I , inhibitory spike train δ I and inhibitory input rate m I are defined. Thus the intra-cortical inhibitory conductance is

[24] 54w Given g E ; g I ; m E ; m I from the ring model data, the fitted connectivity matrices C E ; C I can be obtained using a leastsquares method by Eq. ( 12) and (13). Normally, the fitted C E ; C I can be calculated in the form as

[25] 30w However, simulation data often leads to singularity of matrices m E m E T and m I m I T . In this case, we substitute pseudoinverses, giving for inverses

[26] 82w where M E and M I are orthogonal matrices obtained from SVDs performed on m E and m I . The numbers of eigenvectors in M E and M I are the ranks of matrices m E m E T and m I m I T . We can calculate m E and m I efficiently and accurately. Given the sampling frequency v and excitatory spike trains δ E k , we first integrate passively after one sampling interval Δt=1/v and get

[27] 26w r -e -Δt=τ d 8 < : If spikes δ E k occur during this interval, we update the variables h and m k E ,

[28] 7w Inhibitory m I may be calculated likewise.