Inference of synaptic connectivity and external variability in neural microcircuits

Abstract

A major goal in neuroscience is to estimate neural connectivity from large scale extracellular recordings of neural activity in vivo. This is challenging in part because any such activity is modulated by the unmeasured external synaptic input to the network, known as the common input problem. Many different measures of functional connectivity have been proposed in the literature, but their direct relationship to synaptic connectivity is often assumed or ignored. For in vivo data, measurements of this relationship would require a knowledge of ground truth connectivity, which is nearly always unavailable. Instead, many studies use in silico simulations as benchmarks for investigation, but such approaches necessarily rely upon a variety of simplifying assumptions about the simulated network and can depend on numerous simulation parameters. We combine neuronal network simulations, mathematical analysis, and calcium imaging data to address the question of when and how functional connectivity, synaptic connectivity, and latent external input variability can be untangled. We show numerically and analytically that, even though the precision matrix of recorded spiking activity does not uniquely determine synaptic connectivity, it is in practice often closely related to synaptic connectivity. This relation becomes more pronounced when the spatial structure of neuronal variability is jointly considered.

Appendix

1.1 Experimental methods

All procedures were carried out in accordance with the ethical guidelines of the National Institutes of Health and were approved by the Institutional Animal Care and Use Committee (IACUC) of Baylor College of Medicine.

The animals (n = 5 PV-Cre/Ai9 crosses on a C57Bl/6 background, labeled with the fluorescent marker tdTomato) with an age range of p40 to p60 were initially anesthetized with Isoflurane (3%) and then anesthesia was maintained by either Isoflurane (2%) or a mixture of Fentanyl (0.05 mg/kg), Midazolam (5 mg/kg), and Medetomidin (0.5 mg/kg) with anesthesia boosts consisting half of the initial dose every three hours. The body temperature of the animal was maintained at 37C throughout the surgery using a homeothermic blanket system (Harvard Instruments). In some experiments we applied eye oil ointment (polydimethylsiloxane) to prevent dehydration of the cornea. Surgery and dye injections of the Oregon Green 488 BAPTA-1 AM (OGB1, Invitrogen) calcium indicator were performed as previously described (Garaschuk et al. 2006).

We used stereotactic information to locate our recordings to the primary visual cortex of the mouse (V1). In some experiments we used intrinsic imaging to verify the location of V1 (Kalatsky and Stryker 2003). We recorded calcium traces using a custom built two-photon microscope equipped with a Chameleon Ti-Sapphire laser (Coherent) tuned at 800 nm and a 20x, 1.0 NA Olympus objective. Scanning was controlled by a custom built acousto-optic deflector system (AODs) (Cotton et al. 2013). The average power out of the objective was kept less than 120mW. Calcium activity was typically sampled at a mean of 260Hz (min/max: 78-450 Hz). We recorded data from depths of 100-540 μ m below the cortical surface.

The measured fluorescent traces were preprocessed in order to reduce common mode noise related to small cardiovascular movements (Cotton et al. 2013) and the firing rates were estimated using by nonnegative deconvolution (Vogelstein et al. 2012).

1.2 OU identity for normal $\bar {\mathbf {K}}$

For the OU system defined in Eq. (1), if KK^T = K^TK, G = gI, and QQ^T = σI we have

$$ \begin{array}{@{}rcl@{}} \mathcal{R}(0) &=& \frac{\sigma^{2} g^{2}}{{\tau^{2}_{r}}} \int\limits_{-\infty}^{0} {\exp}\left[ {\frac{1}{\tau_{r}} (\mathbf{I} - g \bar{\mathbf{K}}) s} \right]\\ &&\times {\exp} \left[ {\frac{1}{\tau_{r}} (\mathbf{I} - g \bar{\mathbf{K}})^{\mathrm{T}} s} \right] ds \\ &=& \frac{\sigma^{2} g^{2}}{{\tau^{2}_{r}}} \int\limits_{-\infty}^{0} {\exp} \left[ {\frac{1}{\tau_{r}} (2 \mathbf{I} - g \bar{\mathbf{K}} - g \bar{\mathbf{K}}^{\mathrm{T}}) s} \right] ds \\ &=& \sigma^{2} g^{2} \tau_{r} (2 \mathbf{I} - g \bar{\mathbf{K}} - g \bar{\mathbf{K}}^{\mathrm{T}})^{-1} \end{array} $$

$$ \Rightarrow \qquad \bar{\mathbf{K}} + \bar{\mathbf{K}}^{\mathrm{T}} = \frac{2}{g} \mathbf{I} - \sigma^{2} g \tau_{r} \mathcal{R}^{-1}(0) $$

1.3 Relation between AUROC and discriminability for normal distributions

It is well known that the area under an ROC curve may be parameterized and solved to yield the identity

$$ \text{AUROC} = \mathbb{P}(X_{1} > X_{0}) $$

for normally distributed scores in the positive and negative classes X₁ and X₀, respectively. If these scores are normally distributed then closure properties imply

