Spike-time reliability of layered neural oscillator networks

Abstract

We study the reliability of layered networks of coupled “type I” neural oscillators in response to fluctuating input signals. Reliability means that a signal elicits essentially identical responses upon repeated presentations, regardless of the network’s initial condition. We study reliability on two distinct scales: neuronal reliability, which concerns the repeatability of spike times of individual neurons embedded within a network, and pooled-response reliability, which concerns the repeatability of total synaptic outputs from a subpopulation of the neurons in a network. We find that neuronal reliability depends strongly both on the overall architecture of a network, such as whether it is arranged into one or two layers, and on the strengths of the synaptic connections. Specifically, for the type of single-neuron dynamics and coupling considered, single-layer networks are found to be very reliable, while two-layer networks lose their reliability with the introduction of even a small amount of feedback. As expected, pooled responses for large enough populations become more reliable, even when individual neurons are not. We also study the effects of noise on reliability, and find that noise that affects all neurons similarly has much greater impact on reliability than noise that affects each neuron differently. Qualitative explanations are proposed for the phenomena observed.

Appendices

Appendix A Review of random dynamical systems theory

In this appendix, we review some relevant mathematical theory that justifies the use of Lyapunov exponents in determining the reliability of a system. “Reliability” here refers exclusively to “neuronal reliability”. As these results are very general and can potentially be used elsewhere, we will present them in a context considerably more general than the system defined by Eq. (1).

Consider a stochastic differential equation (SDE) of the form

$$ \label{eq:sde} dx_t = a(x_t)\ dt + \sum\limits_{i=1}^k b_i(x_t)\circ dW^i_t\ . $$

(7)

Here x _t  ∈ M where M is a compact Riemannian manifold of any dimension d ≥ 1, a(·) and b(·) are smooth functions on M, and $(W^1_t, \cdots, W^k_t)$ is a k-dimensional standard Brownian motion. In general, the equation is assumed to be of Stratonovich type, but when $M\,=\,{\mathbb T}^N \equiv {\mathbb{S}}^1 \times {\mathbb{S}}^1 \times \cdots \times {\mathbb{S}}^1$, we have the choice between the Itô and Stratonovich integrals.⁵ Equation (1) is a special case of this setting with $M={\mathbb T}^N$.

1.1 A.1 Stochastic flows associated with SDEs (see e.g. Kunita 1990; Baxendale 1992)

In general, one fixes an initial x ₀, and looks at the distribution of x _t for t > 0. Under fairly general conditions, these distributions converge to the unique stationary measure μ, the density of which is given by the Fokker-Planck equation. For our purposes, however, this is not the most relevant viewpoint. Since reliability is about a system’s reaction to a single realization of Brownian motion at a time, and concerns the simultaneous evolution of all or large ensembles of initial conditions, of relevance to us are not the distributions of x _t but flow-maps of the form $F_{t_1,t_2;\omega}$. Here t ₁ < t ₂ are two points in time, ω is a sample Brownian path, and $F_{t_1,t_2;\omega}(x_{t_1})=x_{t_2}$ where x _t is the solution of (7) corresponding to ω. A well known theorem states that such stochastic flows of diffeomorphisms are well defined if the functions a(x) and b(x) in Eq. (7) are sufficiently smooth (see Kunita 1990). More precisely, the maps $F_{t_1,t_2;\omega}$ are well defined for almost every ω, and they are invertible, smooth transformations with smooth inverses. Moreover, $F_{t_1,t_2;\omega}$ and $F_{t_3,t_4;\omega}$ are independent for t ₁ < t ₂ < t ₃ < t ₄. These results allow us to treat the evolution of systems described by (7) as compositions of random, i.i.d., smooth maps.

Since reliability questions involve one ω at a time, the stationary measure μ, which gives the steady-state distribution averaged over all ω, is not the object of direct interest. Of relevance are the sample measures {μ _ω }, which are the conditional measures of μ given the past. More precisely, we think of ω as defined for all t ∈ ( − ∞ , ∞ ) and not just for t > 0. Then μ _ω describes what one sees at t = 0 given that the system has experienced the input defined by ω for all t < 0. Two useful facts about these sample measures are

• (a) (F _− t,0;ω )_* μ→μ _ω as t → ∞, where (F _− t,0;ω )_* μ is the measure obtained by transporting μ forward by F _− t,0;ω , and

• (b) the family {μ _ω } is invariant in the sense that $(F_{0,t;\omega})_*(\mu_\omega) = \mu_{\sigma_t(\omega)}$ where σ _t (ω) is the time-shift of the sample path ω by t.

Thus in the context of a reliability study, if our initial distribution is given by a probability density ρ and we apply the stimulus corresponding to ω, then the distribution at time t is (F _0,t;ω )_* ρ. For t sufficiently large, one expects in most situations that (F _0,t;ω )_* ρ is very close to (F _0,t;ω )_* μ, which by (a) above is essentially given by $\mu_{\sigma_t(\omega)}$. The time-shift by t of ω is necessary because by definition, μ _ω is the conditional distribution of μ at time 0.

1.2 A.2 Lyapunov exponents of random dynamical systems (see e.g. Arnold 2003)

The fact that the evolution of systems described by (7) can be represented as compositions of random, i.i.d., smooth maps allows us to tap into a large part of dynamical systems theory, namely the theory of random dynamical systems (RDS). Many of the techniques for analyzing smooth deterministic systems have been extended to this random setting, including the notion of Lyapunov exponents. For the stochastic flows above, the largest Lyapunov exponent is defined to be

$$ {\lambda_{\rm max}}(x, \omega) = \lim_{t \to \infty} \frac{1}{t} \log \|DF_{0,t;\omega}(x)\|\ . $$

These numbers are known to be defined for μ-a.e. x ∈ M and a.e. ω. Moreover, they are nonrandom, i.e., they do not depend on ω, and when μ is ergodic, λ _max does not depend on x either, i.e., λ _max is equal to a single (fixed) number for almost every initial condition in the phase space and for almost every sample path.

Numerical calculation of Lyapunov exponents

Lyapunov exponents can be computed numerically by solving the variational equations associated with Eq. (7); the largest Lyapunov λ _max is given by the logarithmic growth rate of a typical tangent vector. This is what we have done in this paper, using the Euler method for SDEs to solve the variational equations.

1.3 A.3 Implications of the sign of λ _max

The two results below are in fact valid in greater generality, but let us restrict ourselves to the SDE setting at the beginning of this subsection.

Theorem 1

Let μ be an ergodic stationary measure of the RDS defined by Eq. (7).

• (1) (Random sinks) (Le Jan 1987) If λ _max  < 0, then for a.e. ω , μ _ω is supported on a finite set of points.

• (2) (Random strange attractors) (Ledrappier and Young 1988) If μ has a density and λ _max  > 0, then for a.e. ω , μ _ω is a random SRB measure.

In Part (1) of the theorem, if in addition to λ _max  < 0, mild conditions (on the relative motions of two points) are assumed, then almost surely μ _ω is supported on a single point (Baxendale 1992). From the discussion above, μ _ω being supported on a single point corresponds to the collapse of trajectories starting from almost all initial conditions to a single trajectory. In the context of Eq.(1), this is exactly what it means for the system to be neuronally reliable as explained in Section 3.1.

The conclusion of Part (2) requires clarification: In deterministic dynamical systems theory, SRB measures are natural invariant measures that describe the asymptotic dynamics of chaotic dissipative systems (in the same way that Liouville measures are the natural invariant measures for Hamiltonian systems). SRB measures are typically singular. They are concentrated on unstable manifolds, which are families of curves, surfaces etc. that wind around in a complicated way in the phase space (Eckmann and Ruelle 1985). Part (2) of Theorem 1 generalizes these ideas to random dynamical systems. Here, random (meaning ω-dependent) SRB measures live on random unstable manifolds, which are complicated families of curves, surfaces, etc. that evolve with time. In particular, in a system with random SRB measures, different initial conditions lead to very different outcomes at time t when acted on by the same stimulus; this is true for all t > 0, however large. In the context of a reliability study, therefore, it is natural to regard the distinctive geometry of random SRB measures as a signature of unreliability.

We do not claim here that mathematically, the results in Theorem 1 apply to Eq. (1). To formally apply these results, conditions of ergodicity, invariant density etc. have to be verified. Evidence—both analytic and numerical—point to an affirmative answer when the coupling constants a, a _i , a _ff etc. are nonzero.

Appendix B A 2-D toy model of two-layer networks

We provide here more detail on how two-neuron models can be used to shed light on two-layer networks as suggested in Section 5.2. Specifically, we will explain how the shapes of the phase distributions P ₁₂ and P ₂₁ are predicted.

Consider a system comprised of two neurons whose dynamics obey Eq. (1). For definiteness, we set ω ₁ = ω ₂ = 1, A _ff = 2.8, and A _fb = 0.8 to mimic the parameters in the two-layer networks considered in Section 5.2, with neurons 1 and 2 representing layers 1 and 2 in the two-layer system. The phase space of this system is the 2-torus, which we identify with the square [0.1]² with periodic boundary conditions; the coordinates are denoted by (θ ₁,θ ₂). In Fig. 12, we show a few trajectories of the undriven system, i.e., with ε = 0. Away from the edges, they are northeasterly with slope 1; near the edges, they are bent due to the coupling. We now turn on the stimulus, setting ε = 2.5 as in Section 5.2. Because only neuron 1 hears the stimulus, it perturbs trajectories only in the θ ₁ direction. When the stimulus is turned on, trajectories will, for the most part, continue to go in roughly the same directions as those shown in Fig. 12, but they become “wriggly”, being driven randomly to the left and right by the white-noise stimulus.

Fig. 12

A few trajectories for a two-oscillator toy model with ε = 0. The trajectories are drawn on a “lift” of the 2-torus to the plane. The parameters are ω ₁ = ω ₂ = 1, A _ff = 2.8, and A _fb = 0.8

We will refer to the top and bottom edges of the square (which are identified with each other) as Σ. Every time a trajectory crosses Σ, neuron 2 spikes, and the location in Σ tells us the phase of neuron 1 when this spiking occurs. We view Σ as a cross-section to the flow, and consider the induced return map Φ: Σ→Σ. In the case of two neurons with feedback, the distribution of trajectories of Φ on Σ tells us the phase distribution of neuron 1 when it receives a synaptic input from neuron 2. In our analogy with the two-layer system, this is the analog of P ₂₁. Similarly, the distribution of returns to the left and right edges identified (we call that Σ′) represents the phases of neuron 2 when it receives an input from neuron 1., i.e., the distribution analogous to P ₁₂.

To understand these distributions, let us view the return dynamics to Σ as the result of two “moves” in succession (this is not entirely accurate but will suffice for present purposes): The first is the return map for the flow with ε = 0, and the second is a “smearing”, obtained by, e.g., taking a convolution with a Gaussian, to simulate the perturbations experienced by the trajectory between returns to Σ.

The return dynamics of the undriven flow are very simple: From the geometry of the flowlines, one sees that starting from any point in Σ, there is a leftward displacement due to the fact that the upward kick along the vertical edges are stronger than the rightward kicks along the horizontal edges (i.e. A _ff > A _fb). This leftward displacement is quite substantial away from θ ₁ = 0, reaching a maximum at θ ₁ ≈ 0.75. Because of the character of the phase response function, this displacement is very small (but strictly positive) near θ ₁ ≈ 0. It is so weak there that with ε = 0, all trajectories spend most of their time near the diagonal, with only brief excursions in between. In other words, when ε = 0, the phase distributions on Σ peak sharply at θ ₁ = 0.

With ε = 2.5, the “smearing” is nontrivial. Immediately, one sees that it causes the distribution to be more spread out. It is also easy to see that some concentration near θ ₁ = 0 will be retained, only that the peak will be more rounded. We now explain why one should expect the peak to be shifted to the right: Suppose we start with a roughly constant distribution on an interval centered at θ ₁ = 0. Since the return map is nearly the identity in this region, we may assume it does not change this distribution substantially. Next we take a convolution, which causes the distribution to have wider support. Now the part of the distribution that is pushed to the right of the original interval will get pushed back in when we apply the return map of the undriven flow again, due to the leftward displacement discussed earlier, whereas the part that is pushed to the left will be pushed further away from θ ₁ = 0 by the ε = 0 dynamics. The result is an obvious asymmetry in the distribution, one that is reinforced in subsequent iterations. (The argument we have presented does not constitute a proof, but a proof is probably not out of reach.)

Once we have convinced ourselves of the geometry of P ₂₁, it is easy to see from the northeasterly direction of the flowlines that if P ₂₁ peaks to the right of θ ₁ = 0 on Σ, then P ₁₂ must peak just below θ ₂ = 0 on Σ′. This completes the explanation promised in Section 5.2.

We finish with a couple of remarks that may be illuminating:

• (1) The phenomenon described above occurs for both excitatory and inhibitory couplings: with A _ff < 0 and |A _fb| < |A _ff|, the displacement of the return map Φ: Σ→Σ is to the right. But in the inhibitory situation, there are other forces shaping the phase distribution, making the picture there more complicated. (In case the reader wonders how to interpret our reasoning in the pure-feedforward case: statements about phase distributions on Σ are valid, except that spikings of neuron 2 do not impact neuron 1!)

• (2) We may also represent the single-layer system by two neurons. Here by far the biggest difference is that both neurons receive the stimulus in the same way, and that translates into perturbations that are in the direction of the diagonal. Such perturbations are not very effective in spreading out distributions, especially when the trajectories are concentrated near the diagonal. These observations provide a geometric understanding for the material in Section 4.2.

A detailed analysis of the two-neuron model with different emphasis is carried out in Lin et al. (2009a).