Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC5191923

Formats

Article sections

Authors

Related links

Nat Neurosci. Author manuscript; available in PMC 2017 April 30.

Published in final edited form as:

PMCID: PMC5191923

NIHMSID: NIHMS821164

Robert Rosenbaum,^{1,}^{2} Matthew A. Smith,^{3,}^{4,}^{5} Adam Kohn,^{6,}^{7} Jonathan E. Rubin,^{5,}^{8} and Brent Doiron^{5,}^{8}

Users may view, print, copy, and download text and data-mine the content in such documents, for the purposes of academic research, subject always to the full Conditions of use:
http://www.nature.com/authors/editorial_policies/license.html#terms

The publisher's final edited version of this article is available at Nat Neurosci

See other articles in PMC that cite the published article.

Shared neural variability is ubiquitous in cortical populations. While this variability is presumed to arise from overlapping synaptic input, its precise relationship to local circuit architecture remains unclear. We combine computational models and *in vivo* recordings to study the relationship between the spatial structure of connectivity and correlated variability in neural circuits. Extending the theory of networks with balanced excitation and inhibition we find that spatially localized lateral projections promote weakly correlated spiking, but broader lateral projections produce a distinctive spatial correlation structure: Nearby neuron pairs are positively correlated, pairs at intermediate distances are negatively correlated and distant pairs are weakly correlated. This non-monotonic dependence of correlation on distance is revealed in a new analysis of recordings from superficial layers of macaque primary visual cortex. Our findings show that incorporating distance-dependent connectivity improves the extent to which balanced network theory can explain correlated neural variability.

The spiking activity of cortical neurons is often characterized by their average response over a large number of trials, prompting a wealth of theoretical studies relating the structure of neuronal networks to their trial averaged firing rate dynamics^{1}. However, trial averages do not capture the stochastic and irregular dynamics characteristic of cortical populations and the nervous system in general^{2}. Indeed, trial-to-trial fluctuations are central to contemporary theories of cortical computation^{3,4}. A deep mechanistic understanding of neuronal variability remains an open challenge.

Early theoretical studies deduced that variable spiking activity could arise through a balancing of strong, yet opposing, excitatory and inhibitory synaptic inputs^{5,6}. Expanding on this conjecture, van Vreeswijk and Sompolinsky^{7} showed that networks of recurrently coupled model neurons robustly create a state where strong excitation is approximately balanced by inhibition, creating a push-pull dynamic that generates irregular spiking activity. More recently, balanced networks have been implicated in theories of optimal coding^{8}, working memory^{9}, and stimulus tuning^{10}. Numerous experimental studies have established that excitation is often approximately balanced by inhibition in cortical circuits^{11,12,13,14,15,16,17}. In sum, balanced networks provide a parsimonious model of the irregular spiking activity observed in cortical circuits.

Early balanced network models produced asynchronous activity through sparse connectivity^{7,18}. However, several experimental studies reveal that local cortical networks are densely connected, with connection probabilities between nearby neurons sometimes exceeding 40 percent^{19,20,21,22}. These data imply substantial overlap between local synaptic inputs, which could, in principle, synchronize cortical networks. However – counter to intuition – balanced networks with dense connectivity show weak spike train correlations^{23}. This “asynchronous state” results from the correlated excitatory (e) or inhibitory (i) afferents to neuron pairs being actively cancelled by a strong negative e-i correlation establishing weak correlations even when connectivity is not sparse^{23}.

Consistent with the predicted asynchronous state, some multi-unit extracellular recordings show noise correlations that are nearly zero on average^{24}. However, a majority of population recordings in cortex reveal comparatively large correlations^{25,26}. Several studies suggest that the magnitude of noise correlations is dependent on many factors^{27}, including arousal^{28}, attention^{29}, anesthetic state^{23,24,30,31} and cortical layer^{32,33}. Finally, while *in vivo* whole cell recordings reveal strong positive e-e and i-i correlations coexisting with strong e-i correlations^{13}, these correlation sources do not always perfectly cancel as predicted by some theoretical models^{28}. Taken together, these studies show that cortical circuits can exhibit both weak and moderate noise correlations, at odds with predictions from the current theory of balanced networks^{23}.

In this study, we generalize the theory of correlations in densely connected balanced networks to include the widely observed dependence of synaptic connection probability on distance^{34,21}. We show that spatially broad recurrent projections disrupt the asynchronous state, producing a signature spatial correlation structure: Nearby pairs of neurons are positively correlated on average, pairs at intermediate distances are negatively correlated and distant pairs are weakly correlated. These positive and negative correlations cancel so that the average correlation between pairs sampled randomly over a large range of distances is nearly zero. We uncover this non-monotonic dependence of correlation on distance in recordings from superficial layers of macaque primary visual cortex, but only after correcting for a latent source of shared fluctuations. Our findings decouple balanced excitation and inhibition from asynchronous network activity, greatly extending the applicability of balanced network theory to explain cortical dynamics.

We consider a network of excitatory and inhibitory exponential integrate-and-fire model neurons. Neurons provide recurrent, lateral synaptic input to one another and receive feedforward synaptic input from a non-local presynaptic population. A detailed mathematical analysis of correlations in the limit of large network size is provided in Supplementary Note S.1 and Supplementary Figures 1–3. Below, we provide an outline of these theoretical results and confirm their predictions using computer simulations. We first use a simplified network model to demonstrate how the asynchronous state considered in previous theoretical work^{23} is broken by heterogeneous input correlations^{35}. We then consider a more realistic model where neurons belong to a continuous spatial domain and connection probability depends on distance.

To demonstrate the mechanisms affecting correlations in recurrent networks, we first simulated a simplified network of *N*=20,000 neurons (half excitatory and half inhibitory) that all receive the same fluctuating feedforward input and are each connected with probability 0.25 (Figure 1a). Despite the fact that neuron pairs share all of their feedforward input and 25% of their recurrent synaptic input on average, spiking activity was asynchronous, with an average pairwise spike count correlation of 6.8×10^{−5} (Figure 1b,c).

This small average correlation is a defining characteristic of the asynchronous state. Mathematically, this state is realized when spike count covariances in the network satisfy^{23}

$${C}_{SS}~O(\mathit{1}/N)$$

where *C _{SS}* denotes the average spike count covariance between pairs of neurons in the recurrent network,

