Possible solutions to the symmetric PME problem take the form of exponential functions characterized by two parameters, λ₁ and λ₂, which serve as Lagrange multipliers for the constraints:

The factor Z normalizes P to be a probability distribution.

By combining individual probabilities of events as given by Equation (2) the following relationship must be satisfied by any symmetric PME solution:

This is equivalent to the condition that the strain measure of Ohiorhenuan and Victor (2010) be zero (in particular, the strain is negative whenever p₃/p₀ − (p₂/p₁)³ < 0, a condition identified in Ohiorhenuan and Victor (2010) as corresponding to sparsity in the neural code).

For three-cell, symmetric networks, models that exactly satisfy Equation (3) will also be exactly described via PME. Moreover, note that probability models that meet this constraint fall on a surface in the space of (normalized) histograms, given by the probabilities p_j. One can verify by straightforward calculations (see Appendix) that—given fixed lower order moments—D_KL(P, P˜) is a convex function of the probabilities p_j. This has interesting consequences for predicting when large vs. small values of D_KL(P, P˜) will be found (see Appendix).

It is not necessary to assume permutation symmetry when deriving the PME fit P˜ to an observed distribution P, or in computing derived quantities such as D_KL(P, P˜), and we do not do so in this study. However, most of the distributions we study are derived from mechanistic models that are themselves symmetric or near-symmetric. Therefore, we anticipate that the simplified calculations for permutation-symmetric distributions will yield analytical insight into our findings.

We modeled a single ON parasol RGC in two stages (for details see section 4). First, we characterized the light-dependent excitatory and inhibitory synaptic inputs to cell k (g^exc_k(t), g^inh_k(t)) in response to randomly fluctuating light inputs s_k(t) via a linear-nonlinear model, e.g.,:

where N^exc is a static non-linearity, L^exc is a linear filter, and η^exc_k is an effective input noise that captures variability in the response to repetitions of the same time-varying stimulus. These parameters were determined from fits to experimental data collected under conditions similar to those in which PME models have been tested empirically (Shlens et al., 2006, 2009; Trong and Rieke, 2008). The modeled excitatory and inhibitory conductances captured many of the statistical features of the real conductances, particularly the correlation time and skewness (data not shown).

Second, we used Equation (4) and an equivalent expression for g^inh_k(t) as inputs to an integrate-and-fire model incorporating a non-linear voltage and history-dependent term to account for refractory interactions between spikes (Badel et al., 2007, 2008). The voltage evolution equation was of the form

where F(V, t − t_last) was allowed to depend on the time of the last spike t_last. Briefly, we obtained data from a dynamic clamp experiment (Sharpe et al., 1993; Murphy and Rieke, 2006) in which currents corresponding to g^exc(t) and g^inh(t) were injected into a cell and the resulting voltage response measured. The input current I_input injected during one time step was determined by scaling the excitatory and inhibitory conductances by driving forces based on the measured voltage in the previous time step; that is,

We used this data to determine F and C using the procedure described in Badel et al. (2007); details, including values of all fitted parameters, are described in section 4. Recurrent connections were implemented by adding an input current proportional to the voltage difference between the two coupled cells.

The prescription above provided a flexible model that allowed us to study the responses of three-cell RGC networks to a wide range of light inputs and circuit connectivities. Specifically, we simulated RGC responses to light stimuli that were (1) constant, (2) time-varying and spatially uniform, and (3) varying in both space and time. Correlations between cell inputs arose from shared stimuli, from shared noise originating in the retinal circuitry (Trong and Rieke, 2008), or from recurrent connections (Dacey and Brace, 1992; Trong and Rieke, 2008). Shared stimuli were described by correlations among the light inputs s_k. Shared noise arose via correlations in η^exc_k and η^ink_k as described in section 4. The recurrent connections were chosen to be consistent with observed gap-junctional coupling between ON parasol cells. We also investigated how stimulus filtering by L^exc and L^inh influenced network statistics. To compare our results with empirical studies, constant light, and spatially and temporally fluctuating checkerboard stimuli were used as in Shlens et al. (2006, 2009).

We start by considering networks without recurrent connectivity and with constant, full-field (i.e., spatially uniform) light stimuli. Thus, we set s_k(t) = 0 for k = 1, 2, 3, so that the cells received only Gaussian correlated noise η^exc_k and η^inh_k and constant excitatory and inhibitory conductances. Time-dependent conductances were generated and used as inputs to a simulation of three model RGCs. Simulation length was sufficient to ensure significance of all reported deviations from PME fits (see section 4). We found that the spiking distributions were strikingly well-modeled by a PME fit, as shown in the righthand panel of Figure 2A; D_KL(P, P˜) is 2.90 × 10⁻⁵ bits. This result is consistent with the very good fits found experimentally in Shlens et al. (2006) under constant light stimulation.

Next, we introduce temporal modulation into the full-field light stimuli such that each cell received the same stimulus, s_k(t) = s(t), where s(t) refreshed every few milliseconds with an independently chosen value from one of several marginal distributions. For our initial set of experiments, the marginal distribution was either Gaussian (as in Ganmor et al., 2011) or binary (as used in Shlens et al., 2006). For both choices, we explored inputs with a range of standard deviations (1/16, 1/12, 1/8, 1/6, 1/4, 1/3, or 1/2 of a baseline light intensity) and refresh rates (8, 40, or 100 ms). The shared stimulus produced strong pairwise correlation between conductances of neighboring cells. However, values of D_KL(P, P˜) remained small, under 10⁻² bits in all conditions tested.

We next asked whether PME models capture RGC responses to stimuli with varying spatial scales. We fixed stimulus dynamics to match the two cases that yielded the highest D_KL(P, P˜) under the full-field protocol: for both Gaussian and binary stimuli, we used 8 ms refresh rate and σ = 1/2. The stimulus was generated as a random checkerboard with squares of variable size; each square in the checkerboard, or stixel, was drawn independently from the appropriate marginal distribution and updated at the corresponding refresh rate. The conductance input to each RGC was then given by convolving the light stimulus with its receptive field, where the stimulus was positioned with a fixed rotation and translation relative to the receptive fields. This position was drawn randomly at the beginning of each simulation and held constant throughout (see insets of Figures 3B,C for examples, and section 4 for further details).

The RGC spike patterns remained very well described by PME models for the full range of spatial scales. Figure 3A shows this by plotting D_KL(P, P˜) vs. stixel size. Values of D_KL(P, P˜) increased with spatial scale, sharply rising beyond 128 μm, where a stixel had approximately the same size as a receptive field center, illustrating that introducing spatial scale via stixels produces even closer fits by PME models (the points at 512 μm correspond to the full-field simulations).

Values reported in Figure 3A are averages of D_KL(P, P˜) produced by five random stimulus positions. At stixel sizes of 128 μm and 256 μm, the resulting spiking distributions differed significantly from position to position; in Figure 3B, we show the probabilities of the distinct singlet [e.g., P(1, 0, 0)] and doublet [e.g., P(1, 1, 0)] spiking events produced at 256 μm. Each stimulus position created a “cloud” of dots (identified by color); large dots show the average over 20 sub-simulations. Each sub-simulation was identified by a small dot of the same color; because the simulations were very well-resolved, most of them were contained within the large dots (and hence not visible in the figure). Heterogeneity across stimulus positioning is indicated by the distinct positioning of differently colored dots. At smaller spatial scales, the process of averaging stimuli over the receptive fields resulted in spiking distributions that were largely unchanged with stimulus position, as shown in Figure 3C, where singlet and doublet spiking probabilities are plotted for 60 μm stixels. Thus, filtered light inputs were largely homogeneous from cell to cell, as each receptive field sampled a similar number of independent, statistically identical inputs; the inset of Figure 3C shows the projection of input stixels onto cell receptive fields from an example with 60 μm stixels. The resulting excitatory conductances and spiking patterns were very close to cell-symmetric (see Figures S2B,C).

By contrast, spiking patterns showed significant heterogeneity from cell to cell when the stixel size was large, as illustrated in Figure 3B. This arises because each cell in the population may be located differently with respect to stixel boundaries, and therefore receive a distinct pattern of input activity; this is illustrated by the inset of Figure 3B, which shows the projection of input stixels onto cell receptive fields from one such simulation. However, PME models gave excellent fits to data regardless of heterogeneity in RGC responses (see Figures S2E,F); as seen in Figure 3A, over all 20 sub-simulations, and over all individual stixel positions, we found a maximal D_KL(P, P˜) value of 0.00811.

Intrigued by the consistent finding of low values of D_KL(P, P˜) from the RGC model circuit despite stimulation by a wide variety of highly correlated stimulus classes, we sought to further characterize the processing of light stimuli by this circuit. In particular, we examined the effects of different marginal statistics of light stimuli, standard deviation of full-field flicker, and refresh rate on the marginal distributions of excitatory conductances. We focused on excitatory conductances because they exhibit stronger correlations than inhibitory conductances in ON parasol RGCs (Trong and Rieke, 2008).

With constant light stimulation (no temporal modulation) the excitatory conductances were unimodal and broadly Gaussian (Figure 2A, middle panel). For a short refresh rate (8 ms) or small flicker size (standard deviation 1/6 or 1/4 of baseline light intensity), temporal averaging via the filter L^exc and the approximately linear form of N^exc over these light intensities produced a unimodal, modestly skewed distribution of excitatory conductances, regardless of whether the flicker was drawn from a Gaussian or binary distribution (see Figures 2B,C, center panels). For a slower refresh rate (100 ms) and large flicker size (s.d. 1/3 or 1/2 of baseline light intensity), excitatory conductances had multi-modal and skewed features, again regardless of whether the flicker was drawn from a Gaussian or binary distribution (Figure 2D). Other parameters being equal, binary light input produced more skewed conductances. While some conductance distributions had multiple local maxima, these were never well separated, with the envelope of the distribution still resembling a skewed distribution.

The mechanism that leads to unimodal distributions of conductances, even when light stimuli are binary, is high-pass filtering—a consequence of the differentiating linear filter in Equation (7) and illustrated in Figure 2D. To demonstrate this, we constructed an alternative filter with a more monophasic shape [Equation (9), illustrated in Figure S1] and compared the excitatory conductance distributions side-by-side. We saw a striking difference in the response to long time scale, binary stimuli: the distributions produced by the monophasic filter reflected the bimodal shape of the input. Interestingly, the resulting simulation produced eight-times greater D_KL(P, P˜) (Figure 4). This suggests that greater D_KL(P, P˜) may occur when ganglion cell inputs are primarily characterized via monophasic filters, e.g., at low mean light levels for which the retinal circuit acts to primarily integrate, rather than differentiate over time.

In Figure 4A, we examine this effect over all full-field stimulus conditions by plotting D_KL(P, P˜) from simulations with the monophasic filter, against D_KL(P, P˜) from simulations in which the original filter was used with the same stimulus type. An increase in D_KL(P, P˜) was observed across stimulus conditions, with a markedly larger effect for longer refresh rates. This consistent change could not be attributed to changes in lower order statistics; there was no consistent relationship between the change in pairwise model performance and either firing rate or pairwise correlations (data not shown). Instead, large effects in D_KL were accompanied by a striking increase in the bi- or multi-modality of excitatory conductances (see Figure 4B). In Figure 4C, we show an example stimulus and excitatory current trace taken from the simulation shown in Figure 4B: the monophasic filter allows the excitatory synaptic currents to track a long-timescale, bimodal stimulus with higher fidelity, transferring the bimodality of the stimulus into the synaptic currents. This finding was robust to specifics of the filtering process; we were able to reproduce the same results by designing integrating filters in different ways (data not shown).

We next considered the role of recurrence in shaping higher-order interactions by incorporating gap junction coupling into our simulations. We did this separately for each full-field stimulus condition described earlier. In each case, we added gap junction coupling with strengths from 1 to 16 times an experimentally measured value (Trong and Rieke, 2008), and compared the resulting D_KL with that obtained without recurrent coupling (Figure 5).

At the experimentally measured coupling strength (g^gap = 1.1 nS) itself, the fit of the pairwise model barely changed (Figure 5A) from the model without coupling. At twice the measured coupling strength (g^gap = 2.2 nS), recurrent coupling had increased higher-order interactions, as measured by larger values of D_KL for all tested stimulus conditions. Higher order interactions could be further increased, particularly for long refresh rates (100 ms), by increasing the coupling strength to four or eight times its baseline level (g^gap = 4.4 nS or g^gap = 8.8 nS; see Figures 5B,C). Consistent with the intuition that very strong coupling leads to “all-or-none” spiking patterns, D_KL(P, P˜) decreased as g^gap increased further, often to a level below what was seen in the absence of coupling (Figure 5D). In summary, the impact of coupling on D_KL is maximized at intermediate values of the coupling strength. However, the impact of recurrent coupling on the maximal values of D_KL evoked by visual stimuli is small overall, and almost negligible for experimentally measured coupling strengths.

Finally, we repeated the full-field, recurrent, and alternate filter simulations previously described with light stimuli drawn from either Cauchy or heavy-tailed distributions: such distributions have been found to model the frequency of occurrence of luminance values in photographs of natural scenes (Ruderman and Bialek, 1994). In contrast to previous results with Gaussian and bimodal inputs, here we found very low D_KL(P, P˜) over all stimulus conditions: the largest values found were more than an order of magnitude smaller than those obtained earlier. Specifically, for all conditions, we found D_KL(P, P˜) < 4.5 × 10⁻⁴, over all 42 network realizations; for many simulations, this number did not meet a threshold for statistical significance (see section 4.1.7), indicating that P and P˜ were not statistically distinguishable. Using a more monophasic filter resulted in no apparent consistent change to D_KL(P, P˜). When gap junction coupling was added, D_KL(P, P˜) was maximized at an intermediate value; when g^gap = 8.8, all simulations produced a statistically significant D_KL(P, P˜) ≈ 3 − 4 × 10⁻³. However, overall levels remained relatively low, roughly 1/2 the value achieved with Gaussian or binary stimuli.

To explain these findings, we examined the excitatory input currents: we found that over a broad range of refresh rates and stimulus variances, the marginal distributions of excitatory input conductances produced were remarkably unimodal in shape, and showed little skewness (Figure 6A). By examining the time evolution of the filtered stimuli (see Figure 6B), we see that heavy-tailed distributions allow rare, large events, but at the expense of medium-size events which explore the full range of the linear-nonlinear model used for stimulus processing (compare the blue with the red/green traces). When combined with the Gaussian background noise, this produces near-Gaussian excitatory conductances and, as may be expected from our original full-field simulations, very low D_KL.

We hypothesize that the methodology of averaging over the entire stimulus ensemble may not capture the significance of rare events that may individually be detected with high fidelity: D_KL was low even for full-field, high variance stimuli, which presumably caused (infrequent) global spiking events. Additionally, an important avenue for future work would be to test the ability of our RGC model, which was trained on Gaussian stimuli, to accurately model the response of a ganglion cell to stimuli whose variance is dominated by large events. Recent work examining the adaptation of retinal filtering properties to higher-order input statistics found little evidence of adaptation; however, the stimuli used in this work incorporated significant kurtosis but not heavy tails (Tkacik et al., 2012).

In summary, we probed the spiking response of a small array of RGC models to changes in light stimuli, gap junction coupling, and stimulus filtering properties, and identified two circumstances in which higher-order interactions were robustly generated in the spiking response. First, higher-order interactions were generated when excitatory currents had bimodal structure; we observed such structure when bimodal light stimuli was processed by a relatively monophasic filter. Secondly, higher-order interactions were maximized at an intermediate value of gap junction coupling; this value was, however, much larger (eight times) than the experimentally observed coupling strength.

In the previous section, we developed results for a computational model tuned to a very specific cell type; we now ask whether these findings will hold for a more general class of neural circuits, or whether they are the consequence of system-specific features. To answer this question, we considered a simplified model of neural spiking: a feedforward circuit in which three spiking cells sum their inputs and spike according to whether or not they cross a threshold. Such highly idealized models of spiking have a long history in neuroscience (McCulloch and Pitts, 1943) and have been recently shown to predict the pairwise and higher-order activity of neural groups in both neural recordings and more complex dynamical spiking models (Nowotny and Huerta, 2003; Tchumatchenko et al., 2010; Yu et al., 2011; Leen and Shea-Brown, 2013).

In more detail, each cell j received an independent input I_j and a “triplet”—(global) input I_c that is shared among all three cells. Comparison of the total input S_j = I_c + I_j with a threshold Θ determined whether or not the cell spiked in that random draw. An additional parameter, c, identified the fraction of the total input variance σ² originating from the global input; that is, c ≡ Var[I_c]/Var[I_c + I_j]. The global input was chosen from one of several marginal distributions, which included those used in the RGC model: Gaussian, bimodal, and heavy-tailed. The independent inputs I_j were, in all cases, chosen from a Gaussian distribution, consistent with our RGC model. When the common inputs are Gaussian, our model is equivalent to the Dichotomized Gaussian model previously studied by several groups (Amari et al., 2003; Macke et al., 2009, 2011; Yu et al., 2011), cf. (Tchumatchenko et al., 2010). For further details, see section 4.2.

In the RGC model large effects in D_KL were accompanied by a striking increase in the bi- or multi-modality of excitatory conductances. Why are bimodal inputs, shared across cells, able to produce spiking responses that deviate from the pairwise model? We use our simple thresholding model to provide some intuition for how bimodal common inputs to thresholding cells lead to spiking probabilities that violate the constraints (Equation 3) which must hold for the pairwise model. For example, suppose that the common input I_c can take on values that cluster around two separated values, μ_A < μ_B, but rarely in the interval between; that is, the distribution of I_c is bimodal. If μ_B is large enough to push the cells over threshold but μ_A is not, then we see that any contribution to the right-hand side of Equation (3), p₂/p₁, depends only on the distribution of the independent inputs I_j; if either one or two cells spike, then the common input must have been drawn from the cluster of values around μ_A, because otherwise all three cells would have spiked.

To be concrete, let P[x] refer to the probability of spiking event x = (x₁, x₂, x₃), and P[x | I_c ≈ μ_A] refer to the probability that x occurs, conditioned on the event I_c ≈ μ_A. Then

because P[(1, 0, 0) | I_c ≈ μ_B] = 0: for the same reason,

therefore

On the other hand,

By changing the relative likelihood of drawing the common input from one cluster or the other, without changing the values of μ_A and μ_B themselves (that is, change P[I_c ≈ μ_B] and P[I_c ≈ μ_A] but leave the conditional probabilities (e.g., P[(1, 0, 0) | I_c ≈ μ_A]) fixed) one may change the ratio p₃/p₀ without changing the ratio p₂/p₁. Hence the constraint specifying those network responses exactly describable by PME models can be violated when the common input is bimodal.

In contrast, we may instead consider a unimodal common input, of which a Gaussian is a natural example. Here, the distribution of the common input I_c is completely described by its mean and variance; both parameters can impact the ratio p₃/p₀ (by altering the likelihood that the common input alone can trigger spikes) and the ratio p₂/p₁. Each value of I_c is consistent with both events p₁ and p₂, with the relative likelihood of each event depending on the specific value of I_c; it is no longer clear how to separate the two events. In the following sections, we will confirm this intuition by direct evaluation of the resulting departure from pairwise statistics.

Motivated by our observations of excitatory currents that arose in the RGC model, we chose several input distributions that allow us to explore other salient features, such as symmetry and the probability of large events. A distribution is called sub-Gaussian if the probability of large events decays rapidly with event size, so that it can be bounded above by a scaled Gaussian distribution (see section 4). We considered two sub-Gaussian distributions; the Gaussian itself, and a skewed distribution with a sub-Gaussian tail (hereafter referred to as “skewed”). We also considered the two “heavy-tailed” distributions used as stimuli to the RGC model—the Cauchy distribution, and a skewed distribution with a Cauchy-like tail (hereafter referred to as “heavy-tailed skewed”). In these distributions, the probability of large events decays polynomially rather than exponentially.

For each choice of common input marginal, we varied the input parameters so as to explore a full range of firing rates and pairwise correlations: specifically, we varied the input correlation coefficient c in the range [0, 1], the total input standard deviation σ in the range [0, 4], and the threshold Θ in [0, 3]. In all cases the independent inputs I_j were chosen from a Gaussian distribution [of variance (1 − c) σ²]. For each choice of input parameters, we determine the resulting distribution on spiking states (as described in section 4) and compute the PME approximation.

We first considered common inputs chosen from a unimodal (e.g., Gaussian) distribution. If I_c is Gaussian, then the joint distribution of S = (S₁, S₂, S₃) is multivariate normal, and therefore characterized entirely by its means and covariances. Because the PME fit to a continuous distribution is precisely the multivariate normal that is consistent with the first and second moments, every such input distribution on S exactly coincides with its PME fit. However, even with Gaussian inputs, outputs (which are now in the binary state space {0, 1}³) will deviate from the PME fit (Amari et al., 2003; Macke et al., 2009). As shown below, non-Gaussian unimodal inputs can produce outputs with larger deviations. Nonetheless, these deviations are small for all cases in which inputs were chosen from a sub-Gaussian distribution, and PME models are quite accurate descriptions of circuits with a broad range of unimodal inputs.

We first considered circuits with either Gaussian or skewed common inputs. Over the full range of input parameters, distributions remained well fit by the pairwise model, with a maximum value of D_KL(P, P˜) (of 0.0038 and 0.0035 for Gaussian and skewed inputs, respectively) achieved for high correlation values and σ comparable to threshold. In Figure 7A we illustrate these trends with a contour plot of D_KL(P, P˜) for a fixed value of threshold (here, Θ = 1.5) and Gaussian common inputs (the analogous plot for skewed inputs is qualitatively very similar, Figure S3A).

Clear patterns also emerged when we viewed D_KL(P, P˜) as a function of output spiking statistics rather than input statistics (as in Macke et al., 2011). Non-linear spike generation can produce substantial differences between input and output correlations; this relationship can vary widely based on the specific non-linearity (Moreno et al., 2002; de la Rocha et al., 2007; Marella and Ermentrout, 2008; Shea-Brown et al., 2008; Vilela and Lindner, 2009; Barreiro et al., 2010, 2012; Tchumatchenko et al., 2010; Hong et al., 2012). Figure 7B shows D_KL(P, P˜) and Δ for all threshold values (including the data shown in Figure 7A), but now plotted with respect to the output firing rate. The data were segregated according to the Pearson's correlation coefficient ρ between the responses of cell pairs (ρ≡Cov(xi,xj)Var(xi)Var(xj)=ρ^−μ2μ(1−μ)). For a fixed correlation, there was generally a one-to-one relationship between firing rate and D_KL(P, P˜). For these distributions (Figure 7B, for Gaussian inputs; skewed inputs shown in Figure S3B), D_KL(P, P˜) was maximized at an intermediate firing rate. Additionally, D_KL(P, P˜) had a non-monotonic relationship with spike correlation: it increased from zero for low values of correlation, obtained a maximum for an intermediate value, and then decreased. These limiting behaviors agree with intuition: a spike pattern that is completely uncorrelated can be described by an independent distribution (a special case of PME model), and one that is perfectly correlated can be completely described via (perfect) pairwise interactions alone.

We next considered circuits in which inputs were drawn from one of two heavy-tailed distributions, the Cauchy distribution and a heavy-tailed skewed distribution, defined earlier. Here, distinctly different patterns emerge: for a fixed Θ, D_KL(P, P˜) is maximized in regions of high input correlation and high input variance σ, but relatively high values of D_KL are achievable across a wide range of input values (see Figure 7C for Cauchy inputs; heavy-tailed skewed in Figure S3C). However, the maximum achievable values of D_KL were achieved at intermediate output correlations ρ ≈ 0.4 (see Figure 7D for Cauchy inputs; heavy-tailed skewed shown in Figure S3D); this suggests that high input correlations do not result in high output correlations.

This somewhat unintuitive finding may be explained by the structure of the PDF of a heavy-tailed common input, which favors (infrequent) large events at the expense of medium-size events. For instance, the probability that a Cauchy input is above a given threshold (P[I_c > Θ > E[I_c]]) is often much smaller than for a Gaussian distribution of the same variance. However, an input can trigger at best one single spiking event regardless of size: therefore a Cauchy common input generates fewer correlated spiking events with larger inputs, while a Gaussian common input triggers correlated spiking events with smaller, but more frequent, input values. As a result, heavy-tailed inputs are unable to explore the full range of output firing statistics: Figure 7D shows that high output correlations only occur at very low firing rates. Overall, D_KL(P, P˜) reaches higher numerical values than for sub-Gaussian inputs, possibly reflecting the higher-order statistics in the input. However, the maximal D_KL(P, P˜) attained still falls far short of exploring the full range of possible values (compare with Figure 1B).

Finally, we examine the behavior of the strain, which quantifies both the magnitude and sign of deviation from the pairwise model (see Ohiorhenuan and Victor, 2010). It has been previously observed that the strain is negative for the DG model (Macke et al., 2011), a condition that has been related to sparsity of the neural code and with which our results agree (data not shown). However, we found that any other choice of input marginal statistics, both positive and negative values are seen; for heavy-tailed common inputs, positive values predominated except at very low firing rates.

Having shown that a wide range of unimodal common inputs produced spike patterns that are well-approximated by PME fits, we next examined bimodal common inputs. Such inputs substantially increased departures from PME fits in the ganglion cell models described above. As in the previous section, we varied c, σ, and Θ so as to explore a full range of firing rates and pairwise correlations.

As a function of input parameter values, D_KL(P, P˜) is maximized for large input correlation and moderate input variance σ² [see Figure 7E, which illustrates D_KL(P, P˜) for a fixed threshold Θ = 1.5]. Figure 7F shows D_KL(P, P˜) values as a function of the firing rate and pairwise correlation elicited by the full range of possible bimodal inputs. We see that D_KL(P, P˜) is maximized at an intermediate (but relatively high: ν ≈ 0.4) firing rate, and for intermediate-to-large correlation values (ρ ≈ 0.6 − 0.8).

We find distinctly different results when we view Δ (Equation 1), for these same simulations, as a function of output spiking statistics (right panels of Figures 7B,D,F). For unimodal, sub-Gaussian distributions (Figure 7B), Δ is very close to 1, with the few exceptions at extreme firing rates. For heavy-tailed and bimodal inputs (Figures 7D,F), Δ may be appreciably far from 1 (as small as 0.5) with the smallest numbers (suggesting a poor fit of the pairwise model) occurring for low correlation ρ. This highlights one interesting example where these two metrics for judging the quality of the pairwise model, D_KL(P, P˜) and Δ, yield contrasting results.

Finally, we emphasize that while bimodal inputs can produce greater higher-order interactions than unimodal inputs, the values of D_KL(P, P˜) accessible by feedforward circuits with global inputs remain far below their upper bounds at any given firing rate. The maximal values of D_KL(P, P˜) reached by Cauchy and heavy-tailed skewed inputs were 0.0078 and 0.0153; bimodal common inputs reached a maximal value of 0.091. This is an order of magnitude smaller than possible departures among symmetric spike patterns (compare Figure 1B). The difference is illustrated in Figure S4, which compares the D_KL(P, P˜) values obtained in the thresholding model and those obtained by direct exhaustive search at each firing rate by superposing the datapoints on a single axis.

The central finding above is that circuits with bimodal inputs can generate significantly greater higher-order interactions than circuits with unimodal inputs. To probe this further, we investigated the behavior of D_KL(P, P˜) for the feedforward threshold model with a perturbation expansion in the limit of small common input. We found that as the strength of common input signals increased, circuits with bimodal inputs diverged from the PME fit more rapidly than circuits with unimodal inputs; the full calculation is given in the Appendix. In brief, we determined the leading order behavior of D_KL(P, P˜) in the strength c of (weak) common input. D_KL(P, P˜) depended on c³ for unimodal distributions, i.e., the low order terms in c dropped out; for symmetric unimodal distributions, such as a Gaussian, D_KL(P, P˜) grew as c⁴. For bimodal distributions, D_KL(P, P˜) grew as c². Because of the c² dependence, rather than c³ or c⁴, as the strength of common input signals c increases, circuits with bimodal inputs are predicted to produce greater deviations from their PME fits.

We next modified our thresholding model to incorporate the effects of recurrent coupling among the spiking cells. To mimic gap junction coupling in the RGC circuit, we considered all-to-all, excitatory coupling, and assumed that this coupling occurs on a faster timescale compared with the timescale over which inputs arrive at the cells.

Our previous model was extended as follows: if the inputs arriving at each cell elicited any spikes, there was a second stage at which the input to each neuron receiving a connection from a spiking cell was increased by an amount g. This represented a rapid depolarizing current, assumed for simplicity to add linearly to the input currents. If the second stage resulted in additional spikes, the process was repeated: recipient cells received an additional current g, and their summed inputs were again thresholded. The sequence terminated when no new spikes occurred on a given stage; e.g., for N = 3, there were a maximum of three stages. The spike pattern recorded on a given trial was the total number of spikes generated across all stages.

We then explored the impact of varying g for a single representative value of σ and Θ, and several values of the correlation coefficient c. We found that as g increased D_KL(P, P˜) varied smoothly, reflecting the underlying changes in the spike count distribution. For small c (c = 0.02 shown in Figure 8A), where the variance of common input is very small, the results varied little by input type: for all input types D_KL(P, P˜) reached an interior maximum near g ≈ 1.7. As c increases, the distinctions between inputs types become apparent (Figures 8B,C show c = 0.2, 0.5, respectively): for most input types and values of c, the value of D_KL(P, P˜) reaches an interior maximum that exceeds its value without coupling (i.e., g = 0). However, overall values of D_KL(P, P˜) remained modest, never exceeding 0.01 across the values explored here.

We examined a highly idealized model of neural spiking, so as to explore the generality of our earlier findings in a small array of RGC models. We found that our main results from the RGC model—that higher-order interactions were most significant when inputs had bimodal structure, and that when fast excitatory recurrence was added to the circuit, higher-order interactions were maximized at an intermediate value of the recurrence strength—persisted in this simplified model. Moreover, we were able to show that the first of these findings is general, in that it holds over a complete exploration of parameter space.

To model the relationship between the light stimulus and synaptic conductances, the current measurements I^exc and I^inh were fit to a linear-nonlinear model:

where s is the stimulus, L^exc (L^inh) is a linear filter, N^exc (N^inh) is a non-linear function, and η^exc (η^inh) is a noise term. The linear filter was fit by the function

and the non-linear filter by the polynomial

Fits minimized the mean-square distance between model and data. L^inh and N^inh were fit using the same parametrization.

The noise terms η^exc_k, η^inh_k were fit to reproduce the statistical characteristics of the residuals from this fitting. We simulated the noise terms η^exc and η^inh using Ornstein–Uhlenbeck processes with the appropriate parameters; these were entirely characterized by the mean, standard deviation, and time constant of autocorrelation τ_η,exc (τ_{η, inh}), as well as pairwise correlation coefficients for noise terms entering neighboring cells. The noise correlation coefficients were estimated from the dual recordings of Trong and Rieke (2008).

Linear filter parameters computed (also listed in Table 1) were P_exc = −8 × 10⁴ s⁻¹, n_exc = 3.6, τ_exc = 12 ms, T_exc = 105 ms, and P_inh = −1.8 × 10⁵ s⁻¹, n_inh = 3.0, τ_inh = 16 ms, T_inh = 120 ms. Non-linearity parameters were A_exc = −8.3 × 10⁻⁷ nS, B_exc = 7 × 10⁻³ nS, C_exc = −0.95 nS, and A_inh = 1.67 × 10⁻⁶ nS, B_inh = 6.2 × 10⁻³ nS, C_inh = 4.17 nS. Noise parameters were measured to be mean(η^exc_k) = 30, std(η^exc_k) = 500, τ_η,exc = 22 ms, and mean(η^inh_k) = −1200, std (η^inh_k) = 780, τ_η,inh = 33 ms. In addition, excitatory (inhibitory) noise to different cells η^exc_k, η^exc_j (η^inh_k, η^inh_j) had a correlation coefficient of 0.3 (0.15).

For the filter demonstrated in Figure 4, we added a cosine component to the previous filter, i.e.,

Here P_exc,M = −3.2 × 10⁵ s⁻¹, n_exc,M = 2, τ_exc,M = 12 ms, T_exc,M,S = 120 ms and T_exc,M,C = 100 ms, and P_inh,M = −3.5 × 10⁵ s⁻¹, n_inh,M = 2, τ_inh,M = 13.2 ms, T_inh,M,S = 132 ms and T_inh,M,C = 110 ms, while R_exc,M = R_inh,M = 0.8.

We create a model of the cell as a non-linear integrate-and-fire model using the method of Badel et al. (2007), in which the membrane voltage is assumed to respond as

where C is the cell capacitance, t_last is the time of the last spike before time t, and I_input(t) is a time-dependent input current. We use the current-clamp data, which yields cell voltage in response to the input current I_input(t) = −g^exc(t)(V − V_E) − g^inh(V − V_I), to fit a function F(V, t). When voltage data is segregated according to the time since the last spike t − t_last, the I − V curve is well fit by a function of the form

where parameters are the membrane time constant τ_m, resting potential (E_L), spike width Δ_T and knee of the exponential curve V_T.

The values of these constants differed in each bin of voltage data; to estimate these constants, we first extracted their values from each mean I − V curve. We found that these constants, as a function of t − t_last, were well fit by either a single exponential or a difference of two exponentials, with relaxation to a baseline rate (as in Badel et al., 2007, Figure 3). Specifically, we chose:

We obtained the coefficients by least-squares fitting to the above functional forms: specifically, we found that (up to four digits): (c_τm,1, c_τm,2, c_τm,3) = (0.3719 ms⁻¹, 0.5412 ms⁻¹, 13.2726 ms), (c_EL,1, c_EL,2, c_EL,3, c_EL,4) = (−59.4858 mV, 5.8966 mV, 8.3076 ms, 233.1114 ms), (c_ΔT,1, c_ΔT,2, c_ΔT,3, c_ΔT,4) = (20.0487 ms, 19.0560 ms, 3.6280 ms, 2.4304 s), and (c_VT,1, c_VT,2, c_VT,3) = (−44.3323 mV, 25.1812 mV, 4.7653 ms). Coefficients are also listed in Table 2.

The capacitance was inferred from the voltage trace data by finding, at a voltage value where the voltage/membrane current relationship is approximately Ohmic, the value of C that minimizes error in the relation Equation (10) (Badel et al., 2007). The estimated value was C = 28 pF.

For simulations without electronic coupling, our model neuron comprises Equations (10, 11) for V < V_threshold; a spike was detected when V reached V_threshold = −30 mV; voltage was then reset to V_reset = −55 mV. The cell was then unable to spike for an absolute refractory period of τ_abs = 3 ms.

All simulations presented here were done in a three-cell network.

Gap junction coupling was introduced as an additional current on the right-hand side of Equation (10):

The coupling strength g^gap was held constant during a simulation. When coupling was present (i.e., when g^gap ≠ 0), g^gap was varied from the measured level (1.1 nS) (Trong and Rieke, 2008) to 16 times this value (17.6 nS) between simulations. When present, coupling was all-to-all.

As in the feedforward model, Equations (10, 11) were integrated for V < V_threshold, and a spike was detected when V reached V_threshold = −30 mV. To model the voltage trajectory immediately following a spike, an averaged spike waveform was extracted from voltage traces of the same ON parasol cell used to fit Equations (10, 11). This spike waveform was then used to replace 1 ms of the membrane voltage trajectory during and after a spike; at the end of the 1 ms, the voltage was released at approximately −58 mV. The cell was unable to spike for an absolute refractory period of τ_abs = 3 ms. A relative refractory period was induced by introducing a declining threshold for the period of 3–6 ms following a spike, after which V_threshold returns to −30 mV.