$$ \mathbb{P}(X_{1} > X_{0}) = \mathbb{P}(X_{1} - X_{0} > 0) = \mathbb{P}(Y > 0) $$

$$ = 1 - \mathbb{P}(Y \leq 0) = 1 - F_{Y}(y) $$

$$ Y \sim \mathcal{N}(\mu_{1} - \mu_{0}, {\sigma^{2}_{1}} + {\sigma^{2}_{0}}) $$

and so

$$ \text{AUROC} = 1 - \frac{1}{2} \left[ 1 + \text{erf} \left( -\frac{\mu_{1} - \mu_{0}}{\sqrt{2} \sqrt{{\sigma^{2}_{1}} + {\sigma^{2}_{0}}}} \right) \right] $$

and so by appropriately defining the discriminability D and re-arranging we obtain

$$ \text{AUROC} = \frac{1}{2} \text{erfc} \left( -\frac{|D|}{\sqrt{2}} \right) $$

where the absolute value restricts the AUROC to $\left [ \frac {1}{2},1 \right ]$ which corrects for anti-classifiers. The result is the bijective mapping from D to AUROC, leading to the natural invertibility between the two which is necessary for our arguments.

1.4 Central limit of precision with independent external input

Precision values from the case of independent external input can be expressed element-wise as

$$ \mathbf{P}^{a b}_{\alpha \beta} = \mathbf{W}^{a b}_{\alpha \beta} + \mathbf{W}^{b a}_{\beta \alpha} - {\sum}_{\gamma = 1}^{N} \mathbf{W}^{c a}_{\gamma \alpha} {\mathbf{W}}_{\gamma \beta}^{c b} $$

which in the thermodynamic limit of network size ($N \rightarrow \infty$) leads the summation term to converge to the limiting normal distribution, so long as all elements of W are independent with finite variance. So long as we enforce normal distributions on and zero-variance gains on the synaptic strengths, and since the normal distribution is closed under summation, elements of P will also be normally distributed. All these assumptions hold under the standard ER case as well as the CER case, whereby correlations are specially constructed so that pairs of the form $\mathbf {W}^{c a}_{\gamma \alpha } {\mathbf {W}}_{\gamma \beta }^{c b}$ are independent across γ though not across pairs of α, β. The only breakdown of the CLT independence conditions occurs in the HD_out case, as all such pairs are now highly correlated across γ, giving an apparent power law limiting distribution instead.

While the above requirements on J hold exactly only if it is normally distributed, as long as Ω is Erdos-Renyi this theory will still hold as a fair approximation since the $\mathcal {O}(1)$ part of the precision distributions are determined by the summation term with the direct strengths offering only a $\mathcal {O}(1 \slash \sqrt {N})$ deviation. Even if synaptic weights are specified from a one-sided distribution of finite variance, the induced asymmetries against the limiting Gaussian will dissipate in the large N limit. Any variance present in the gains will also lead to small errors with an exact Gaussian, but these effects will again decay for large N and so the discriminability theory outlined in the paper should still hold as a good approximation for the expected AUROC.

1.5 Simulation and figure parameters

For all simulations, we use an alternative parameterization of synaptic strength which makes modulation easier to control. We express the average synaptic weights as

$$ \left[\begin{array}{ll} j_{ee} & j_{ei} \\ j_{ie} & j_{ii} \end{array}\right] = k_{1} \left[\begin{array}{ll} \psi_{ee} & \psi_{ei} \\ \psi_{ie} & \psi_{ii} \end{array}\right] $$

where the synaptic proportionsψ are normalized by the inhibitory component as ψ_ab = j_ab/j_ii and the mean synaptic magnitudek₁ is then modulated while holding the proportions fixed. Similarly, the synaptic variance is parameterized as

$$ \left[\begin{array}{ll} v_{ee} & v_{ei} \\ v_{ie} & v_{ii} \end{array}\right] = k_{2} \left[\begin{array}{ll} \psi^{2}_{ee} & \psi^{2}_{ei} \\ \psi^{2}_{ie} & \psi^{2}_{ii} \end{array}\right] $$

where the magnitude of synaptic variancek₂ is modulated. For all figures and simulations, we use a version of the synaptic proportions used in Pyle and Rosenbaum (2016) that have been perturbed in order to give non-zero real part to the eigenvalues. These proportions are

$$ \left[\begin{array}{ll} \psi_{ee} & \psi_{ei} \\ \psi_{ie} & \psi_{ii} \end{array}\right] = \left[\begin{array}{ll} 0.1 & 0.6 \\ 0.45 & 1 \end{array}\right] $$

as well as the fixed synaptic ratios of q_e = 0.8, q_i = 0.2 and the same fixed density of p_ab = p = 1.

In Fig. 1, we used σ_a = σ = 1, g_a = g = 1, and u_a = u = 0 for both rows in order to simplify the interpretation. For the top row we took k₁ = 2.5, k₂ = 6.25 × 10^− 5 and for the bottom row, k₁ = 12.5 and k₂ = 1.25. The same exact networks used in Fig. 1 are used in Fig. 2, though they are partitioned according to the more informative mask.

In Fig. 3, all shared network parameters are identical to Fig. 1, except k₁ = 2.5, k₂ = 1.25. The model specific parameters are ρ = 0.2 for the CER model and μ = 5, ξ = 0.25 for the HD models. As mentioned, σ for the HD model is estimated numerically using the bisection method to find the fixed point.

In Fig. 4, the same recurrent networks were used across all cases. Shared recurrent parameters are consistent with the top row of Fig. 1. Shared external input parameters are $p_{ax} = p_{x} = 0.1, j_{a x} = \psi _{a x} k_{x,1}, v_{a x} = \psi ^{2}_{a x} k_{x,2}, \psi _{e x} = 1.333, \psi _{i x} = 1$. Individually modified parameters are k_x,1 = 0.2571, k_x,2 = 0.6571,〈s_x, s_x〉_αα = r_x = 10,〈s_x, s_x〉_αβ = c = 0.1 for α≠β in both the full-rank correlated state and the correlated part of the first combination. The first combination also had Γ_ind = I for the independent part. The second combination used k_x,1 = 250, k_x,2 = 0, r_x = 10, c = 0, q_x = 0.2, Γ_ind = 10I. Full-rank effects in the correlated case are induced by taking q_x = 7 to give N_x = 14, 000 total external neurons, making Γ invertible with high probability even without regularization from the independent external source.

In Fig. 5, shared parameters are consistent with the top row of Fig. 1 and used a spatial width of ς_ab = ς = 0.2. The additional spatial variability inherent to these networks accounts for the difference in AUROC using only precision as a marginal metric between Fig. 2 and Table 1.

The membrane potential dynamics for the AdEx model are

$$ C_{m} \frac{d \mathbf{V}}{d t} = -g_{L} (\mathbf{V} - E_{L}) + g_{L} {\Delta}_{T} e^{\frac{\mathbf{V} - \mathbf{V}_{T}}{{\Delta}_{T}}} + \mathbf{T} - \mathbf{w} $$

$$ \tau_{w} \frac{d \mathbf{w}}{d t} = a (\mathbf{V} - E_{L}) - \mathbf{w} $$

subject to the rule that if a voltage exceeds the threshold ($\mathbf {V}^{a}_{\alpha }(t) \geq V_{th}$) then it returns to reset ($\mathbf {V}^{a}_{\alpha }(t) \rightarrow V_{re}$), its adaptation current is incremented ($\mathbf {w}^{a}_{\alpha } \rightarrow \mathbf {w} + b$), and a spike is recorded. The input terms stemming from recurrent sources R and external sources X may be expressed as

$$ \mathbf{T} = C_{m} (\mathbf{X} + \mathbf{R}) + \mathbf{Q} \ \frac{d\mathcal{W}_{t}}{dt} $$

$$ \mathbf{R}^{a}_{\alpha} = {\sum}_{c = e,i} {\sum}_{\gamma \in c} \mathbf{K}^{a c}_{\alpha \gamma} {\sum}_{n} \eta_{c} (t - t^{c,\gamma}_{n}) $$

$$ \mathbf{X}^{a}_{\alpha} = {\sum}_{\gamma \in x} (\mathbf{K}_{x})^{a}_{\alpha \gamma} {\sum}_{n} \eta_{x} (t - t^{x,\gamma}_{n}) $$

where $t_{n}^{c,\gamma }$ is the n^th spike time of neuron γ in population c = {e, i, x} and synaptic kinetics are modeled by the filter $\eta _{c}(t) = \frac {1}{\tau _{c}} e^{-\frac {t}{\tau _{c}}}{\Theta }(t)$ where Θ(t) is the Heaviside step function.

In Fig. 7, network parameters for the OU model are identical to the top row of Fig. 1 and have independent noise level σ = 0.1 and timescale τ_r = 1, discretized at the level of dt = 0.1. Shared network parameters for the EIF model are the same except for k₁ = 250 and k₂ = 0, as well as the fact that the statistics of the gains are no longer specifiable and are a consequence of the non-linearity of the system. AdEx-specific parameters are as follows: C_m = 1, g_L = 0.0667, E_L = − 72, V_th = − 50, V_re = − 75,Δ_T = 1, V_T = − 55, τ_w = 150, τ_e = 8, τ_i = 4, τ_x = 10, a = 0, and b = 0.1. External input followed: q_x = 0.2, k_x,1 = 250, k_x,2 = 0, and identical r_x, c, ψ, p_x as from Fig. 4 purple. The balanced network exhibited less than 20% relative error to both the balanced rate approximation and the mean-field covariance approximation from Baker et al. (2019).

All code pertaining to simulations and analysis may be found at https://github.com/cb239/Inference-of-Synaptic-Connectivity.