Our theoretical analysis proceeds by noting that spike count covariance is inherited from synaptic input covariance^{26} and therefore the two scale similarly with *N* in the asynchronous state,

$${C}_{II}~O(\mathit{1}/N)$$

where *C _{II}* denotes the average covariance between neurons’ synaptic inputs.

Synaptic inputs can be decomposed into their feedforward and recurrent sources, *I*=*F*+*R*, so that neurons’ input covariances decompose as

$${C}_{II}={C}_{FF}+{C}_{RR}+2{C}_{RF},$$

where *C _{FF}* is the average covariance between neurons’ feedforward input currents,

Shared input fluctuations are visualized by averaging the inputs to several neurons, so that the unshared contributions average out (Figure 1d). Overlapping inputs cause *C _{FF}* and

To study the impact of heterogeneity on correlations in balanced networks, we modified the above model by dividing the neurons into two populations. Each population received a separate feedforward input (Figure 1e). The two feedforward input sources were statistically identical but uncorrelated. Recurrent connectivity was not changed: Neurons were randomly connected without respect to population membership (identical to Figure 1a). This input heterogeneity dramatically changed the structure of correlations in the network. Pairs of neurons in the same population had strongly positive spike count correlations on average (0.34), while neuron pairs from opposite populations were negatively correlated with a nearly identical correlation magnitude (−0.34), and the average correlation between all pairs was nearly zero (4.2×10^{−4}; Figure 1f,g).

The mechanism responsible for this change in correlations can be understood by again separating the synaptic input covariance into its recurrent and feedforward sources, but generalizing the decomposition to account for neuron “distance” to obtain

$${C}_{II}(d)={C}_{FF}(d)+{C}_{RR}(d)+2{C}_{RF}(d).$$

(1)