We defined each cell's stimulus as the linear convolution of an image with its receptive field. The receptive fields include an ON center and an OFF surround, as in Chichilnisky and Kalmar (2002):

where the parameters k and 1/r give the relative strength and size of the surround. Q specifies the shape of the center and was chosen to have a 1 standard deviation (SD) radius of 50 μm and to be perfectly circular. The receptive field locations x→₁, x→₂, and x→₃ were chosen so that the 1 SD outlines of the receptive field centers will tile the plane (i.e., they just touch). Other parameters used were k = 0.3, r = 0.675.

Stimulation images were defined on a 512 μm × 512 μm grid that overlapped all three receptive fields. For full-field stimuli, light intensity was chosen be spatially constant and refreshed every 8, 40, or 100 ms by choosing independently from the specified stimulus distribution (Gaussian, binary, Cauchy, or heavy-tailed skew). For spatially variable stimuli, a checkerboard pattern was imposed on the stimulation image: the intensity value in each checkerboard square was chosen independently and refreshed at the appropriate interval. The checkerboard pattern was first given a random rotation and translation relative to the receptive fields: this was chosen at the outset of each batch of stixel simulations (for a total of five rotation/translation pairs per stixel size, refresh rate, and stimulus distribution). Two example placements are shown in Figures S2A,D for 256 μm and 60 μm pixels respectively.

All simulations and data analysis were performed using MATLAB. Equations (10, 11) were integrated using the Euler method for >10⁵ ms with a time step of 0.1 ms. The synaptic noise terms, η^exc_k and η^inh_k, as well as the light input, were generated independently for each simulation. In response to uniform light stimuli, firing rates were 11.51 ± 0.38 Hz (standard deviations given across a total of 60 cells; 3 cells each from 20 10⁵ ms simulations); 10 ms bins were used to discretize the spiking output. Firing rates were higher for full-field stimuli, ranging from 12 to 43 Hz (firing rates increased with stimulus variance); therefore shorter (5 ms) bins were used to discretize spike output for all other simulations. With this range of firing rates and bin size, multiple spikes were very rare (occurring in <1% of occupied bins). Empirical spiking distributions were computed from the binned spike data.

For each stimulus condition, 20 simulations (or sub-simulations) were run, for a total integration time of > 20 × 10⁵ ms. These 20 sub-simulations were used to estimate standard errors in both the probability distribution over spiking events and D_KL(P, P˜). Numbers reported in section 2 are, unless specified otherwise, produced by collating the data from the 20 simulations.

To fit a maximum entropy model P˜ to an empirical probability distribution P, we used standard methods that have been explained elsewhere (Malouf, 2002). Briefly, we minimized the negative log-likelihood function:

where

Z_λ is the partition function, f_k, k = 1, …, M is a set of functions or “features” of the spiking state, and λ is a vector of parameters, each of which serves as a Lagrange multiplier enforcing the constraint E_P˜[f_k]. For the pairwise (PME) model on N cells, λ corresponds to N firing rates and N(N − 1)/2 covariances, and the sum is over all possible spiking states of the system. For N = 3 there are six such parameters, and

The function in Equation (15) is a convex function of the parameters λ which will be minimized precisely (and uniquely) when P˜ matches the desired moments from P: e.g., E_P[x₁] = E_P˜[x₁]. Since P˜ is in log-linear form, the result will be the maximum entropy distribution that matches the desired moments (Malouf, 2002). In principle any unconstrained gradient descent method may be used; we used an implementation of the non-linear conjugate gradient method. The Kullback Leibler divergence D_KL(P, P˜) was computed using the identity D_KL(P, P˜) = S(P˜) − S(P), where S(P) is the entropy of P, i.e., S(P) = −∑_x P(x) log P(x).

To test our finding that the observed distributions were well-modeled by the PME fit, we also performed the PME analysis on each of the 20 simulations for each stimulus condition. While in general D_KL(P, P˜) can be quite sensitive to perturbations in P, the numbers remained small under this analysis. To confirm that our results for D_KL(P, P˜) are sufficiently resolved to remove bias from sampling, we performed an analysis in which we collect the 20 simulations in subgroups of 1, 2, 4, 5, 10, and 20, and plot the mean D_KL with estimated standard errors. As expected (e.g., Paninski, 2003), bias decreases as the length of subgroup increases and asymptotes at—or before—the full simulation length.

To provide a cross-validation test for the significance of our reported D_KL(P, P˜) values, we divided our data into halves (which we denote P₁ and P₂, each including data from 10 sub-simulations) and performed the PME analysis on one half (say P₁) to yield a model P˜₁. We then computed D_KL(P₂, P˜₁) and D_KL(P₂, P₁) (as in Yu et al., 2011), which we refer to the cross-validated and empirical likelihood, respectively. The former tests whether the PME fit is robust to over-fitting; the latter tests how well-resolved our “true” distribution is in the first place. Most cross-validated likelihoods fall on or near the identity line; most empirical likelihoods are close to zero [and importantly, significantly smaller than either D_KL(P, P˜) or D_KL(P₂, P˜₁), indicating that D_KL(P, P˜) is accurately resolved]. We conclude that the deviations that we observe when these conditions are met can not be accounted for by the differences in testing and training data.