Near-Optimal Decoding of Transient Stimuli from Coupled Neuronal Subpopulations

Model of the VS network.

Our study is based on a model of the VS tangential cells closely related to that of Borst and Weber (2011). We briefly describe the model, and note the differences between the specific implementations. Parameters not explicitly stated, and details of the model not discussed are identical to those given by Borst and Weber (2011).

To mimic stimulation of the fly visual system under a variety of conditions, including natural flight, we started by projecting a random image onto the surface of a sphere. We considered images that consist of randomly arranged bars of varying sizes, as well as random checkerboard images and compositions of natural scenes. Spherical images were rotated about a horizontal axis with varying azimuthal angle, thereby generating a pattern of optic flow (Fig. 2A). Details on image and optic flow stimulus generation are given in the next section.

Processing of these stimuli by the fly visual system is captured by several successive computational steps. The rotated image sequences (“optic flow stimuli”) were first filtered by an array of vertically oriented local motion detectors (LMDs or “Reichardt detectors”; Reichardt, 1987; Borst et al., 2003; Haag et al., 2004). The LMDs were spaced approximately evenly on the surface of the sphere. There were 5000 detectors per hemisphere, corresponding approximately to the number of facets on the left and right eyes (Hengstenberg, 1992). The input to a single detector was composed of luminance signals from two vertically aligned pixels separated by an elevation of 2°. First-order filtered low- and high-pass versions of the input from the two pixels were cross-multiplied and then subtracted (Fig. 2B). A negative (respectively positive) detector response reflected upward (respectively downward) motion. The downward and upward components were separately weighted according to the dendritic receptive field (RF) of each cell, and activated the dendritic compartment of the VS model neurons as excitatory and inhibitory conductances, respectively. Dendritic RFs were vertically centered Gaussians, with horizontal width 15° and vertical width 60° (Fig. 2C). Hence, each cell effectively sampled the entire vertical surround above and below its RF center, given in Figure 1A. Maximal vertical velocity within a cell's RF was generated by azimuthal rotation angles approximately orthogonal to the centers of the RF. Such rotations therefore resulted in maximal excitation or inhibition.

The axonal compartments of adjacent, ipsilateral VS neurons are electrically coupled to each other. Figure 2D shows a schematic of the model processing stages. Figure 2E shows the response of each VS neuron to downward stimulation in a narrow (10° wide) vertical strip which was swept across the visual field. The vertical strip contained a square-wave horizontal grating with a spatial frequency of 22.5° drifting downward at a constant velocity of 125°/s, corresponding to a temporal frequency of 5.5 Hz.

The axonal and dendritic membrane potentials for the VS neurons in each hemisphere evolve according to the following: Here V_Ax and V_De are vectors whose entries correspond to the 10 axonal and dendritic voltages, respectively. The full VS model consists of two copies of this system, representing the activity of the system in the left and right hemispheres. The two differ only in their RF centers (Figs. 1A, 2C). The parameter g_Ax-Den sets the conductance for the coupling of the axonal and dendritic compartments of each neuron, whereas C_m is the membrane capacitance and g_L,Ax, g_L,De are the leak conductances of each compartment. The membrane time constant of each compartment is: τ_X = C_m / g_L,X, X = Ax, De. The resting potential of each compartment is zero. Following a perturbation from rest, the membrane potential decays exponentially back to the resting potential with a characteristic timescale τ_X. Intrinsic variability is modeled by standard white noise processes ξ_Ax(t), ξ_De(t), and σ_Ax, σ_De, set the noise intensities.

The input currents to the dendrite of cell i result in the term I_i(t) = E_Eg_E,i(t), where g_E,i(t) is the excitatory conductance to cell i induced by the optic flow stimulus, E_E is the associated reversal potential, and likewise for the inhibitory quantities g_I,i(t) and E_I. The matrix G_De(t) is diagonal with entries describing the leak conductance, axonal coupling, and input currents as follows: The matrix G_Ax has entries given by: Here, g_gap sets the strength of the axo-axonal gap junction coupling between adjacent, ipsilateral VS neurons. One difference between our simulation protocol and that of Borst and Weber (2011) is that we generated visual inputs at time steps of 1 ms, but integrated Equation 1 at a smaller time step of 0.01 ms to guarantee numerical accuracy. We first calculated the conductances elicited by the optic flow stimulus at the coarser time step, then linearly interpolated to obtain a realization of the conductance at the finer timescale. Typical responses of the uncoupled and coupled systems to rotation of a random bar image at θ_stim = 90° are shown in Figure 3A,B, respectively. Central to the ability of the VS population to encode the axis of rotation is the strong, nonlinear dependence of a VS neuron's response on the rotational velocity of the visual stimulus within its RF.

When we consider the encoding of the rotation axis in the VS axonal responses, we take the output of the system to be temporal averages of the axonal membrane potential. For transient responses, the VS output is as follows: In particular, when considering transient responses, we assume the system starts from rest (0 mV) at the beginning of the period over which we average. Similarly, steady-state responses were computed as follows: In contrast to the transient response defined above, the steady-state response is defined so that, at the beginning of the integration period (τ_ss), the entire VS system is (approximately) in steady-state. Despite the fast time constants of VS model neurons, they do not immediately reach steady-state, because it takes some time for the motion detector-filtered stimulus received by the VS dendrites to equilibrate. In Figure 3, the shaded boxes indicate the periods over which we calculated the transient and steady-state responses.

It has been observed that there is a mutually inhibitory interaction between the end cells (VS1 and VS10) in each hemisphere. This may be implemented by electrical coupling of VS7–10 to an inhibitory cell Vi which forms a chemical synapse onto the ipsilateral VS1 cell, and electrical coupling of VS1 to an inhibitory cell Vi2 which forms a chemical synapse onto the ipsilateral VS10 cell (Haag and Borst, 2007; Borst and Weber, 2011). Following Weber et al. (2008), we implemented this “repulsive” coupling using a negative-conductance gap-junction between VS1 and VS10 (g_inh in Eq. 2). This repulsive coupling was scaled when we changed the strength of axo-axonal gap junction coupling between VS neurons. Unless otherwise specified, we set g_inh = −0.06 g_gap. Our findings do not depend qualitatively on the presence of this connection (results not shown).

For simplicity, we did not model several known functional and anatomical properties of VS cells, such as the rotational structure of their RFs (Krapp and Hengstenberg, 1996; Krapp et al., 1998) or dendrodendritic connections with the dCH neuron (Haag and Borst, 2007). These properties of the VS network are the subject of current experimental investigations (Hopp et al., 2014). Although we do not expect them to significantly affect our conclusions, investigating the impact of such VS cell features is an important avenue for future investigation. Previous computational studies of the VS network have made similar simplifying assumptions (Karmeier et al., 2005; Cuntz et al., 2007; Weber et al., 2008; Elyada et al., 2009, 2013).

Generation of images and optic flow patterns.

Optic flow patterns were generated by first projecting various types of random images onto the surface of the unit sphere. We considered three classes of random images: random bars, random checkerboards, and natural scenes (Fig. 1C shows examples of each type of image). In the first two cases, images were binary, pixel intensities were either 0 or 1, and for natural scenes, pixel intensities varied continuously between 0 and 1. All images were discretized at 1° increments in spherical coordinates. Throughout, images were generated independently across trials.

Random bar images were parameterized by the number of bars, as well as bar width and length. Each bar was generated by first randomly placing an initial line segment of the specified bar width on the surface of the unit sphere. We then expanded this segment along the direction of the length of the bar by rotation about the appropriate axis, turning “on” all pixels along the path touched by the rotating segment. Bar images used in Results (see Figs. 5⇓–7) contained 25 bars of length 40° and width 5°. Choosing bars with different dimensions or changing the number of bars did not affect the results qualitatively (data not shown).

For checkerboard images, we defined a coarse discretization of the image consisting of 4° × 4° squares, and randomly set all pixels within a square to be zero or one, independently across squares. Last, for natural images, we first took a subcollection of one hundred natural scenes from the van Hateren and Schilstra (1999) dataset. These images were selected to exclude man-made objects, sharp edges, and gratings. We then randomly selected six (with replacement) of these 100 images, projected them onto the sides of a cube, which was itself then projected onto a sphere, mimicking the approach of Borst and Weber (2011). For natural scene compositions, we also included an initial rotation of random magnitude about a randomly chosen (not generally horizontal) axis to control for the effects of edges between different natural scenes. We also confirmed our results with images chosen randomly from the van Hateren and Schilstra (1999) dataset without restrictions. Results were nearly identical (data not shown).

The sequences of images comprising optic flow patterns were generated by rotation of the sphere about an axis in the horizontal plane (we did not consider translatory motion). At a given positive time, the value of a pixel was set equal to the value of the pixel obtained by a reversal of the rotation applied to the original image. The rotational velocity was constant across simulations and set to 500°/s, falling well within the parameters of typical motion of the fly during flight (Egelhaaf et al., 2012). This value is also consistent with values considered in previous computational studies of the VS network (Karmeier et al., 2005; Cuntz et al., 2007; Weber et al., 2008; Elyada et al., 2009). Increasing or decreasing the rotational velocity to 250°/s or 750°/s did not affect the results quantitatively, nor change our general conclusions (data not shown).

Optimal linear and zero-crossing estimators.

To assess the ability of the VS cell network to encode the direction of rotation, we considered several estimators based on the axonal membrane potential of VS cells. In the following three sections, the reader should keep in mind that the random vector V̄ is a surrogate for the time-averaged axonal response of VS cells. We will define the optimal linear estimator of the rotation axis based on the steady-state averaged axonal response of VS cells, V̄_Ax^ss, defined in Equation 4. In subsequent sections, we will also define the minimum mean-square estimator (MMSE) using the transient averaged axonal response, V̄_Ax^tr, defined in Equation 3.

For our analysis of the steady-state encoding of the axis of rotation, we applied an optimal linear estimator (OLE), the linear estimator which minimizes the expected value of the squared error over all stimuli and responses (Salinas and Abbott, 1994). We considered a linear, rather than an affine estimator because of the (near) rotational symmetry of the system. The OLE is simple and intuitive, and we used it to demonstrate the impact of gap junction coupling between VS cells. For a more detailed analysis we used the MMSE defined in the section entitled “Computation of the MMSE.”

We are interested in estimating the axis of rotation which is characterized by the unit vector, s, pointing along its direction. This vector is in the horizontal plane of the fly, with azimuthal angle θ_stim. The OLE is a linear combination of the responses of the N neurons, V̄ = (V̄₁,…,V̄_N)^T (where T denotes transposition). To obtain the OLE, we denote the joint probability density of stimuli and responses by P(V̄,s) = P(V̄|s)P(s). We assume throughout the following that the prior distribution over the stimuli, P(s), is flat. The tuning curve for the i′th neuron is then: We also set L_i = E[sμ_i(s)] = ∫sμ_i(s)P(s) ds, and let L = [L₁|…|L_N] be the matrix with columns L₁,…,L_N. We denote the second moment of the responses of the i′th and j′th neurons (averaged across stimulus values) by Σ_ij = ∫ V̄_iV̄_jP(V̄|s)P(s)dsdV̄. Given an observed response, V̄_obs, the OLE then has the following form: where Σ is the matrix of second moments, Σ_ij. Note that the OLE requires measurement of only first and second moments of the response. It is thus considerably simpler to obtain than the true MMSE defined in the next section, which requires knowledge of the full distribution, P(V̄ | s). The optimal linear estimator and the MMSE coincide only when the joint distribution P(V̄, s) is Gaussian, which is not generally true for the system we consider.

The stimulus s estimated by the vector ŝ was a unit vector pointing along the axis of rotation in the horizontal plane. However, the azimuthal angle, θ_stim, of s is the behaviorally relevant variable for optomotor control by the fly. We therefore report the angle θ̂ that the vector ŝ makes with the direction that the fly is facing, i.e., the azimuthal angle of ŝ, as in earlier analyses of similar directional stimuli (Georgopoulos et al., 1988; Salinas and Abbott, 1994; Lewis and Kristan, 1998).

We also implemented an estimator proposed by Cuntz et al. (2007) and Elyada et al. (2009) who suggested that the axis of rotation can be defined as a “zero crossing” of the population response. To obtain this estimator we defined for each cell a “zero angle”. Intuitively, this is the azimuthal rotation angle that is most likely to elicit no response for a given cell. Because rotation about the azimuthal angle coinciding with the center of a VS cell RF will elicit little downward and upward motion within the RF, the zero angles coincide precisely with the centers of the RFs shown in Figure 1A. Then, to estimate the “zero crossing” rotation angle given the response to a rotation about the stimulus angle, θ_stim, we search for consecutive pairs of VS neurons which exhibit a sign change in their responses (i.e., their membrane potential responses lie below and above their resting values, respectively). These neurons have the smallest responses to the stimulus, and hence the rotation angle is likely to lay between their respective zero angles. The axis of rotation is then estimated by linearly interpolating between the two zero angles based on the responses of these two VS neurons.

More precisely, we define the zero angle of the i′th cell on each side, θ_i⁰, as the angle that maximizes the likelihood of giving a zero axonal response: The posterior distribution of the azimuthal angle of the stimulus conditioned on the zero response of a VS cell will generally have two relative maxima, about 180° apart. To resolve this ambiguity we choose the angle at which voltage is increasing with increasing θ.

To implement the zero-crossing estimator, we first search for a pair of consecutive VS neurons which exhibit responses V̄₁ < 0, V̄₂ > 0, having zero angles θ₁⁰, θ₂⁰, where we assume the two cells are labeled so that θ₂⁰ − θ₁⁰ mod 360° ∈ [0°, 180°). A sufficient condition for such a pair to exist is existence of at least one pair of VS neurons exhibiting positive and negative responses, respectively. Once such a pair is located, the zero-crossing estimate θ_est^ZC is the angle associated with the zero potential level for the line connecting (θ₁⁰,V̄₁) and (θ₂⁰,V̄₂),

Modeling the joint distribution of VS axonal responses using copulas.

For our analysis of the encoding of the axis of rotation in the transient state, we applied the true MMSE. Obtaining this estimator of the rotation axis requires estimating the joint probability distribution of axonal membrane potential of VS cells given the stimulus. However, even with the benefit of modern computational power, it is not feasible to directly estimate probability distributions for continuous variates in more than a few dimensions. For this reason, we must first formulate an approximation of the joint probability distribution of VS axonal responses to implement the MMSE.

Two approaches to this problem are to fit a maximum entropy distribution that matches a set of empirical statistics of the data (Jaynes, 1957; Roudi et al., 2009; Shlens et al., 2009; Ohiorhenuan et al., 2010; Fairhall et al., 2012), or to apply copulas (Nelsen, 2006). We chose the latter approach, common in the valuation of financial derivatives, but not widely applied in neuroscience (cf. Berkes et al., 2009; Onken et al., 2009a,b). One advantage of the copula approach is that it allows us to use the empirical marginal probability distributions in the fit.

To fix ideas, we remind the reader that we will fit a copula to the probability distribution of the time-averaged transient VS axonal response, defined in Equation 3. Thus, we consider a random vector V̄ = (V̄₁,…,V̄_N)^T with cumulative distribution function F(v̄₁,…,v̄_N). A copula for the distribution function F is a function C : ℝ^N→ℝ such that: Here F_i(·) is the marginal cumulative probability distribution function for the variable V̄_i, i = 1, …, N. The copula C exists for any distribution F with marginals {F_i}_i=1^N (cf. Nelsen, 2006, his Theorem 2.10.9). From Equation 5, it is clear that C determines completely the intervariable dependence structure contained in the distribution F in terms of the marginal distributions, F_i.

An N-dimensional copula is equivalent to a distribution function on the N-dimensional unit hypercube [0,1]^N with uniform marginals: define the random vector U = (U₁, …, U_N) where U_i = F_i(V̄_i). The probability integral transform implies that each U_i is a marginally uniformly distributed random variable (Gabbiani and Cox, 2010, their Sect. 11.8). The copula C(u₁, …, u_N) for V̄ has an equivalent definition as the distribution function of U.

The “curse of dimensionality” prevents us from directly approximating the corresponding copula C. A common approach is to select a parameterized copula family which can then be fit via the maximum likelihood principle (Yan, 2007). We applied the Gaussian copula (Xue-Kun Song, 2000), which takes the following form: Here Φ_Σ is the joint Gaussian distribution function with correlation matrix and Σ; (i.e., Σ_ii = 1 for each i) and Φ is the standard univariate Gaussian distribution function. Given independent identically distributed samples of a random vector V̄ = (V̄₁, …, V̄_N)^T with marginals [F_i], the correlation Σ_ij equals (Bouyé et al., 2000): Here corr(x,y) denotes the correlation coefficient of x and y (Xue-Kun Song, 2000).

Given a general copula distribution function C as defined in Equation 5, the copula density is defined as follows: c(u₁,…,u_N) = ∂∂u1…∂∂uNC(u1,…uN).

The joint density f corresponding to the distribution F may be written as follows: where each f_i is the marginal density corresponding to the distribution F_i. The density of the Gaussian copula may be expressed in closed form as follows: We fit a Gaussian copula to the transient response of the system, V̄_Ax^tr(T). To test the goodness of the fit distribution, we selected 20 random subsets of three left-side VS neurons, and associated with each subset a random stimulus rotation angle. We then sampled the marginal copula for each subset [i.e., C(u_i, u_j, u_k) with u_i = F_i(V̄_i) in Eq. 5], comparing it with the fit copula used in the MMSE calculations using Equations 6 and 7. We compared the empirically observed and fit copula distribution values at 1000 equally spaced points in the unit cube of the form (0.1i, 0.1j, 0.1k), 1 ≤ i, j, k ≤ 10. In Figure 4A, we present a probability-probability (P–P) plot of the true (empirical) copula values against the fit values. In other words, letting C^true and C^fit indicate the true and fit copula distribution functions, Figure 4 presents a scatter plot of the 1000 points: We also computed for each of the 20 subsets and for all 1000 sample points the relative error between the probabilities from the true and fit copula distributions. In particular, we defined the following for each point: Because the probabilities lie within [0,1] by definition, the relative errors also lie within [0,1], with the value 0 indicating a perfect match. In Figure 4B, we plot a histogram of the relative errors for all 20 random subsets (thus comprising a total of 20,000 data points). We found the true and fit copula distributions generally agreed quite well. The average relative error across all 20,000 points was ≈0.0438, and 90.3% of relative errors were <0.1.

Computation of the MMSE.

To obtain the MMSE (Kay, 1993), we first simulated the response of the VS network to determine an empirical estimate of the marginal distribution functions F_i. We did not assume a parametric form for the marginal distributions, but obtained a discrete estimate by binning values of the membrane potential integral at a sufficiently fine resolution. We then fit the Gaussian copula to the joint responses. Additional details are given in the last section of Materials and Methods.

Both the marginal distributions and the copula were determined as a function of the stimulus rotation angle θ_stim at a resolution of 1°. Marginal distribution histograms were approximated from 10,000 samples at each rotation angle, and the copula from 1000 samples. The MMSE was then computed based on 1600 samples taken at 5° increments.

The MMSE of the axis of rotation is calculated as the mean of the posterior distribution of the rotation vector given the axonal membrane potentials. To be precise, given an observed response V̄_Ax^tr(T), we associated each stimulus angle value θ_stim with a corresponding two-dimensional unit rotation vector s(θ_stim) = [cos(θ_stim), sin(θ_stim)]^T. The estimate, θ̂_stim^MMSE, was obtained by first computing an estimate of s by averaging over the posterior distribution, i.e.: and then reporting the azimuthal angle of this estimate, θ̂_stim^MMSE = arg ŝ^MMSE, as done above for the OLE and zero-crossing estimators. Values of the posterior density P(θ_stim|V̄_Ax^tr(T)) were determined using the fit copula and the measured marginal distributions, along with Equations 8 and 9. The integral over the posterior density was calculated via simple Riemann integration at a 1° discretization.

In much of the Results we obtain the MMSE from a partial readout of the VS response. The estimation procedure is the same regardless of the size of the response vector V̄_Ax^tr(T). However, we will sometimes use the notation V̄_Ax^subpop,tr(T) to emphasize that an estimate is based on the response of a subpopulation.

An approximating Ornstein-Uhlenbeck model.

Assessing the generality of our results required determining whether our observations depended on the details of the model. In addition, we wished to isolate the essential features governing the modeling results. To do so, we derived a simplified model that shares the essential characteristics of the full model defined in Equation 1. The first change was to replace the time-dependent parameter G_De(t) by a constant, so that the simplified model became an Ornstein–Uhlenbeck (OU) process. Furthermore, we replaced the optic flow generated input by spatially correlated noise. However, cells retained the sinusoidal tuning curves of the full model. Finally, the correlations between the inputs to different cells decayed with the distance between them to capture the effect of overlapping RFs. In this simplified model, described by Equation 12, we again changed the magnitude of the diffusive coupling between the cells to examine its impact on encoding.

In form, the corresponding Langevin equations for the evolution of the coupled OU processes X_Ax, X_De are similar to those of the full model (Eq. 1) as follows: The matrix A_Ax has entries given by Thus, X_Ax, X_De are 10-dimensional processes, with a copy of the system for each hemisphere. The two copies are uncoupled, and differ parametrically only in their RF centers, as in the full model. We labeled the parameters so as to facilitate comparisons with their counterparts in the full model. For example, the parameter a_gap captures the coupling between neighboring compartments in the OU model, and corresponds to the parameter g_gap in the full model.

The tuning curves, μ(θ) in Equation 12, describe the steady-state response in the absence of fluctuations. To approximate the steady-state responses of the full system, we defined individual tuning curves as sinusoids, μ_i(θ) = sin(θ − ϕ_i), where ϕ_i indicates the RF center of the i′th cell. Therefore, a cell will respond most strongly to stimulation angles orthogonal to the cell's RF center, as in the detailed VS model (Fig. 2E). We again obtained estimators from the time integrals of the cell's responses. In Equation 12, the last term in the equation for X_De modeled the correlated input to the dendrites. The matrix B_De was chosen so that correlations to neighboring dendritic compartments decayed exponentially with a space constant of 10°: if Δφ was the shorter angular distance between the RFs of two neurons, the correlations between the inputs to the two was exp( −Δφ/10°).

In Figure 12, we used the following parameters: τ_Ax,OU = 0.2, τ_De,OU = 40, a_Ax-Den = 1, σ_Ax,OU = σ_De,OU = 0.2. RF centers, ϕ_i, agreed with the full model. The window of integration for responses was T = 10 units of dimensionless time. The coupling strength a_gap varied, and values used are indicated in the legend of Figure 12.

Generation of figures.

In all figures, to generate a single sample from the true distribution of the VS model axonal response we first used Equation 1 to model the response of the VS model to random optic flow stimuli. We then computed for each individual simulation the corresponding temporal average (Eqs. 3, 4) to obtain a single sample from the response distribution.

To compute the approximate minimum mean-square estimate used to generate Figure 10 and similar figures, we first had to approximate the joint distribution of V̄_Ax^tr, the transient averaged VS axonal responses. We did so in two steps: we first estimated the marginal distribution for each V̄_Ax,i^tr by generating samples as described above and binning at a sufficiently fine resolution. For checkerboard and bar images, we generally binned over the interval [−8 mV, 8 mV] at a resolution of 0.038 mV (420 bins). For natural scenes, where responses were weaker, we generally binned over the interval [−1 mV, 1 mV] at a resolution of 0.0048 mV (420 bins). We verified that bin sizes were small enough so that the results did not change with a further decrease in bin size (results not shown). We determined the marginal histograms at a 1° resolution, for each different time window considered. The histograms were approximated using 10,000 data points at each rotation angle. Given this large sample of realizations of the random vector V̄_i, the values F_i(V̄_i) were approximated by applying a rank transformation (Berg, 2009).

In the second step, we fit the Gaussian copula, also by generating samples directly from the model. The correlation matrix, which parameterizes the Gaussian copula, was determined using Equation 7. As with the histograms, we fit the copula at each angle at a 1° resolution and for each different time window considered. We fit the copula to 1000 samples generated independently from those used to fit the histograms.

With these approximations of the true histograms and maximum likelihood estimates of the parameters for the fit of the Gaussian copula, we drew 1600 new samples from the true model, at a spatial resolution of 5° for θ_stim. For each sample, we computed the approximate MMSE. By averaging across these samples at each value of θ_stim, we obtained estimates of the relation between θ_stim and the mean-square error of the MMSE, as plotted on the left of each panel in Figure 10 and similar figures. The averaged mean-square error shown in the bar plots on the right of each panel in Figure 10 was obtained by averaging these data over all values of θ_stim. The averaged mean-square error calculated for each subset in Figure 9 was determined likewise.

Given the intensive nature of the simulations involved in obtaining even a single data point in these figures, all computations had to be performed on supercomputer clusters. For instance, the sum of all computations performed to generate a single curve in the left panel of Figure 10A amounted to >2000 CPU hours on 2.2 GHz AMD Barcelona processors. The computer code necessary to implement the VS cell system is deposited in the ModelDB repository (accessible at http://senselab.med.yale.edu/modeldb).