Here, *C _{II}*(

In contrast, recurrent connections do not respect population membership and thus neither do the statistics of recurrent input,

$${C}_{RR}(0)={C}_{RR}(1)\phantom{\rule{0.16667em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}{C}_{RF}(0)={C}_{RF}(1).$$

Since covariances *C _{RR}*(

For illustrative purposes, we considered a simplified network model with discrete subpopulations, but correlations and connectivity in many cortical circuits depend on continuous quantities like physical distance or tuning similarity^{36,19,21}. Next, we generalize these findings to more biologically realistic networks with connection probabilities that depend on neuron distance.

We next considered a network of *N*_{e}=40,000 excitatory and *N*_{i}=10,000 inhibitory model neurons arranged on a square-shaped domain modeling a portion of a cortical layer. The neurons receive feedforward synaptic input from a separate layer of Poisson-spiking neurons and are connected with a probability that decays with distance (Figure 2). Specifically, the probability of a connection between two neurons in the recurrent network obeys

$$Pr(\text{connection})~g(d;{{\alpha}^{2}}_{\text{rec}})$$

where *d* is the distance between the neurons measured along the two-dimensional network, ~ denotes proportionality, *g*(*d;σ*^{2})~*exp*(−*d*^{2}/(*2σ*^{2})) is a Gaussian-shaped function and *α*_{rec} approximately represents the average length of a recurrent synaptic projection. Similarly, the probability of a synaptic projection from a neuron in the feedforward layer to a neuron in the recurrent layer decays with distance like a Gaussian with width parameter *α*_{ffwd}, where distance is measured parallel to the cortical surface (Figure 2). We next show that the asynchronous state is realized when *α*_{rec}<*α*_{ffwd}, then show that the asynchronous state cannot be realized when *α*_{rec}>*α*_{ffwd}.

As above, asynchronous spiking requires a cancellation between input covariances, *cf.*
Eq. (1), except that *d* now represents continuous instead of binary distance. Therefore, conditions on asynchrony require first an understanding of how input covariances depend on pairwise neuron distance.

Overlapping feedforward synaptic projections introduce *O*(1) correlations between the feedforward inputs to neuron pairs. Since nearby pairs share more feedforward inputs, these correlations are distance-dependent. Specifically, synaptic divergence causes feedforward input correlations to be twice as broad as synaptic projection widths (Figure 2b),

$${C}_{FF}(d)~g(d;{{\alpha}^{2}}_{\text{ffwd}})$$

The fact that *C _{FF}*(

$${C}_{RF}(d),{C}_{RR}(d)~g(d;{{\alpha}^{2}}_{\text{ffwd}})$$

(2)

The implications of this requirement on the spatial profile of spiking correlations are clarified by noting that recurrent synaptic input is generated by spike trains in the. Synaptic divergence causes the correlations between neurons’ recurrent synaptic inputs to be broader in space than the correlations between spike trains according to (Figure 2c)

$${{\sigma}^{\mathit{2}}}_{RR}={{\sigma}^{\mathit{2}}}_{SS}+{{\alpha}^{2}}_{\text{rec}}.$$

Here, *σ _{RR}* is the width of correlations between neurons’ recurrent synaptic input currents and

Correlations between recurrent inputs are constrained by the cancellation required in the asynchronous state. Specifically, Eq. (2) requires that the width of correlations between recurrent synaptic inputs satisfies

$${{\sigma}^{\mathit{2}}}_{RR}=\mathit{2}{{\alpha}^{2}}_{\text{ffwd}}.$$

Combining the two expressions for *σ _{RR}* above yields

$${{\sigma}^{\mathit{2}}}_{SS}=\mathit{2}({{\alpha}^{2}}_{\text{ffwd}}-{{\alpha}^{2}}_{\text{rec}})$$

(3)

The existence of a real solution to Eq. (3) requires that *α*_{ffwd}*>α*_{rec}; in other words, the spatial width of the recurrent projections must be narrower than the width of feedforward projections for the asynchronous state to exist. Further, Eq. (3) implies that *σ ^{2}_{SS}*>

To test these theoretical findings, we performed network simulations with feedforward synaptic projections broader than recurrent projections (Figure 3a). The simulations confirmed that *C _{RR}*(

As noted above, the cancellation between positive and negative correlations necessary for the asynchronous state cannot be realized when recurrent projections are broader than feedforward (*α*_{rec}>*α*_{ffwd}) since Eq. (3) cannot be solved in this case. Instead, neuron pairs inherit correlations from overlapping feedforward inputs so that *C _{SS}*(

This prediction was confirmed by numerical simulations identical to those discussed above, but with recurrent projections broader than feedforward (Figure 4a). As predicted, recurrent input correlations (Figure 4b, red and purple) were too spatially broad to cancel with the more sharply decaying feedforward correlations (Figure 4b, blue) so that the total input correlation between nearby neurons was large (Figure 4b, black; compare to Figure 3b). This effect introduced moderately strong correlations between nearby spike trains (Figure 4c,d,e) that do not decay to zero at increasing network size (Supplementary Note S.2 and Supplementary Figure 4e,f). Nevertheless, the network maintained excitatory-inhibitory balance (Figure 4f and Supplementary Figure 4a,b).

Since recurrent inputs must cancel feedforward inputs in balanced networks, *C _{RF}*(

Following the same argument made for the homogeneous network, the spike count correlations averaged over neuron pairs at all distances is *O(1/N)* (Figure 4e, dashed gray and Supplementary Note S.1). However, as noted above, the average correlation over each distance cannot be *O(1/N)*. Hence, there must be a cancellation between positive and negative correlations at different distances. As in Figure 1e–h, a competitive dynamic causes nearby neurons to be positively correlated and more distant neurons to be negatively correlated. This competitive dynamic does not extend beyond the reach of recurrent projections, so sufficiently distant neurons are weakly correlated. Hence, correlation decreases then increases with distance. This non-monotonicity can be explained more precisely using a mathematical theory of correlation transfer, (Figure 4e red dashed curve and Supplementary Note S.1). The heterogeneity of positive and negative correlations at different distances increases the standard deviation of pairwise correlations, but only modestly (stdev=0.16; Figure 4d, compare to Figure 3d).

In summary, when recurrent projections are spatially narrower than feedforward projections (*α*_{rec}<*α*_{ffwd} as in Figure 3), correlations are weak between pairs of neurons at all distances. When recurrent projections are broader than feedforward (*α*_{rec}>*α*_{ffwd} as in Figure 4), nearby neurons are positively correlated, neurons at moderate distances are negatively correlated and distant neurons are weakly correlated. Moreover, the average correlation between pairs of neurons sampled randomly at all distances is small. The non-monotonic dependence of correlation on distance is a distinct signature of correlations arising from broad recurrent projections. We next investigate whether this correlation structure predicted by our theory is present in cortical recordings.

We next asked whether our theoretical characterization of correlations in spatially extended networks can explain correlations in a cortical circuit. Layers 2/3 and layer 4C of macaque primary visual cortex (L2/3 and L4C) provide an ideal circuit for testing our predictions. Pairs of neurons in L2/3 exhibit moderately large noise correlations that decay with distance, but neurons in L4C, which are a primary source of inter-laminar input to L2/3, exhibit extremely weak pairwise noise correlations^{32,33} (Figure 5 with data from previous studies^{36,33}).

Neurons in L4C receive much of their feedforward input from thalamic projections, which form spatially broad synaptic fields, around 1 mm in diameter, but lateral projections within macaque L4C form narrower, sub-millimeter, synaptic fields^{34}. Thus, our theoretical prediction that correlations are weak when *α*_{ffwd}>*α*_{rec} is consistent with the weak pairwise correlations observed between L4C neurons in vivo (as in Figure 3).

Inter-laminar projections from L4C to L2/3 have a similar sub-millimeter width as excitatory intra-laminar projections within L4C and lateral projections from inhibitory basket cells in L2/3 form sub-millimeter synaptic fields, similar to those in L4C^{34}. Excitatory neurons in L2/3, however, form long-range lateral synaptic projections with synaptic fields spanning several millimeters in diameter^{34}. Our theoretical results can be generalized to this setting where inhibitory and excitatory projections have different spatial profiles (Supplementary Note S.1). This extension predicts the same correlation structure reported in Figure 4. However, correlations measured in L2/3 are positive on average over a broad range of distances^{36} (Figure 5b), in disagreement with this prediction.

We hypothesized that this inconsistency is explained by recent studies showing that much of the correlated variability measured in L2/3 arises from a low-dimensional shared source of latent variability^{38,31,30,39,40}. We conjectured that this shared variability increases pairwise correlations in L2/3 at all distances, thereby “washing out” the negative correlations predicted by our theory.

To search for low-dimensional variability in our data, we used Gaussian process factor analysis^{41,31} (GPFA), a statistical algorithm that extracts shared fluctuations from a population of spike trains (see Materials and Methods). Applying this algorithm to our L2/3 recordings reveals a source of one-dimensional covariability that decays with distance (Figure 5c). This distance-dependence implies that nearby neurons are affected similarly by the latent variable.

To test whether one-dimensional latent variability explains the discrepancy between our theoretical predictions and data, we built a two-layer network model representing a 10 mm by 10 mm square of cortex (Figure 6a). The first layer, representing L4C, was similar to the model in Figure 3 with the profile of feedforward and recurrent projections chosen to match experimentally constrained thalamic and lateral projection widths^{34}. The second layer, representing L2/3, was similar to Figure 4, with feedforward synaptic input from excitatory neurons in the L4C model and recurrent projection widths also chosen to match anatomical measurements. To capture latent variability in L2/3, the feedforward synaptic input to each neuron in the second layer was modulated by a time-varying, multiplicative gain modulation. We chose a multiplicative source of variability to be consistent with the properties of low-dimensional variability previously reported in macaque V1^{31}, but an additive source of latent variability would produce similar overall results. The gain modulation contributes an *O*(√*N*) source of covariance to the feedforward inputs that the recurrent network cannot cancel^{23}. To capture the distance-dependence of latent variability (Figure 5c), the magnitude of the gain modulation was heterogeneous across the network in such a way that nearby neurons received similar modulations and more distant neurons received less similar modulations.

Simulations of this two-layer model revealed that correlations between neurons in L4C were extremely small on average (Figure 6b,c gray), consistent with our theoretical predictions (Figure 3) and consistent with *in vivo* recordings (Figure 5). Correlations in the model L2/3 layer were moderately large and positive over all distances (Figure 6b,c black), comparable to in vivo recordings (Figure 5).

Thus, our model recovers the coarse structure of correlations in L4C and L2/3. However, our explanation of positive correlations in L2/3 is unsatisfying because the addition of globally shared variability destroyed the distinct non-monotonic relationship between correlation and distance predicted by our theory (compare Figure 4e to Figure 6b). We next asked whether this structure could be recovered by filtering out globally shared variability. To accomplish this, we computed the “residual correlation” matrix estimated by GPFA. Residual correlations approximate the spike count correlations with the contribution from low-dimensional variability removed^{31}.

Residual correlations computed between the simulated L2/3 spike trains exhibit the predicted non-monotonic dependence on distance, corroborating the ability of the GPFA algorithm to extract low-dimensional variability and leave the structure of residual correlations intact. We next computed the mean residual correlation in macaque L2/3 as a function of electrode distance. In doing so, we observed the same non-monotonic dependence of residual correlation on distance predicted by our theory (Figure 7b; further statistical analysis in Supplementary Note S.3 and Supplementary Figure 5).

In summary, combining theoretical analysis and computer simulations of a multi-layer network reveals a parsimonious model of the sources of shared variability in a visual cortical circuit in vivo. Under this model, positive correlations introduced by shared thalamic inputs to L4C neurons are actively canceled by negative correlations arising from recurrent circuitry so that pairs of L4C neurons at all distances exhibit weak average spike count correlations^{32,33}. Correlations between neurons in L2/3 are introduced by overlapping feedforward inputs from L4C and a low-dimensional source of variability. Correlations arising from overlapping L4C projections are filtered by recurrent circuitry in L2/3 to promote a non-monotonic dependence of correlation on distance. This non- monotonic correlation structure is washed out by low-dimensional latent variability, but can be recovered using GPFA to estimate and remove this variability.

Previous theoretical work on spatially homogeneous balanced networks with dense connectivity shows that they produce very weak spike train correlations ^{23}. We have generalized this theory to account for heterogeneous inputs and distance-dependent connection probability. In this framework we have made two important discoveries.

First, in agreement with the original findings, when lateral synaptic projections are spatially narrower than incoming feedforward projections, correlations are extremely weak on average at all distances. This theoretical finding can explain the weak pairwise correlations observed between neurons in middle layers of macaque primary visual cortex^{32,33}. However, correlations measured in cortical recordings are not always weak^{25}. Our second finding is that networks with broader lateral than feedforward projections produce correlations that do not decay to zero at increasing network size.

In previous studies of balanced networks with spatially homogeneous or clustered connectivity^{23,42}, the asynchrony condition (*C _{SS}~O(1/N)*) is satisfied and population averaged pairwise correlations vanish in the large network limit. In contrast, spatially extended networks with broad lateral projections violate the asynchrony condition, and consequently the expected pairwise correlations at a specific distance do not vanish. Nonetheless, mean excitatory and inhibitory currents balance and firing rates are moderate even when the asynchrony condition is violated (Supplementary Figure 4a,b and Supplementary Note S.1). This represents a novel solution for balanced networks that, for the first time, formally decouples network-wide asynchrony from excitatory-inhibitory balance.

We focused on the dependence of correlations on distance, but correlations also depend on tuning similarity. Partitioning L2/3 neuron pairs by tuning similarity reveals that correlations are strongest between similarly tuned neurons^{36}(Supplementary Figure 6a,b). Modifying our computational model to capture tuning-dependent correlations produces a non-monotonic dependence of residual correlation on tuning similarity in some parameter regimes, but the relevant parameters have not been measured experimentally (Supplementary Figure 6c–e and Supplementary Note S.4). Nevertheless, the modified theory could explain negative correlations previously observed in computer simulations of networks with tuning-specific connectivity^{32} and the finding that negative correlations are more frequent between disparately tuned neurons in V1^{43}.

As with nearly any computational model, many of the parameters used in our simulations may not reflect their corresponding values in specific cortical areas of specific species. However, our theoretical analysis does not depend on the precise values of these parameters. Our finding that the asynchronous state requires *α*_{rec}<*α*_{ffwd} is a fundamental property of networks with balanced excitation and inhibition.

We used a simplified model of a visual cortical circuit. In reality, pyramidal neurons in V1 form both local and long-range projections, connection probability in primate V1 depends on both distance and tuning similarity, and these dimensions are coupled^{44}. Moreover, connectivity properties of inhibitory neurons depend on their sub-type ^{45}. We modeled unidirectional connections from L4 to L2/3, but L4 also receives indirect feedback from L2/3 through deeper cortical layers. Spike trains in our model feedforward layer were modeled by homogeneous Poisson processes, in contrast to the oscillatory firing rates evoked by drifting grating stimuli in the data we analyzed. Our model can be extended to account for these additional features without affecting our overall conclusions.

Our findings have important implications for the interpretation of correlations in neural recordings. The average (residual) correlation between cell pairs sampled across a large range of distances could be extremely small, even when nearby pairs are positively correlated with moderate magnitude (Figures 4 and and7).7). Hence, subtracting low-dimensional latent variability and partitioning neuron pairs by distance can reveal correlation structure that would otherwise not be apparent. A previous study^{31} computed residual correlations as a function of distance in primate V1, but did not report a non-monotonic dependence. While we cannot be certain why their findings differ from ours, the accurate estimation of residual correlations with GPFA depends on the amount of data used to estimate shared variability. Our data are well-suited for this purpose, as they contain over 800 pairs of units per recording on average.

There is a long history of computational models of cortical circuits that consider either networks with spatially dependent coupling^{1} or balanced excitation and inhibition in spatially homogeneous networks^{7,23}. Only recently has the spatial structure of cortical connectivity been included in networks with balanced excitation and inhibition^{46,9,37}, and guiding theoretical principles are lacking. Our theory has taken this spatial structure into account and produced two core predictions for cortical circuits with long range lateral connections: First, nearby neurons exhibit significant positive correlations; second, the dependence of pairwise correlation on pairwise distance is non-monotonic. These predictions are clearly falsifiable and hence represent strong tests of our theory. The superficial layers of visual cortex have long-range lateral connections^{34}, making them a suitable test bed for our theory of correlations. After accounting for a source of global variability, both of our predictions were verified from population recordings in macaque V1 (Figure 7). Further, a similar noise correlation structure was reported in recordings from mouse V1^{47}. The successful validation of our predictions marks our theory as a promising framework for studying the structure of neural variability in cortical circuits. Nevertheless, there are many aspects of cortical dynamics that remain unexplained by balanced networks, such infrequent yet large membrane fluctuations during spontaneous dynamics^{15,48}. Capturing these dynamics in cortical models with balanced architectures remains an open challenge.

We modeled a square of cortex with *N* neurons, *N*_{e} of which are excitatory and *N*_{i} inhibitory. The membrane potential of neuron *j* from the excitatory (*a* = e) or inhibitory (*a* = i) population obeyed exponential integrate-and-fire (EIF) dynamics,

$${\mathit{C}}_{\mathit{m}}\frac{\mathit{d}{\mathit{V}}_{\mathit{j}}^{\mathit{a}}}{\mathit{dt}}={\mathit{I}}_{\mathit{L}}({\mathit{V}}_{\mathit{j}}^{\mathit{a}})+\mathit{f}({\mathit{V}}_{\mathit{j}}^{\mathit{a}})+{\mathit{I}}_{\mathit{j}}^{\mathit{a}}(\mathit{t});$$

each time that
${\mathit{V}}_{\mathit{j}}^{\mathit{a}}$ exceeds a threshold at *V*_{th}, the neuron spikes and the membrane potential is held for a refractory period *τ*_{ref} then reset to a fixed value *V*_{re}. The leak current is given by *I _{L}*(

Synaptic input currents were defined by

$${C}_{m}^{-1}\phantom{\rule{0.16667em}{0ex}}{I}_{j}^{a}(\mathit{t})={F}_{j}^{a}(t)+{R}_{j}^{a}(t)$$

where *F ^{a}*(

$${R}_{j}^{a}(t)={\sum}_{b=e,i}{\sum}_{k=1}^{{N}_{b}}\frac{{J}_{jk}^{ab}}{\sqrt{N}}{\sum}_{n}{\eta}_{b}(t-{t}_{n}^{b,k})$$

(4)

where
${t}_{n}^{b,k}$ is the *n*th spike time of neuron *k* in population *b* = e, i. The *1*/√*N* scaling of synaptic weights is a defining feature of the balanced network formalism and captures the balance between excitatory and inhibitory currents as well as intrinsically generated temporal variability for large *N*
^{7}. Each term
${J}_{jk}^{ab}$ represents the synaptic weight from presynaptic neuron *k* in population *b* to postsynaptic neuron *j* in population *a*. For all simulations, we modeled synaptic kinetics using *η _{b}*(

For the model in Figure 1, there were *N* = 20,000 neurons, half of which were excitatory and half inhibitory. For each (presynaptic) neuron, we randomly and uniformly chose 2500 excitatory and 2500 inhibitory postsynaptic neurons in the network. Postsynaptic neurons were sampled with replacement, so that a single presynaptic neuron can make multiple contacts with a postsynaptic neuron. The synaptic weight of each connection depended on the pre- and postsynaptic neuron types (excitatory or inhibitory). Specifically,

$${J}_{jk}^{ab}=(\#\phantom{\rule{0.16667em}{0ex}}\text{contacts})\times {j}_{ab}$$

(5)

where *j*_{ee} = 12.5mV, *j*_{ie} = 20 mV and *j*_{ii} = *j*_{ei} = −50 mV. Note that synaptic weights were scaled by √*N*= 141 in Eq. (4), so that the actual synaptic weight of each contact was on the order of 0.1mV.

For Figure 1a–d, the feedforward input to each neuron was given by the sum of an input bias and a smoothly varying signal,
${F}_{j}^{a}(t)=\sqrt{N}{m}_{a}+{\sigma}_{s}s(t).$Here, *s*(*t*) is a shared source of smooth, unbiased Gaussian noise defined by its auto-covariance function,

$$\mathit{cov}(s(t),s(t+\tau ))=exp(-{\tau}^{2}/2({\tau}_{s}^{2}))$$

*τ _{s}* = 40ms sets the correlation timescale and

To model the spatially extended recurrent network in Figures 3 and and4,4, we arranged *N*_{e} = 40,000 excitatory and *N*_{i} =10,000 inhibitory EIF model neurons on a uniform grid covering a two-dimensional square domain. The feedforward layer was modeled by a population of *N _{F}* = 5625 excitatory Poisson-spiking neurons on a uniform grid covering a square that is parallel to the recurrent network. Feedforward input to the recurrent layer was defined by

$${F}_{j}^{a}(t)=\sum _{k=1}^{{N}_{F}}\frac{{J}_{jk}^{aF}}{\sqrt{N}}\sum _{n}{\eta}_{e}(t-{t}_{n}^{F,k})$$

where
${t}_{n}^{F,k}$ is the *n*th spike time of neuron *k* in the feedforward layer. Each spike train in the feedforward layer was modeled as independent Poisson processes with rate *r _{F}* = 5 Hz.

To simplify calculations, we measured distance in units of the side-length of the square domain. In these units, the domain is represented as the unit square, Γ = [0, 1] × [0, 1]. Neurons were connected randomly and the probability that two neurons were connected depended on their distance measured periodically on Γ. The precise algorithm for generating connections is described in Supplementary Note S.5. This algorithm assures that the expected number of synaptic contacts from a presynaptic neuron at coordinates **y** = (*y*_{1}, *y*_{2}) in population *b* to a postsynaptic neuron at **x** = (*x*_{1}, *x*_{2}) in population *a* is given by

$${p}_{ab}(\mathbf{x}-\mathbf{y})=\frac{{K}_{ab}^{\mathit{out}}}{{N}_{a}}g({x}_{1}-{y}_{1};{\alpha}_{b})g({x}_{2}-{y}_{2};{\alpha}_{b})$$

where *g*(*u; α*) is a wrapped Gaussian distribution^{37}. Out-degrees were
${K}_{ee}^{\mathit{out}}={K}_{ei}^{\mathit{out}}=2000,{K}_{ie}^{\mathit{out}}={K}_{ii}^{\mathit{out}}=500,{K}_{eF}^{\mathit{out}}=10,000$ and
${K}_{eF}^{\mathit{out}}=800$. It follows that the network-wide average number of synaptic inputs to excitatory and inhibitory neurons in the recurrent network was *K* = 3715. Synaptic weights were determined by Eq. (5) where *j*_{ee} = 40 mV, *j*_{ie} = 120 mV, *j*_{ei} = *j*_{ii} = −400 mV, and *j*_{e}* _{F}* =

The first layer (L4C) in the model in Figure 6 was identical to the model in Figure 3 except that *α*_{ffwd} = 0.1, *α*_{e} = 0.05 and *α*_{i} = 0.03. The length units used in Figures 6 and and7a7a were determined by interpreting the network domain, Γ, as a 10 mm by 10 mm square. Thus, in physical dimensions, *α*_{ffwd} = 1 mm, *α*_{e} = 0.5 mm and *α*_{i} = 0.3 mm. Average firing rates in the L4C layer were 3.7 Hz for excitatory and 6.1 Hz for inhibitory neurons.

Connectivity in the second layer (L2/3) in Figure 6 was identical to the L4C layer except that *α*_{e} = 0.15 (or 1.5 mm), *α*_{ffwd} = 0.05 (or 0.5 mm). The spike times from Eq. (6) for the L2/3 layer were given by the spike times of neurons in the L4C layer, so that *N _{F}* = 50,000. The magnitude of feedforward connectivity was also modified by setting

$${F}_{j}^{a}(t)\leftarrow {F}_{j}^{a}(t)[1+{w}_{L}(\mathbf{x})L(t)].$$

The shared gain modulation, *L*(*t*), is a realization of unbiased Gaussian noise defined by its auto-covariance function,

$$\mathit{cov}(L(t),L(t+\tau ))=exp(-{\tau}^{2}/2({\tau}_{L}^{2}))$$

with correlation timescale *τ _{L}* = 40 ms. The dimensionless weight factor,

$${w}_{L}({x}_{1},{x}_{2})=0.5g({x}_{1}-c;{\sigma}_{L})g({x}_{2}-c;{\sigma}_{L})$$

where *c* = 0.5 and *σ _{L}* = 0.25 (in physical dimensions,

All simulations and numerical computations were performed on a MacBook Pro running OS X 10.9.5 with a 2.3 GHz Intel Core i7 processor. All simulations were written in a combination of C and Matlab (Matlab R 2014b, Mathworks). The differential equations defining the neuron model were solved using a simple forward Euler method with time step 0.1 ms and all simulations were run for a 22 s duration.

Anesthesia was induced with ketamine (10 mg/kg) and maintained during preparatory surgery with isoflurane (1.0–2.0% in 95% O2). Anesthesia during recordings was maintained with sufentanil citrate (6–18 *μ*g·kg^{−1}·h^{−1}). Vecuronium or pancuronium bromide (0.1–0.15 mg·kg^{−1}·h^{−1}) was used to suppress eye movements. Physiological signs were monitored to ensure adequate anesthesia and animal well-being. Vital signs (EEG, ECG, blood pressure, end-tidal P_{CO2}, temperature, and airway pressure) were monitored continuously. We used supplementary lenses to bring the retinal image into focus. At the end of the recording session, animals were killed and tissue was processed histologically to verify recording locations. All procedures were approved by the Institutional Animal Care and Use Committee of the Albert Einstein College of Medicine.

The pairwise correlations used to make the distributions in Figure 5a were computed using the same methods to those reported previously^{33} and include 1,613 L2/3 pairs and 469 L4 pairs. The data were recorded from nine anesthetized, adult male macaque monkeys (Macaca fascicularis). Recording were made with a group of five to seven linearly arranged (305 *μ*m spacing) platinum-tungsten electrodes or tetrodes (Thomas Recording), inserted normally to the cortical surface. The electrodes were advanced together through cortex, sampling in 200 *μ*m intervals until all electrodes had exited into white matter. Location of the middle layer sites involved a number of criteria including nominal depth, histology, and CSD analysis, as detailed previously^{33}.

Data for distance-dependent correlations reported in Figures 5b, 5c and and7b7b were from recordings described previously^{36}. The data were recorded from seven hemispheres of four adult male macaque monkeys (Macaca fascicularis). There was one experimental group of normal animals. The array consisted of a 10×10 grid of silicon microelectrodes (1 mm in length) spaced 400 *μ*m apart. We inserted the array partially into cortex, resulting in recordings confined mostly to layers 2–3. In two cases, we recorded simultaneously with a group of seven linearly arranged (2 mm extent) platinum-tungsten microelectrodes or tetrodes (Thomas Recording), positioned so that the nearest electrode was ~5 mm anterior to one edge of the multielectrode array. In this configuration, the distances between these electrodes and the array ranged from ~5 to 10mm. Neuronal receptive fields in all data sets were within 5° of the fovea. Waveform segments that exceeded a user-defined threshold were sorted offline (Plexon Offline Sorter). We quantified sort quality using the signal-to-noise ratio (SNR) of each candidate unit^{49}, keeping units with an SNR of at least 2.3 (for the data in Figure 5a) or at least 2.75 (for Figures 5b and and7b).7b). In both data sets, we eliminated neurons for which the best grating stimulus did not evoke a response of at least 2 spikes/s. Changing the SNR or responsivity threshold did not qualitatively change any of the results described herein.

Visual stimuli were generated with EXPO and displayed on a linearized CRT monitor (mean luminance 40 cd/m^{2}) with a resolution of 1,024 by 768 pixels and a refresh rate of 100 Hz. Stimuli were presented in a circular aperture surrounded by a gray field of average luminance. We mapped the spatial receptive field (RF) of units by presenting small (0.6°) drifting gratings at a range of spatial positions. We then centered our stimuli on the aggregate RF of the recorded units.

Stimuli were viewed binocularly and presented for 1.28 s, separated by 1.5 s intervals of isoluminant gray screen (except in one penetration, for which the interval was 10 s). We presented full-contrast drifting sinusoidal gratings at 8 or 12 orientations spaced equally (22.5 or 30° increments). The spatial frequency (1–1.3 cpd) and temporal frequency (3 or 6.25 Hz) values were chosen to correspond to the typical preference of parafoveal V1 neurons^{50}. The position and size (3.9–10°) of the grating were sufficient to cover the receptive fields of all the neurons. Stimulus orientation was block randomized, each stimulus orientation was repeated 100–200 times per recording session and each of the eight recording sessions yielded 21–68 units (210–2278 pairs) that met our SNR and responsivity thresholds. This yielded a total of 318 units and 6,907 simultaneously recorded pairs.

To compute spike count correlations for Figures 1, ,3,3, ,44 and and66 we first randomly sampled a subset of the excitatory neurons in the recurrent network. The number of neurons sampled is indicated in each figure caption. Neurons with firing rates less than 1 Hz were excluded from the correlation analysis. Pearson correlation coefficients were computed by counting spikes over 250 ms windows then computing the Pearson correlation coefficients between all pairs. The first 2 s of each simulation was excluded from the correlation analysis.

To obtain the residual correlations in Figure 7a, we randomly sampled 500 neurons with firing rates above 2 Hz from the L2/3 layer simulations. After omitting the first 2 s of the simulation, spike counts were computed using 50 ms bins, then partitioned into a sequence of 250 ms windows. To each window, we directly applied the GPFA algorithm used in a previous study^{31} and downloaded from http://toliaslab.org/publications/ecker-et-al-2014/. This is a modified version of the algorithm introduced previously^{41}. Briefly, GPFA extracts shared fluctuations from a population of spike trains by fitting spike count covariance matrices to a model of the form

$$cov(\mathbf{y})={\mathbf{cc}}^{T}+Q.$$

Here, **c** is a column vector of coefficients quantifying the contribution of latent variability to each neuron’s spiking activity. Latent variability is modeled by a one-dimensional Gaussian stochastic process. Hence, **cc ^{T}** is a rank-one matrix of latent covariances and

Spike count correlations in Figure 5b were computed across trials, counting spikes over 250ms time bins and subtracting the trial-averaged rate of each unit over each stimulus orientation. Hence, all reported spike count correlations represent correlated trial-to-trial variability (noise correlations) with stimulus correlations removed. Latent covariance in Figure 5c and residual correlation in Figure 7b were computed from the mean-subtracted spike counts using 50ms time bins and 250ms window sizes (treating each 250ms window as a separate trial) then using the same GPFA algorithms and procedures that were used for the simulated data in Figure 7a (see description above). The three p-values reported in Figure 7b were computed using one-sided, unpaired t-tests with t-values 8.1, 2.1 and 2.5 and degrees of freedom 3,937, 2,404 and 2,508 respectively.

The four error bars in Figure 3e were computed from n=880,057, 2,646,645, 4,420,219 and 4,550,579 pairs respectively. The five error bars in Figure 5b were computed from n=1,688, 2,251, 1948, 458 and 562 pairs respectively, and the same for Figures 5c and and7b.7b. The six error bars in Figure 7a were computed from n=3,885, 11,696, 19,493, 27,521, 35,336 and 26,819 pairs respectively. No statistical methods were used to pre-determine sample sizes but our sample sizes are the same as those reported in previous publications^{33,36}. Data collection and analysis were performed blind to the conditions of the experiments. No animals were excluded from the analysis.

Computer code for all simulations and analysis of the resulting data is included in Supplementary Software.

The data that support the findings of this study are available from the corresponding author upon request.

Click here to view.^{(559K, doc)}

Click here to view.^{(678K, pdf)}

Click here to view.^{(135K, pdf)}

This work was supported by National Science Foundation grants NSF-DMS-1517828 (RR), NSF-DMS-1313225 (BD), NSF-DMS-1517082 (BD), NSF-DMS-1612913 (JER), NSF-DMS-1516288 (JER) and NSF-DMS-1312508 (JER); National Institute of Health grants R01NS070865 (BD, JER), CRCNS-R01DC015139 (BD), R01EY016774 (AK), R01EY022928 (MAS) and P30EY008098 (MAS); two grants from the Simons Foundation collaboration on the global brain (SCGB#325293MC;BD, BD and 364994AK, AK); by the Eye and Ear Foundation of Pittsburgh (MAS); and by Research to Prevent Blindness (AK, MAS). We are grateful to Tai Sing Lee and Anthony Movshon for research support.

**Competing financial interests**

The authors declare no competing financial interests.

1. Ermentrout B. Neural networks as spatio-temporal pattern-forming systems. Reports on progress in physics. 1998;61:353.

2. Faisal A, Selen L, Wolpert DM. Noise in the nervous system. Nat Rev Neurosci. 2008;9:292–303. [PMC free article] [PubMed]

3. Pouget A, Beck JM, Ma WJ, Latham PE. Probabilistic brains: knowns and unknowns. Nature neuroscience. 2013;16:1170–1178. [PMC free article] [PubMed]

4. Shamir M. Emerging principles of population coding: in search for the neural code. Current opinion in neurobiology. 2014;25:140–148. [PubMed]

5. Softky WR, Koch C. The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J Neurosci. 1993;13:334–50. [PubMed]

6. Shadlen MN, Newsome WT. Noise, neural codes and cortical organization. Curr Op Neurobiol. 1994;4:569–79. [PubMed]

7. van Vreeswijk C, Sompolinsky H. Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science. 1996;274:1724–1726. [PubMed]

8. Deneve S, Machens CK. Efficient codes and balanced networks. Nature Neuroscience. 2016;19:375–382. [PubMed]

9. Lim S, Goldman MS. Balanced cortical microcircuitry for spatial working memory based on corrective feedback control. J Neurosci. 2014;34:6790–6806. [PMC free article] [PubMed]

10. Hansel D, van Vreeswijk C. The mechanism of orientation selectivity in primary visual cortex without a functional map. The Journal of Neuroscience. 2012;32:4049–4064. [PubMed]

11. Shu Y, Hasenstaub A, McCormick DA. Turning on and off recurrent balanced cortical activity. Nature. 2003;423:288–293. [PubMed]

12. Haider B, Duque A, Hasenstaub AR, McCormick DA. Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. J Neurosci. 2006;26:4535–4545. [PubMed]

13. Okun M, Lampl I. Instantaneous correlation of excitation and inhibition during ongoing and sensory-evoked activities. Nat Neurosci. 2008;11:535–537. [PubMed]

14. Dorrn AL, Yuan K, Barker AJ, Schreiner CE, Froemke RC. Developmental sensory experience balances cortical excitation and inhibition. Nature. 2010;465:932–936. [PMC free article] [PubMed]

15. Graupner M, Reyes AD. Synaptic input correlations leading to membrane potential decorrelation of spontaneous activity in cortex. The Journal of Neuroscience. 2013;33:15075–15085. [PMC free article] [PubMed]

16. Zhou M, et al. Scaling down of balanced excitation and inhibition by active behavioral states in auditory cortex. Nat Neurosci. 2014;17:841–50. [PMC free article] [PubMed]

17. Xue M, Atallah BV, Scanziani M. Equalizing excitation? inhibition ratios across visual cortical neurons. Nature. 2014;511:596–600. [PMC free article] [PubMed]

18. Amit D, Brunel N. Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb Cortex. 1997;7:237–252. [PubMed]

19. Ko H, et al. Functional specificity of local synaptic connections in neocortical networks. Nature. 2011;473:87–91. [PMC free article] [PubMed]

20. Fino E, Yuste R. Dense inhibitory connectivity in neocortex. Neuron. 2011;69:1188–1203. [PMC free article] [PubMed]

21. Levy RB, Reyes AD. Spatial profile of excitatory and inhibitory synaptic connectivity in mouse primary auditory cortex. J Neurosci. 2012;32:5609–5619. [PMC free article] [PubMed]

22. Oswald AM, Doiron B, Rinzel J, Reyes AD. Spatial profile and differential recruitment of GABAB modulate oscillatory activity in auditory cortex. J Neurosci. 2009;29:10321–10334. [PMC free article] [PubMed]

23. Renart A, et al. The Asynchronous State in Cortical Circuits. Science. 2010;327:587–590. [PMC free article] [PubMed]

24. Ecker A, et al. Decorrelated Neuronal Firing in Cortical Microcircuits. Science. 2010;327:584–587. [PubMed]

25. Cohen MR, Kohn A. Measuring and interpreting neuronal correlations. Nat Neurosci. 2011;14:811–819. [PMC free article] [PubMed]

26. Doiron B, Litwin-Kumar A, Rosenbaum R, Ocker GK, Josíc K. The mechanics of state-dependent neural correlations. Nature Neuroscience. 2016;19:383–393. [PubMed]

27. Kohn A, Zandvakili A, Smith MA. Correlations and brain states: from electrophysiology to functional imaging. Curr Opin Neurobiol. 2009;19:434–438. [PMC free article] [PubMed]

28. Poulet JF, Petersen CC. Internal brain state regulates membrane potential synchrony in barrel cortex of behaving mice. Nature. 2008;454:881–885. [PubMed]

29. Cohen MR, Maunsell JH. Attention improves performance primarily by reducing interneuronal correlations. Nat Neurosci. 2009;12:1594–1600. [PMC free article] [PubMed]

30. Mochol G, Hermoso-Mendizabal A, Sakata S, Harris KD, de la Rocha J. Stochastic transitions into silence cause noise correlations in cortical circuits. Proc Natl Acad Sci USA. 2015;112:201410509. [PubMed]

31. Ecker AS, et al. State dependence of noise correlations in macaque primary visual cortex. Neuron. 2014;82:235–248. [PMC free article] [PubMed]

32. Hansen BJ, Chelaru MI, Dragoi V. Correlated variability in laminar cortical circuits. Neuron. 2012;76:590–602. [PMC free article] [PubMed]

33. Smith MA, Jia X, Zandvakili A, Kohn A. Laminar dependence of neuronal correlations in visual cortex. J Neurophysiol. 2013;109:940–947. [PubMed]

34. Lund JS, Angelucci A, Bressloff PC. Anatomical substrates for functional columns in macaque monkey primary visual cortex. Cereb Cortex. 2003;13:15–24. [PubMed]

35. Wimmer K, et al. The dynamics of sensory integration in a hierarchical network explains choice probabilities in MT. Nat Commun. 2015;6:1–13.

36. Smith MA, Kohn A. Spatial and temporal scales of neuronal correlation in primary visual cortex. J Neurosci. 2008;28:12591–603. [PMC free article] [PubMed]

37. Rosenbaum R, Doiron B. Balanced Networks of Spiking Neurons with Spatially Dependent Recurrent Connections. Phys Rev X. 2014;4:021039.

38. Yu B, Kohn A, Smith MA. Estimating shared firing rate fluctuations in neural populations. Society for Neuroscience Meeting (SfN 2011); Washington, DC, USA. November 12–16, 2011; 2011.

39. Lin IC, Okun M, Carandini M, Harris K. The Nature of Shared Cortical Variability. Neuron. 2015:1–13. [PMC free article] [PubMed]

40. Scholvinck ML, Saleem AB, Benucci A, Harris KD, Carandini M. Cortical State Determines Global Variability and Correlations in Visual Cortex. J Neurosci. 2015;35:170–178. [PMC free article] [PubMed]

41. Yu BM, et al. Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. J Neurophysiol. 2009;102:614–635. [PubMed]

42. Litwin-Kumar A, Doiron B. Slow dynamics and high variability in balanced cortical networks with clustered connections. Nat Neurosci. 2012;15:1498–1505. [PMC free article] [PubMed]

43. Chelaru MI, Dragoi V. Negative correlations in visual cortical networks. Cerebral Cortex. 2016;26:246–256. [PMC free article] [PubMed]

44. Bosking WH, Zhang Y, Schofield B, Fitzpatrick D. Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex. J Neurosci. 1997;17:2112–27. [PubMed]

45. Pfeffer CK, Xue M, He M, Huang ZJ, Scanziani M. Inhibition of inhibition in visual cortex: the logic of connections between molecularly distinct interneurons. Nat Neurosci. 2013;16:1068–76. [PMC free article] [PubMed]

46. Kriener B, Helias M, Rotter S, Diesmann M, Einevoll GT. How pattern formation in ring networks of excitatory and inhibitory spiking neurons depends on the input current regime. Frontiers in Computational Neuroscience. 2014;7 [PMC free article] [PubMed]

47. Rikhye XRV, Sur M. Spatial Correlations in Natural Scenes Modulate Response Reliability in Mouse Visual Cortex. J Neurosci. 2015;35:14661–14680. [PMC free article] [PubMed]

48. Tan AY, Chen Y, Scholl B, Seidemann E, Priebe NJ. Sensory stimulation shifts visual cortex from synchronous to asynchronous states. Nature. 2014;509:226. [PMC free article] [PubMed]

49. Kelly RC, et al. Comparison of recordings from microelectrode arrays and single electrodes in the visual cortex. J Neurosci. 2007;27:261–264. [PMC free article] [PubMed]

50. De Valois RL, Albrecht DG, Thorell LG. Spatial frequency selectivity of cells in macaque visual cortex. Vis Res. 1982;22:545–559. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |