Diagrammatic expansion for spike train statistics

We will show that any coupled point process model, even one with nonlinearities or negative interactions, has an associated expansion for all spike train cumulants organized by the strength of coupling of higher statistical moments with lower ones (e.g., the influence of the two-point correlation on the mean). The full model we aim to describe is one where each neuron generates a spike train which is conditionally renewal with intensity: (1) Here g_ij(t) is a matrix of interaction filters, λ_i(t) is the baseline point process drive to neuron i and * denotes convolution: (with the integral starting at the initial time for the realization). ϕ_i is the transfer function of neuron i. Neuron j’s spike train is , a sum over Dirac deltas at each of the k spike times. We will take the spike trains to be conditionally Poisson given the input, so that in each time window (t, t + dt), the probability of neuron i generating m spikes is (r_i(t)dt)^m/m! exp (−r_i(t)dt). This corresponds to a point process generalized linear model (GLM), or nonlinear multivariate Hawkes process [32]. In contrast to biophysical or integrate-and-fire models in which spike trains are generated deterministically given the membrane potential (which might, however, depend on noisy input), this model with escape noise generates spike trains stochastically with a rate that depends on the “free” membrane potential (i.e., with action potentials removed) [33].

Current methods for the analysis of single-neuron and joint spiking statistics rely on linear response techniques: using a self-consistent mean field theory to compute mean firing rates and then linearizing the spiking around those rates to determine the stability and correlations. We begin with a simple example highlighting the need to account for nonlinear effects. We take an excitatory-inhibitory network of N_E = 200 excitatory (E) neurons and N_I = 40 inhibitory (I) neurons, all with threshold-quadratic transfer functions . In this example we took network connectivity to be totally random (Erdös-Rényi), with connection probabilities p_EE = 0.2 and p_EI = p_IE = p_II = 0.5. For simplicity, we took the magnitude of all connections of a given type (E − E, etc.) was taken to be the same. Furthermore, the time course of all network interactions is governed by the same filter (with τ = 10 ms), so that g_ij(t) = W_ijg(t). W is a matrix of synaptic weights with units of mV, so that the input to ϕ can be interpreted as the free membrane potential. We set the strength of interactions such that the net inhibitory input weight on to a neuron was, on average, twice that of the net excitatory input weight so that for sufficiently strong interactions, the network was in an inhibitory-stabilized regime [34].

We examined the magnitude and stability of firing rates as we increased the strength of synaptic coupling. We used mean field theory to predict the firing rates, and predicted their linear stability by the standard measure of the spectral radius of the stability matrix . ( denotes the first derivative of neuron i’s transfer function with respect to its input.) As the strength of interactions increases, the mean field prediction for the firing rates loses accuracy (Fig 1A). This occurs well before the mean field theory crosses the stability boundary |Ψ| = 1 (Fig 1B). Examining simulations as the weights are increased reveals an even more fundamental failure of the theory: before the synaptic weights are strong enough for the mean field theory to be unstable, the simulations reveal divergent firing rates (Fig 1C; the raster stops when the instantaneous rates diverge).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Dynamics approaching the firing-rate instability in threshold-quadratic networks.

A) Average firing rate of the excitatory neurons as synaptic weights are scaled. While the ordinate axis shows the excitatory-excitatory synaptic weight, all other weights are scaled with it. Solid lines: prediction of mean field theory. Dots: result of simulation. Inset: threshold-quadratic transfer function. B) Spectral radius of the stability matrix of mean field theory as synaptic weights are scaled. Stars indicate the weight values for the simulations below. C) Example realizations of activity for three different interaction strengths. As synapses become stronger, correlated activity becomes apparent. When synapses are strong enough the activity becomes unstable, even though the mean field theory is stable. All plotted firing rates in A) are averaged over the time period before the rates diverged (if they did). Left: (W_EE, W_EI, W_IE, W_II) = (.025, -.1, .01, -.1) mV. C). Middle: (W_EE, W_EI, W_IE, W_II) = (1, −4, .4, −4) mV. Right: (W_EE, W_EI, W_IE, W_II) = (1.5, −6, .6, −6) mV.

https://doi.org/10.1371/journal.pcbi.1005583.g001

Rather than restricting theoretical analysis to regimes of extremely low coupling strength or linear models, we here develop a framework for activity cumulants that can apply to models with nonlinear input-spike transfer. This will allow us to properly predict both the spiking cumulants and the location of the rate instability in the nonlinear network above. Thus, we develop a framework for activity cumulants that can apply to models with nonlinear input-spike transfer for strongly coupled networks. The mean field and linear response predictions for spiking cumulants correspond to the lowest order terms of this expansion and are good whenever that lowest-order truncation is valid.

We will build this framework systematically: We begin with statistics of the drive λ_i(t), then consider a filtered point process, g * λ(t). In these simple models we will introduce the method and terminology that we will use to understand the more complicated models. We continue by considering the linearly self-exciting Hawkes process, taking a single neuron so g = g and ϕ(x) = x, before proceeding to arbitrary nonlinearities ϕ. Finally, we introduce an arbitrary network structure g. This model is intimately related to the popular GLMs for spiking activity, where the nonlinearity ϕ is commonly taken to be exponential, refractory dynamics can be embedded in the diagonal elements of g, and λ corresponds to the filtered stimulus [26]. A use of GLMs it to fit them to recorded spike trains and then ask about the structure of the inferred functional connectivity g. In contrast, we interpret g as reflecting the structural connectivity and synaptic and membrane dynamics of a specified model network and aim to compute statistics of its neurons’ spike trains. The derivation given here will be heuristic. A more rigorous derivation of the expansion is given in Methods: Path integral representation.

Introduction to the general framework: Poisson process

An inhomogeneous Poisson process generates counts within a window dt independently with a Poisson distribution at rate λ(t). A spike train produced by this process is (2) where t_k is the kth spike time, and N(t) is the spike count. The mean and autocovariance for this process are given by the familiar formulas: (3) (4) where angular brackets denote the expectation across realizations and the subscript c denotes a cumulant, not the moment (i.e., we have subtracted all terms which factor into products of lower moments) [35]. The delta function arises because the process is independent at each time step, so that there is no correlation between events from one time t and any other time t′. In fact, because the events are generated independently at each time point, all of the cumulants of this process can be written as (5) where integrating out one of the delta functions puts the second cumulant in the above form. We can interpret this equation as describing a source of events appearing at rate λ(t) at time t that propagate to times t_i. In this case, because the events are generated independently at each t, the events only propagate to the same time t. For a general point process, cumulants measured at the collection of times {t_i} could be affected by events occurring at any past time, so that we would have to account for how events propagate across time.

The expansion for cumulants we will develop has a natural construction in terms of graphs (Feynman diagrams), wherein components of the graph represent factors in each term. A set of defined rules dictate how terms in the expansion are formed from those graphs. While this graphical representation is not necessary to understand the inhomogeneous Poisson process, we describe it in detail in this simple case to develop intuition and introduce terminology. We use cumulants in this construction because they provide the fundamental building blocks of moments; any n-th order moment can be constructed from up to n-th order cumulants). This also simplifies the developed expansion.

To begin, the nth cumulant has n arguments {t_i, i = 1, …, n}, one for each factor of in Eq 5. We represent each of these by an open white vertex in the graph labeled with the time point, t_i, and we represent the source term λ(t) with a gray vertex. The white vertices are called “external” vertices whereas the gray vertex is called “internal.”

The internal gray vertex represents the intensity of the underlying stochastic process generating events with rate λ(t). The white external vertices represent the spike trains whose statistics we measure at times {t_i}. For each delta function, δ(t − t_i), in Eq 5, we place an edge drawn as a dotted line from the gray vertex to the white vertex, i.e., from the event-generating process to the spike train’s measurement times (Fig 2). More generally, events generated by the source propagate through the graph to affect the external vertices.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Feynman diagrams for the first three cumulants of the inhomogeneous Poisson process.

Each dotted edge corresponds to a delta-function connecting the time indices of its two nodes. White nodes denote the measurement times, while gray nodes denote the times at which spikes are generated. The cumulants are constructed by convolving the term corresponding to the gray node with the product of all outgoing edges’ terms (Eq 5).

https://doi.org/10.1371/journal.pcbi.1005583.g002

In order to construct the term associated with each diagram, we multiply the factors corresponding to edges (delta functions linking t and t_i) or the internal vertex (λ(t)), and then integrate over the time t associated with the internal vertex. This links the generation of events by λ(t) to their joint measurement at times {t_i} through the propagator (here δ(t − t_i)). For the diagrams shown in Fig 2, these rules reproduce the cumulant terms in Eq 5. Note that these graphs are directed, since we only consider causal systems where measured cumulants are only influenced by past events.

In general, a given moment will be the sum of terms associated with many different graphs. For example the second moment is given by (6) Each term on the right-hand side will have a corresponding graph. Moreover, the graph for the second term will include two disconnected components, one for each factor of the mean rate, which appears as in Fig 2. The graphs for the cumulant will always be described by connected graphs.

Filtered Poisson process

We proceed to a simple model of synaptic input: a presynaptic Poisson spike train, with count N(t) and intensity λ(t), drives postsynaptic potentials with shape g(t): (7) where * denotes convolution: (with the integral starting at the initial time for the realization). We assume g is normalized, , so that ϵ gives the magnitude of the filtering.

The cumulants of the postsynaptic potential ν(t) can be calculated directly. In general, they are given by: (8) where the input spikes are generated at times t, arrive at times given by the delta functions and influence the cumulant measured at {t_i} through g. Eq 8 is the same as that for the inhomogeneous Poisson process but with factors of g*. This provides a simple interpretation of Eq 8: cumulants of the filtered Poisson process are given by taking the cumulants of the underlying Poisson process and examining how they can be filtered through the system at hand.

Similarly to the case for the Poisson process, we can represent the cumulants graphically. We again represent each measurement time on the left-hand side of Eq 8 by an external vertex (Fig 3a). The convolution of δ and g in Eq 8 corresponds to an internal time point which we integrate over (denoted by primes). We also represent these internal time points with white vertices that carry a factor of ϵ, the magnitude of the filter. We represent the source term λ(t) with a gray vertex. All vertices that are not “external” are again called internal. Every internal vertex also carries its own time index, t′.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Feynman diagrams for the first three cumulants of the filtered inhomogeneous Poisson process.

A) Cumulant corresponding to the graph. B) Diagrammatic expressions using the filter and the underlying Poisson process. Each dotted edge corresponds to a delta-function connecting the time indices of its two nodes. Each wavy edge corresponds to the filter g connecting the time indices of its two nodes. C) Diagrammatic expressions using the propagator. In all graphs, external white nodes (leaves of the graph) denote measurement times. Gray nodes denote the times at which spikes are generated in the input spike train. Internal white nodes (with time indices t′) denote the times at which input spikes arrive at the postsynaptic neuron. The cumulants are constructed by convolving the term corresponding to the gray node with the product of all outgoing edges’ terms (Eq 8).

https://doi.org/10.1371/journal.pcbi.1005583.g003

The internal gray vertex again represents the intensity of the underlying stochastic process, λ(t). The white external vertices represent the processes whose statistics we measure at times {t_i}. For each delta function, δ(t′ − t), in Eq 8, we place an edge drawn as a dotted line from the gray vertex to the white vertex, i.e., from the event-generating process to the arrival time of the event t′. In this example an event “arrives” as soon as it is generated. A wavy edge corresponds to the filter, g, and represents the effect of a spike arriving at time t′ on the output process measured at time t_i (Fig 3b). Events generated by the source thus propagate through the graph to affect the observed, external vertices.

In order to construct the expression associated with each diagram, we again multiply the factors corresponding to each edge in the diagram (e.g., δ(t′ − t) or g(t_i − t′)) or internal vertex (ϵ or λ(t)) and then integrate over the times associated with the internal vertices. Note that integration over the internal times t′, t′′, etc. results in the convolutions ϵ(g * δ)(t_i − t). Integration over the time t associated with the source term corresponds to the outermost integral in Eq 8 This links the generation of events by λ(t) to their joint measurement at times {t_i} through their arrival times (via δ(t − t′)) and temporal filtering (g(t_i − t′)). For the diagrams shown in Fig 3, these rules reproduce the cumulant terms in Eq 8. Note that the graphs are directed, as for the expansion we describe the “propagator” term will be causal.

We can simplify the cumulants of this process (and the corresponding diagrammatic representations) by considering the propagator of ν(t) (also known as the linear response or impulse response). The propagator measures the change in 〈ν(t)〉 in response to the addition of one input spike in N(t). We can compute it by taking a functional derivative with respect to the input intensity λ(t): (9)

Since the dynamics are linear, this is also equivalent to the change of the expected rate with the addition of one spike to the input, i.e., taking λ(t) ← λ(t) + δ(t′ − t) and 〈ν(t)〉_c ← 〈ν(t)〉_c + Δ * δ(t) (or equivalently the Green’s function of the expected rate). This allows us to rewrite the cumulants in terms of the input rate and the propagator: (10) which can be represented graphically by introducing a solid, directed edge for Δ(t, t′) (Fig 3c). The propagator will be a central feature of the expansion for cumulants in more complicated models involving connections among neurons.

Impact of self-excitation on activity statistics of any order: Linearly self-exciting process

In order to generalize the graphical representation of Poisson cumulants, we begin with a linearly self-exciting process as considered by Hawkes [19]. Let the rate be a linear function of the instantaneous event rate (that is to say the firing rate conditioned on a particular realization of the event history) (11) We assume that g(τ) and λ(t) are such that r(t) > 0, and . If ϵ < 1, then an event will generate less than one event on average, and the rate will not diverge. The history dependence of the firing rate will now enter into our calculations. We can compute the expected rate using the self-consistency equation: (12)

We provide an alternate derivation of this result that will prove useful below: We construct a perturbative expansion of the mean firing rate and show how this expansion can be re-summed to yield the full rate of the self-exciting process. This procedure can also be applied to obtain cumulants of arbitrary order for this process.

We will begin with a recursive formulation of the self-exciting process. In contrast to the filtered Poisson process of the previous section, here the process with count N generates events, which then influence its own rate, dN/dt. Each event can thus generate others in turn. In the case of a linear filter, g, the following approach is equivalent to the Poisson cluster expansion [36–38] and similar to the construction of previous linear response theories for spike train covariances [22]. Define the nth order self-exciting process, N_n(t), to be the inhomogeneous Poisson process given by: (13) where N₀(t) and M_n(t) are inhomogeneous Poisson processes with rates λ(t) and ν_n(t), respectively, where (14) so . M_n(t) is a process with intensity that depends on a stochastic realization of N_n(t), making M₀(t) a “doubly stochastic” process. We can generate these processes recursively: To generate N_n(t), we use a realization of N_n−1(t) to compute the rate ν_n−1 and generate a realization of M_n−1(t). These are added to events generated from the original inhomogeneous Poisson process with rate λ(t) to produce N_n(t). We can use this recursive procedure to develop an expansion for the cumulants of the process at a given order in ϵ (thus a given order in the self-convolution of g).

Let us compute the value of in powers of ϵ using our recursive approach. The zeroth order solution, , is the rate of the inhomogeneous Poisson process λ(t). At order n, we compute using the (n − 1)st order solution in the right-hand side of Eq 13. At first order, using the Poisson solution for we get (15) (16) (17) (18) At second order we similarly arrive at (19) (20) (21) At higher orders we would obtain further terms with additional convolutions with g.

It will be useful to write these expansions in another way, which will allow their form to generalize to nonlinear processes: we will construct the cumulants from the baseline rate and the propagator. We can always replace (22) resulting in (23)

Fig 4a shows the graphical representation of this expansion. As before, the order of the moment is given by the number of external vertices and each external vertex carries a measurement time t_i. We have three types of internal vertices: two open white vertices that carry factors of ϵ (one type has one wavy incoming and one wavy outgoing line; the other has one incoming dotted line and one wavy outgoing line) and one gray vertex (that has one outgoing dotted line). As before, each gray internal vertex corresponds to the source term, and thus represents the factor λ(t). The white internal vertices and their edges represent how the events generated by the source are propagated through the filter g. Each white vertex corresponds to a possible past event time, t′. To construct the cumulant corresponding to a diagram, we integrate over all these possible internal times, weighting each by their influence on the current spiking rate. These weights are given by the filters, g, represented by the wavy edges. The graphical representation of (using the delta function as in Eq 23) is shown in Fig 4a.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Diagrammatic expansion for the mean firing rate and linear response of the self-exciting process.

A) First-order approximation of the firing rate. B) Diagrams corresponding to the re-summing of the expansion of the mean field rate (Eq 25), which is represented by the black dot. C) Diagrams corresponding to the re-summing calculation of the propagator (Eq 28), which is represented by the solid edge. In all diagrams, time indices associated with internal vertices have been suppressed.

https://doi.org/10.1371/journal.pcbi.1005583.g004

We can compute the firing rate of the self-exciting process as the limit of the nth order self-exciting processes, continuing the process outlined for Eq 13: (24) where g⁽ⁿ⁾ is the n-fold convolution of g with itself and g⁽⁰⁾(t) = δ(t). Indeed, we can see that this expression for yields the same recursive self-consistency condition as above: (25)

We can also represent this recursive relation graphically as in Fig 4b, using a black vertex to denote the mean-field rate . The infinite sum defined by Eq 13 has a specific graphical representation: the leftmost vertex and wavy line in the right-hand side of Fig 4b (top) can be detached and factored, with the remaining series of diagrams corresponding exactly to those of the mean. This series of subgraphs on the right-hand side sums to , leading to the recursion relation in Eq 25 (Fig 4b). This graphical representation is equivalent to the recursion relation.

The propagator, Δ(t, t′), measures the fluctuation in the expected rate (around the mean-field value) in response to the addition of one spike at time t′ to the drive λ(t). Setting λ(t) ← λ(t) + δ(t − t′) and in Eq 12 gives: (26) where for convolutions involving Δ(t, t′), we use the notation (f * Δ)(t, t′) = ∫ dt′′ f(t − t′′)Δ(t′′, t′) and (Δ * f)(t, t′) = ∫ dt′′ Δ(t, t′′)f(t′′ − t′)

As with the expected rate , we can examine the propagators of the n-th order self-exciting processes. For the first-order process N₁(t), (27) The first term is the propagator of the inhomogeneous Poisson process. The second term of Δ₁ is the propagator of the filtered Poisson process, Eq 9. This equation can be represented by the same type of graphs as for the expected rate (Fig 4c top), but stand for functions between two time points: the initial time t′ of the perturbation, and the later time t, at which we are computing the rate of the process. We don’t represent these initial and final points as vertices, because the propagator is a function that connects two vertices. However, we still integrate over the times corresponding to the internal vertices since the propagator accounts for the total effect of a perturbation of the source on the observed activity.

In general, the propagator for the nth-order self-exciting process can be computed by taking a functional derivate of the rate with respect to the input rate λ: (28) This recursion relation can be expressed graphically just as for the mean rate (Fig 4c, top). Factoring out ϵg* corresponds to popping off an edge and vertex from the series (Fig 4c, middle). Taking the limit n → ∞ in Eq 28 yields the self-consistency condition for the full propagator Δ(t, t′) given by Eq 26, and indicated by the solid black line in Fig 4c (bottom).

These diagrammatic expansions may seem cumbersome for so simple a model. Even for the self-exciting Hawkes process, however, they allow the fast calculation of any order of spike train cumulant. Let us begin with the second cumulant of the instantaneous firing rate. Again we will construct an expansion in ϵ, i.e., powers of g. To zeroth order, this is the inhomogeneous Poisson solution. To first order in ϵ we have (29) The first term on the second line is the second cumulant of the inhomogenous Poisson process. The other terms arise from the dependency of the processes M₀(t) and N₀(t). The expectation over M₀(t) must be performed first, followed by that over N₀(t), because the process M₀(t) is conditionally dependent on the realization of N₀(t), having intensity (Eq 13). This decomposition relies on the linearity of the expectation operator.

We can construct diagrams for each of these terms using the same rules as before, with the addition of two new internal vertices (Fig 5a). These new vertices are distinguished by their edges. The first has two outgoing dotted lines representing the zeroth-order propagator δ(t − t′), as in the second cumulant of the inhomogeneous Poisson process. It represents events that are generated by the drive λ(t) and propagate jointly to the two measurement time points. The second new vertex has the same two outgoing lines and one incoming wavy line for the filter g(t, t′)–it represents the fourth term on the right-hand side of Eq 29. This vertex carries a factor of ϵ and represents the filtering of past events that then propagate to the two measurement time points.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Diagrammatic expansion for the second cumulant for the self-exciting process.

A) First-order approximation of the second cumulant. B) Re-summing to obtain the full second cumulant. Compare the expansions within the square brackets adjacent to external vertices to the expansion of the propagator, Fig 4c, and compare the expansion of the source term to that of the mean field rate, Fig 4b.

https://doi.org/10.1371/journal.pcbi.1005583.g005

Continuing the computation of the second cumulant to any order in ϵ will result in higher order terms of the linear response and expected rate being added to the corresponding legs of the graph. At a given order n, one leg of each diagram will be associated with a particular term in the expansion, to order n, of the expected rate or the linear response. The second cumulant of dN₂/dt would thus add diagrams with two filter edges to the diagrams of Fig 5a, either both on the same leg of the graph or distributed among the graph’s three legs.

As with the filtered Poisson process, we can simplify this sum of four diagrams for the second cumulant of the first-order self-exciting process. Examining subgraphs of each term on the right-hand side of Fig 5A reveals a connection to the linear response and mean rate of the first-order self-exciting processes. On each leg emanating from the internal branching vertex, the four terms sum to the product of two copies of the linear responses of the first-order self-exciting process N₁(t) (compare subgraphs on the top branch of the diagrams in Fig 5a with Fig 4a). Similarly, the sum of the legs coming into the branching vertex is the firing rate of N₁(t) (compare to Fig 4b). So, we will group the terms on the legs of the graph into contributions to the linear response and the mean (Fig 5b middle).

When we add the diagrams of up to order n together, we can separately re-sum each of these expansions because of the distributivity of the expectation. So, we can replace the entire series to all orders in ϵ with simpler diagrams using the full representations for the linear response and expected rate (Fig 5b). This can be proved inductively, or by rigorously deriving the Feynman rules from the cumulant generating functional (Methods: Path integral representation). This yields the following result for the second cumulant, which corresponds to the final graph at the bottom far right of Fig 5b: (30) This is the complete analytic result for the second cumulant of the self-exciting process for fluctuations around the mean field solution [19]. It can be represented by the single term on the right-hand side of Eq 30 and the corresponding single diagram (Fig 5b, right). Compare this with the filtered Poisson process, which has a diagram of the same topology but with different constituent factors (Fig 3C, middle row). The Feynman diagrams capture the form of the re-summed perturbative expansions for the cumulants, while the definitions of the vertices and edges capture the model-specific rate, , and propagator, Δ(t, t′).

One might think that the higher cumulants are generated as simply by replacing each leg of the filtered inhomogeneous Poisson process with the correct propagator, along with the rate . This would mean that the general cumulant term would be given by: (31) This is incorrect, as many important terms arising from the self-interaction would be lost.

The reason this naive generalization fails is that it neglects the higher-order responses to perturbations in the event rate. For example, the second cumulant responds to perturbations in the rate; this quadratic response impacts the third cumulant. We can see this in the third cumulant of the first-order self-exciting process: (32) The first term is the third cumulant of the inhomogeneous Poisson process. The second and third are generalizations of the terms found in the second cumulant (we have used (t ↔ t′ ↔ t′′) to denote “all other permutations of t, t′, t′′”). These terms are part of the naive expression in Eq 31. The last term is the novel one that arises due to the “quadratic response.” It appears when we compute (33) We have to take into account that the process is correlated with the rate of the process (since one is a linear function of the other!). This produces a “cascade” effect that results in the quadratic response. For the first-order process, only one step in the cascade is allowed. By introducing branching internal vertices similar to those in Fig 5, we can express these somewhat unwieldy terms with diagrams. These are shown in Fig 6. The cascade of one source spike producing three spikes in the first-order process is represented by the second diagram of Fig 6a and the cascade of one source spike producing two spikes, one of which then produces another two spikes in the first-order process, is represented by the last diagram of Fig 6a.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Diagrams corresponding to third order cumulants.

A) Diagrams corresponding to the third cumulant of the first-order self-exciting process. B) Diagrams corresponding to the third cumulant of the self-exciting process, after resumming the perturbative expansion. Nodes and edges correspond to the same terms as in Fig 5.

https://doi.org/10.1371/journal.pcbi.1005583.g006

As before, continuing to higher orders in the recursively self-exciting process would add diagrams with additional filter edges along the legs of the graphs in Fig 6a, corresponding to additional steps in the potential cascades of induced spikes. For example, the third cumulant of the second-order process, , would add diagrams with two filter edges to those of Fig 6a, with those two filter edges appearing either sequentially on the same leg of the graph or distributed among the legs of the graph. We can then use the same ideas that allowed us to re-sum the graphs representing the second cumulant. As before, we identify the expansions of the mean-field rate, , and the linear response, Δ, along individual legs of the graph and use the multilinearity of cumulants to resum those expansions to give the diagrams at the bottom of Fig 6.

Considering the re-summed graphs, we have the following result for the third cumulant: (34)

The types of diagram developed for up to the third cumulant encompass all the features that occur in the diagrammatic computations of cumulants of linearly self-exciting processes. The general rules for diagrammatically computing cumulants of this process are given in Fig 7. They are derived in general in Methods: Path integral representation. The graphs generated with this algorithm correspond to the re-summed diagrams we computed above.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Feynman rules for the self exciting process.

These rules provide an algorithm for computing the expansion of the cumulants around the mean field solution . The dots between the legs of the first two vertices indicate that there are such vertices with any number of outgoing legs greater than or equal to two.

https://doi.org/10.1371/journal.pcbi.1005583.g007

For the nth cumulant, , begin with n white external vertices labelled t_i for each i. Construct all fully connected, directed graphs with the vertex and edge elements shown in Fig 7. For each such fully connected directed graph constructed with the component vertices and edges, the associated mathematical term is constructed by taking the product of each associated factor, then integrating over the time points of internal vertices. The nth cumulant is the sum of these terms. This produces cumulants of up to third order, as recently shown by Jovanović, Hertz & Rotter [38], as well as cumulants of any order. As we show next, this procedure can also be generalized to calculate cumulants in the presence of a nonlinearity, including both thresholds enforcing positive activity (as commonly disregarded in studies of the Hawkes process) and any nonlinear input-rate transfer function.

Nonlinearities impose bidirectional coupling between different orders of activity: Nonlinearly self-exciting process

Now we include a nonlinearity in the firing rate, so that the process produces events dN/dt with a rate given by (35) We begin by considering the mean-field solution , which, if it exists, is self-consistently given by . Thus, as always the mean-field solution is given by neglecting second and higher-order cumulants of the spiking process. Next, we consider the propagator, which as above is the linear response of the rate around the mean field, given by expanding Eq 35 around the mean-field solution and examining the gain with respect to a perturbation of the rate. This propagator obeys: (36) where ϕ⁽¹⁾ is the first derivative of ϕ with respect to the input, evaluated at . We will first develop a recursive formulation of the mean-field rate and propagator, which will be required for calculating cumulants of the full process. For an arbitrary nonlinearity ϕ, we would begin by Taylor expanding around 0. For simplicity, we here consider a quadratic ϕ so that: (37) We now develop the point process dN/dt recursively at each order of the nonlinearity: (38) where M_m,n is an inhomogeneous Poisson process with rate and P_m,n is an inhomogeneous Poisson process with rate and N_0,0(t) is an inhomogeneous Poisson process with rate λ(t). To generate a set of events in N at order m in the linear term of ϕ and order n in the quadratic term of ϕ, we take events at each previous order, (m − 1, n) and (m, n − 1) and use those to develop M_m−1,n(t) and P_m,n−1(t). These, together with N_0,0(t), give . In contrast to the linear self-exciting process, the quadratic process here is recursively defined on a lattice.

Similar to the case of the linearly self-exciting process, we can use this recursive definition to develop an expansion for the mean-field rate and propagator in powers of ϵ₁ and ϵ₂. When we calculate higher-order cumulants, we will identify the expansions of the mean-field firing rate and propagator which will allow us to use them to simplify the resulting diagrams. The mean-field rate to finite order in m, n is once again given by neglecting second and higher-order cumulants of N_n,m, which allows us to take an expectation inside the quadratic term of Eq 38. Taking the expectation of both sides of this equation in the mean field approach then yields: (39) For example, (40) where (41) (42)

Similarly, the propagator (for the dynamics of the recursive process, linearized around zero) is, to finite order in m, n: (43)

To zeroth order in ϵ₂, this yields an expansion of the mean-field rate , which takes the same form as the expansion of the rate of the linearly self-exciting process (Eq 13) and admits the same graphical representation (Fig 4b). Similarly, a perturbative expansion of the linear response about the mean-field rate to zeroth order in ϵ₂ takes the same form as for the linearly self-exciting process (Eq 28) and admits the same graphical representation (Fig 4c).

To account for the nonlinear terms arising at first order and greater in ϵ₂, we will need to add another type of internal vertex in diagrammatic descriptions of the cumulants. These vertices, carrying factors of ϵ₂, will have two incoming edges and any number of outgoing edges. Each incoming edge carries the operator g* and the number of incoming edges corresponds to the order in the Taylor expansion of the nonlinearity. (It also corresponds to the order of cumulant influencing that vertex’s activity. The number of outgoing edges corresponds to the order of cumulant being influenced, locally in that subgraph.) The factor of that appears in other vertices is modified to be consistent with the mean firing rate under the quadratic nonlinearity, and will thus obey Eq 35 above.

The mean-field rate and propagator, to first order and greater in ϵ₂, can be represented diagrammatically using the new vertex (e.g., Fig 8a and 8b). Notice that these directed graphs are treelike, but with their leaves in the past. Repeating these calculations to the next order in ϵ₂ can be accomplished by taking the basic structure of Fig 8 and, along each leg entering the new vertex for ϵ₂, inserting the previous-order graphs (Figs 4a and 8a). Including higher-order terms in ϵ₁ would require inserting those graphs along the ϵ-carrying vertices of Fig 4a.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Expansion of the mean firing rate and propagator for the nonlinearly self-exciting process.

A) One of the first nonlinear terms of the expansion of the mean-field firing rate, to first order in the quadratic term of the nonlinearity and zeroth order in the linear term. The two diagrams shown correspond to the two terms in Eq 42. B) First nonlinear terms of the expansion of the propagator around the mean-field firing rate.

https://doi.org/10.1371/journal.pcbi.1005583.g008

In addition to expanding the mean field rate and propagator, we can use Eq 38 to calculate cumulants of to finite order in ϵ₁ and ϵ₂. The first nonlinear correction to the firing rate appears at first order in ϵ₂: (44) which can be represented diagrammatically using the new vertex (Fig 9a). Notice that in contrast to the corresponding graph for the mean field expansion (Fig 8a), this diagram has a “loop” (a cycle were it an undirected graph). This reflects the dependence of the rate on the second cumulant of the baseline process N_0,0. This dependence of the firing rate on higher-order spiking cumulants is a fundamental feature of nonlinearities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Corrections to the mean-field firing rate and propagator of the nonlinearly self-exciting process, to quadratic order in the nonlinearity ϕ.

A) The correction to mean field theory for the firing rate to first order in the quadratic term of ϕ. B) The full one-loop correction to the mean-field rate. C) The full one-loop correction to the propagator.

https://doi.org/10.1371/journal.pcbi.1005583.g009

Proceeding beyond the first order in both ϵ₁ and ϵ₂, we see that the expansion of each term of the nonlinearity depends on the other: (45) so that at each order in ϵ₂ we must insert the solution at the previous order in ϵ₂ and the same order in ϵ₁ (and vice versa). This recursive mixing of expansions between the linear and nonlinear terms of ϕ seems intractable. However, this joint expansion can be re-summed; for a more formal derivation, see Methods: Path integral representation. For the linear model, this re-summing left us with simple expressions for the firing rate and propagator. For the nonlinear model, re-summing leaves us with a new expansion organized by the number of loops appearing the Feynman diagrams, so it is called a loop expansion. The Feynman rules for the re-summed diagrams of the nonlinearly self-exciting process are given in Fig 10. For a quadratic nonlinearity, these rules allow us to write the one-loop contribution to the firing rate (Fig 9b): (46) where ϕ⁽²⁾ is evaluated at . These Feynman rules also allow us write down the one-loop correction to the propagator (Fig 9c): (47)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Feynman rules for the nonlinearly self exciting process.

These rules provide an algorithm for computing the expansion of the cumulants around the mean field solution . The dots between the outgoing legs of the first vertex indicate that there any number of outgoing legs greater than or equal to two. The number b of incoming edges of the second vertex correspond to its factor containing the bth derivative of ϕ, evaluated at the mean field input. The a dots between the outgoing edges of the second vertex indicate that it can have any number of outgoing edges such that a + b ≥ 3.

https://doi.org/10.1371/journal.pcbi.1005583.g010

The appearance of a loop in the Feynman diagram for the mean rate of the quadratically self-exciting process is a general feature of nonlinearities. It represents the influence of higher-order spike train cumulants on lower order ones. In order to measure that dependency, we can count the number of loops in the graphs. To do this, we add a bookkeeping variable, h. We count the number of loops by multiplying factors of h and 1/h. Each internal vertex adds a factor of 1/h and each outgoing edge a factor of h. In this way, every vertex with more than one outgoing edge will contribute a factor of h for every edge beyond the first. h thus effectively counts the order of fluctuations contributed by the vertex. For example, the mean for the linear self-exciting process has a graph with a single internal vertex and a single internal edge, so it is zeroth order in h (Fig 4b). The two-point function, however, having two edges and one internal vertex (Fig 5b), is first order in h. Similarly, the tree-level diagrams will always contribute a total of hⁿ⁻¹, where n is the order of the cumulant.

In terms of powers of h, a graph for a nth order cumulant with one loop will be equivalent to a graph for a n + 1st order cumulant with one less loop. Consider cutting one of the lines that form the loop in Fig 9b at the internal vertex and leaving it hanging. Now the graph for the one-loop correction to the mean rate appears to be a graph for a second cumulant–it has two endpoints. The power counting in terms of h, however, has not changed. The one-loop correction to the mean is of the same order in h as the tree-level second cumulant. In general, we will have that the order h^m will be given by (48) where n is the number of external vertices and l is the number of internal loops. The number of loops thus tracks the successive contributions of the higher order fluctuations. This expansion is called the “loop” expansion and is equivalent to a small-fluctuation expansion. If one can control the size of the fluctuations, one can truncate the loop expansion as an approximation for the statistics of the system. One way of doing this with networks is to insure that the interactions are so that h ∝ 1/N and the expansion becomes a system size expansion.

The Taylor expansion of an arbitrary nonlinearity ϕ could have infinitely many terms. This would lead, in the recursive formulation, to infinitely many processes {M, P, …}. Even after re-summing the recursive formulation, this would leave an infinite number of diagrams corresponding to any given cumulant. There are two ways to approach this. The first is to insist on a perturbative expansion in the nonlinear terms, e.g., only consider terms up to a fixed order in the Taylor expansion of the nonlinearity ϕ.

The second approach to controlling the order of the loop expansion is to consider a regime in which mean field theory is stable as this will also control the fluctuations, limiting the magnitude of the loop contributions [39]. We expect that the bookkeeping variable h could be related to the largest eigenvalue of the tree-level propagator, λ₁, so that m-loop corrections would scale as . The expansion then breaks down in the regime of a bifurcation or “critical point.” In this case, the spectrum of the linear response diverges, causing all loop diagrams to similarly diverge. This is a fluctuation-dominated regime in which mean field theory, along with the fluctuation expansion around it, fails. In that case, renormalization arguments can allow discussion of the scaling behavior of correlations [40].

Interaction between single-neuron nonlinearities and network structure

No new concepts are required in moving from a nonlinear self-exciting process to a network of interacting units. Each external and internal vertex must now be associated with a unique neuron index i and the integrations over time for the internal vertices must now be accompanied by summations over the indices of the internal vertices. In addition, the filter g(τ) must be expanded to include coupling across units. In general, this is given by g_ij(τ) for the coupling from neuron j to neuron i. We will consider the general model of a network of units that generate conditionally Poisson-distributed events, given an input variable. The conditional rate for unit i is given by (49) Similarly, the propagator now obeys (50) These dynamics are qualitatively the same as those of the nonlinearly self-exciting process (Eq 35 but replace the neuron’s own rate with the sum over its presynaptic inputs). Introducing these sums over neuron indices yields the complete set of rules for generating Feynman diagrams for the cumulants of this model (Fig 11).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Feynman rules for networks of stochastically spiking neurons with nonlinear input-rate transfer ϕ.

These rules provide an algorithm for computing the expansion of the cumulants around the mean field solution . The dots between the outgoing legs of the first vertex indicate that there are any number of outgoing legs greater than or equal to two. The number b of incoming edges of the second vertex correspond to its factor containing the bth derivative of ϕ, evaluated at the mean field input. The a dots between the outgoing edges of the second vertex indicate that it can have any number of outgoing edges such that a + b ≥ 3.

https://doi.org/10.1371/journal.pcbi.1005583.g011

Mid-course summary: Diagrammatic expansion reveals interplay of network structure and neural nonlinearities in driving neural activity.

To summarize, we now have in hand a set of tools—a fluctuation expansion—to compute spike train cumulants of arbitrary order in networks of linear-nonlinear-Poisson neurons. This expansion provides a systematic way to account for the synergistic dependence of lower-order activity on higher-order activity through the spiking nonlinearity and naturally incorporates the full microstructure of the neuronal network. The order of the nonlinearity (for non-polynomials, the order of its Taylor expansion) determines the single-neuron transfer gain (Fig 12 left). It also determines how activity propagates through the network (Fig 12).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Fluctuation expansion links single-neuron nonlinearities and network structure to determine network activity.

The linear response for linear neurons depends on the network structure both explicitly and implicitly, through the mean-field rates. The first nonlinear correction brings in additional explicit and implicit dependencies on the connectivity.

https://doi.org/10.1371/journal.pcbi.1005583.g012

We next provide an example of using this fluctuation expansion to compute a cumulant of spiking activity: the first nonlinear correction to the second cumulant. We will compute these using the Feynman rules (Fig 11) to construct the corresponding diagrams, from which we will write the corresponding equation.

We begin by placing external vertices corresponding to the measurement times, each with one propagator edge coming into it. For the two-point cumulant there are two external vertices (Fig 13a). The propagators coming into those external vertices can, according to the Feynman rules, arise from either a source vertex or from an internal vertex with incoming filter edges (Fig 11). If both propagators arise from a source vertex, we arrive at the tree-level diagram of Fig 5b, which provides the linear prediction for the two-point cumulant. To obtain the first nonlinear correction, we will begin by adding an internal vertex. There are two ways we can do this: with one internal vertex providing the propagators for both external vertices or with the internal vertex providing the propagator of just one external vertex (Fig 13b, top and bottom respectively).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Construction of Feynman diagrams for the first nonlinear correction to the two-point cumulant (graphs containing one loop).

In each panel, we add a new layer of vertices to the diagrams, until we arrive at a source vertex. When there are multiple potential ways to add vertices, we add diagrams to account for each of those constructions. A) External vertices corresponding to the two measurement times, with incoming propagator (Δ) edges. B) Diagrams with one internal vertex added. t₁ ↔ t₂ corresponds to switching the two external vertices in the bottom diagram; the top diagram is symmetric with respect to that switch. C) Diagrams with two layers of vertices. The top diagram finishes that of B, top. The second two arise from the second diagram of B, and each also have copies with t₁ ↔ t₂. D) Last diagrams containing one loop. The final diagrams corresponding to the one-loop correction to the second cumulant are the top two of C) and that of D).

https://doi.org/10.1371/journal.pcbi.1005583.g013

We next proceed to add another layer of vertices. The internal vertices added in each diagram of Fig 13b have incoming propagator edges. Those edges could emanate from other internal vertices or from source verticies. We will start by finishing the diagrams where they emanate from source vertices; placing these yields the top two, final diagrams of Fig 13c. We then continue constructing the diagram with an internal vertex providing the two propagators (Fig 13c, bottom). Note that if we added an internal vertex to the propagator hitting the t₂ external vertex, it would require at least two incoming filter edges (in order to obey a + b ≥ 3 per the rules of Fig 11), which would give rise to a second loop in the graph.

This last diagram has a hanging filter edge, which must arise from an internal vertex with one incoming propagator edge. We finish the diagram with that internal vertex and the source vertex providing its propagator (Fig 13d). We could not add additional internal vertices along that path, since they would either violate n + m ≥ 3 or give rise to more than one loop in the diagram (and thus belong at a higher order in the loop expansion).

Following the rules of Fig 11 while restricting ourselves to graphs with a certain number of loops (here one) thus allowed us to construct the diagrams corresponding to the first nonlinear correction to the two-point cumulant. We next write the equation corresponding to these diagrams. For each of the complete diagrams, we begin at the external vertices and, proceeding from left to right in the graph, multiply the factors corresponding to each edge and vertex together. We finish by summing over indices for all internal vertices and integrating over all internal time points. The contributions from each diagram are then added together. This yields: (51)

Demonstrating the interplay between single-neuron transfer, connectivity structure and network dynamics

Our methods predict how spike time statistics of all orders emerge from the interplay of single-neuron input-output properties and the structure of network connectivity. Here we demonstrate how these methods can be used to predict key phenomena in recurrent spiking networks: the fluctuations and stability of population activity. First, we isolate the contributions of nonlinearities in single-neuron dynamics to network activity and coding as a whole. We do so by computing “one-loop” correction terms; these correspond to the first structures in our diagrammatic expansion that arise from nonlinear neural transfer. The one-loop corrections provide for the dependence of nth order spiking cumulants on n + 1st order cumulants. Predictions that would be made by linearizing neural dynamics, as in classic approaches for predicting pairwise correlations [19, 22, 41] and recent ones for higher-order correlations [38], are described as “tree-level.” We show how these one-loop corrections, which give new, explicit links between network structure and dynamics (Fig 12), predict spiking statistics and stability in recurrent networks.

1. Recurrent spike-train correlations drive firing statistics and stability in nonlinear networks

In our analysis of the impact of nonlinear neural transfer on network dynamics, a principal finding was that spike correlations could affect firing rates, as described by the one-loop correction to the mean-field firing rates. In this section we illustrate the importance of this effect in a class of networks under intensive study in neuroscience: randomly connected networks of excitatory and inhibitory cells. We began with a network for which we expect classical theoretical tools to work well, taking the neurons to have threshold-linear transfer functions ϕ(x) = α⌊x⌋. Here, as long as the neurons do not receive input fluctuations that push their rates below this threshold, the tree-level’ theory that takes transfer to be entirely linear should work well. We then move on to consider nonlinear effects.

As in our original motivational example, we took network connectivity to be totally random (Erdös-Rényi), with p_EE = 0.2 and p_EI = p_IE = p_II = 0.5. The magnitude of all connections of a given type (E − E, etc.) was taken to be the same and the time course of all network interactions was governed by the same filter (with τ = 10 ms), so that g_ij(t) = W_ij g(t). (The matrix W contains synaptic weights.) The net inhibitory input weight on to a neuron was, on average, twice that of the net excitatory input weight.

We examined the spiking dynamics as the strength of synaptic weights proportionally increased (Fig 14a) and studied network activity with using both theory and direct simulation. Due to the high relative strength of inhibitory synapses in the network, firing rates decreased with synaptic weight (Fig 14d). The magnitude of spike train covariances (reflected by the integrated autocovariance of the summed excitatory population spike train) increased (Fig 14e). These changes were also visible in raster plots of the network’s activity (Fig 14b and 14c).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14. Dynamics approaching the firing-rate instability in threshold-linear networks.

A) Threshold-linear input-rate transfer function. B,C) Raster plots of 1 second realizations of activity for weak and strong synaptic weights. Neurons 0–199 are excitatory and 200–240 are inhibitory. B) (W_EE, W_EI, W_IE, W_II) = (.025, −.1, .01, −.1) mV. C) (W_EE, W_EI, W_IE, W_II) = (.2, −.8, .08, −0.8) mV. D-F) Average firing rate of the excitatory neurons (D), integral of the auto-covariance function of the summed population spike train (E), and spectral radius of the stability matrix of mean field theory (F) vs. excitatory-excitatory synaptic weight. While excitatory-excitatory weight is plotted on the horizontal axis, all other synaptic weights increase proportionally with it. Black lines: tree-level theory: mean-field firing rates and covariance computed by linearizing dynamics around it, for each value of synaptic weights. Dots: simulation.

https://doi.org/10.1371/journal.pcbi.1005583.g014

At a critical value of the synaptic weights, the mean field theory for the firing rates loses stability (Fig 14f). The location of this critical point is predicted by the linear stability of the dynamics around the mean-field rate; the spectrum of the propagator Δ(ω) = (I − ϕ⁽¹⁾ g(ω))⁻¹ diverges when the spectral radius of ϕ⁽¹⁾ g is ≥1. (This is also the point where the spectral radius of the inverse propagator crosses zero.) Until that critical point, however, the tree-level predictions for both firing rates and spike train covariances (i.e., mean field theory and linear response theory) provided accurate predictions (Fig 14d and 14e).

We next give a simple example of how nonlinearities in neural transfer cause this standard tree-level theory (mean-field formulas for rates and linear response theory for covariances) to fail–and how tractable one-loop corrections from our theory give a much-improved description of network dynamics. We take the same network as above, but replace neurons’ threshold-linear transfer functions with a rectified power law ϕ(x) = α⌊x⌋^p (Fig 15a). This has been suggested as a good description of neural transfer near threshold [42–46]. For simplicity, we take the power law to be quadratic (p = 2). As we increased synaptic weights, the tree-level theory qualitatively failed to predict the magnitude of spike train covariances and firing rates (Fig 15d and 15e black curve vs. dots). This occurred well before the mean-field firing rates lost stability (Fig 15f, black).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 15. Correlation-driven instability in nonlinear networks.

A) Threshold-quadratic input-rate transfer function. B,C) Raster plots of 6 second realizations of activity for weak and strong synaptic weights. Neurons 0–199 are excitatory and 200–240 are inhibitory. B) (W_EE, W_EI, W_IE, W_II) = (.025, −.1, .01, −.1) mV. C) (W_EE, W_EI, W_IE, W_II) = (1.5, −6, .6, −6) mV. D-F) Average firing rate of the excitatory neurons (D), integral of the auto-covariance function of the summed population spike train (E), and spectral radius of the stability matrix of mean field theory (F) vs. excitatory-excitatory synaptic weight. While excitatory-excitatory weight is plotted on the horizontal axis, all other synaptic weights increase proportionally with it. Black line: tree-level theory. Red line: one-loop correction accounting for impact of the next order (pairwise correlations’ influence on mean and triplet correlations’ influence on pairwise). Dots: simulation. All dots after the one-loop spectral radius crosses 1 represent results averaged over the time period before the activity diverges.

https://doi.org/10.1371/journal.pcbi.1005583.g015

Higher-order terms of the loop expansion described above (Nonlinearities impose bidirectional coupling between different orders of activity) provide corrections to mean field theory for both firing rates and spike train correlations. These corrections represent coupling of higher-order spike train cumulants to lower order cumulants. In the presence of an input-rate nonlinearity, for example, synchronous (correlated) presynaptic spike trains will more effectively drive postsynaptic activity [7, 47]. This effect is described by the one-loop correction to the firing rates (Fig 9).

The one-loop correction for the mean field rate of neuron i in a network is given by the same diagram as the one-loop correction for the nonlinearly self-exciting process, Fig 9, but interpreted using the network Feynman rules (Fig 11). This yields: (52) where t₀ is the initial time. This correction was, on average, positive (Fig 15d for excitatory neurons; also true for inhibitory neurons). Similarly to firing rates, the loop expansion provides corrections to higher-order spike train cumulants. The one-loop correction to the spike train covariances (Eq 51, derived in Fig 13) accounts for the impact of triplet correlations (third joint spike train cumulants) on pairwise correlations and provided an improved prediction of the variance of the population spike train as synaptic weights increased (Fig 15e).

Since the one-loop correction to the firing rates could be large, we also asked whether it could impact the stability of the firing rates–that is, whether pairwise correlations could, through their influence on firing rates through the nonlinear transfer function, induce an instability. This is a question of when the eigenvalues of the propagator diverge—or equivalently, when the eigenvalues of the inverse propagator cross zero. The one-loop correction to the inverse propagator is given by the “proper vertex” obtained by amputating the outside propagator edges of the one-loop correction to the propagator; or equivalently, calculated from the Legendre transform of the cumulant-generating functional [39]. We can heuristically derive the one-loop stability correction as follows.

The full propagator, Δ, obeys the expansion (53) where is the tree-level propagator, Δ₁ is the one-loop correction, two-loop corrections are collected in Δ₂, and so on (Fig 16a). Notice from the diagram (Fig 16a, second diagram on the right-hand side) that the one-loop correction to the propagator begins and ends with the tree-level propagator, . We will label the bubble in the middle of the diagram Γ₁. The first two-loop correction is a chain of loops (Fig 16a), and so can also be factored as . We can represent this factorization diagrammatically by pulling out the tree-level propagator and the loop Γ₁ (Fig 16b). Just as at two-loop order we were able to factor out a factor and obtain the expansion of the propagator to one loop, continuing to higher-orders in the loop expansion of the full propagator would all the rest of the full propagator with factors of in front. The remaining terms would have factors starting with the two-loop correction, and so forth.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 16. Calculation of the one-loop stability correction.

A) Loop expansion of the full propagator. B) Factorization of the loop and resumming of the full propagator after that factorization.

https://doi.org/10.1371/journal.pcbi.1005583.g016

Pulling out all terms of Eq 53 that begin with and summing them allows us to write (Fig 16b): (54) We now truncate at one loop and operate on both sides with the inverse of the tree-level propagator: (55) (56) revealing that −Γ₁ is the one-loop correction to the inverse propagator. From the Feynman rules (Fig 11), that factor is: (57) where denotes the second derivative of the transfer function of neuron j, evaluated at its mean-field input (and similar for ). The eigenvalues of this provide a correction to the stability analysis based on the tree-level propagator. This predicts that the firing rates should lose stability significantly before the bifurcation of the mean field theory (Fig 15f, red vs. black). Indeed, we saw in extended simulations that the spiking network could exhibit divergent activity even with synaptic weights that mean field theory predicted should be stable (Fig 1c). In summary, mean field theory can mis-predict the bifurcation of the rate of spiking models since it fails to capture the impact of correlations on firing rates through nonlinear transfer functions.

2. Impact of connectivity structure on correlation-driven instabilities in nonlinear networks

Recent work has shown that cortical networks are more structured than simple Erdős-Rényi networks (e.g., [48–53]). One feature of cortical networks is a broad spread of neurons’ in- and out-degree distributions (i.e., the distributions of the number of synaptic inputs that each neuron receives or sends); another is broadly spread synaptic weights. These network properties, in turn, can have a strong impact on population activity [54–58]. Here, we illustrate the link between network structure and activity in the presence of nonlinear neural transfer. To generate structured networks, we began with the type of excitatory-inhibitory networks discussed in the previous section, but took the excitatory-excitatory coupling to have both heavy-tailed degree and weight distributions. Specifically, we took it to have truncated, correlated power law in- and out-degree distributions (Methods: non-Erdős-Rényi network model). We then took to the synaptic weights to be log-normally distributed [48, 59]. For simplicity, we took the location and scale parameters of the weight distribution to be the same.

We then examined the network dynamics as the location and scale of the excitatory-excitatory synaptic weights increased. For each mean weight, we sampled the excitatory-excitatory weights from a lognormal distribution with that mean and variance. The excitatory-inhibitory, inhibitory-excitatory and inhibitory-inhibitory weights remained delta-distributed. Each such network specified a weight matrix W, which allowed the methods described previously for computing tree-level and one-loop rates, covariances and stability to be straightforwardly applied. For strong and broadly distributed synaptic weights (Fig 17b), the network exhibited a similar correlation-induced instability as observed in the Erdős-Rényi network (Fig 17c) even though mean field theory predicted that the firing rates should be stable (Fig 17f, black vs. red curves). As synaptic weights increased from zero, the mean field theory for firing rates provided a misprediction (Fig 17d, black line vs. dots) and the linear response prediction for the variance of the population spike train also broke down (Fig 17e, black line vs. dots). The one-loop corrections, accounting for the impact of pairwise correlations on mean rates and of triplet correlations on pairwise correlations, yielded improved predictions (Fig 17d and 17e red lines) and a much more accurate prediction for when firing rates would lose stability (Fig 17c and 17f). These effects were similar to those seen in Erdős-Rényi networks (Fig 15), but the transition of the firing rates occurred sooner, both for the mean field (because of the effect of the weight and degree distributions on the eigenvalues of the weight matrix) and one-loop theories (because of the impact of the correlations on the firing rates).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 17. Correlation-driven instability in a non-Erdős-Rényi network with broadly distributed excitatory-excitatory weights.

A) Threshold-quadratic input-rate transfer function. B) Histogram of excitatory-excitatory synaptic weights with location parameter of 1.42 (mean of.29 mV), corresponding to the simulation in panel C. C) Raster plots of 6 second realizations of activity. Neurons 0–199 are excitatory and 200–240 are inhibitory. (W_EE, W_EI, W_IE, W_II) = (1.125, −4.5, .45, −4.5) mV. D-F) Average firing rate of the excitatory neurons (D), integral of the auto-covariance function of the summed population spike train (E), and spectral radius of the stability matrix of mean field theory (F) vs. excitatory-excitatory synaptic weight. While the mean excitatory-excitatory weight is plotted on the horizontal axis, all other synaptic weights increase proportionally with it. Black line: tree-level theory. Red line: one-loop correction. Dots: simulation. If a simulation exhibits divergent activity, the spike train statistics are averaged over the transient time before that divergence for visualization.

https://doi.org/10.1371/journal.pcbi.1005583.g017

3. Exponential single-neuron nonlinearities

In the previous section, we investigated how a non-Erdős-Rènyi network structure could amplify the one-loop corrections by increasing spike train correlations. We now examine a different single-neuron nonlinearity: ϕ(x) = αe^x, which is the canonical link function commonly used to fit GLM point process models to spiking data [26].

As before, we take the mean synaptic weight onto each neuron in the network to be 0. First, we take excitatory and inhibitory neurons to have the same baseline drive, λ_E = λ_I = −1.5. As we scale synaptic weights, we see that the one-loop correction is small compared to the tree-level theory for the firing rates, population variances and stability analysis (Fig 18a–18c, red vs. black lines). It nevertheless provides an improved correction for the variance of the excitatory population spike train (Fig 18b, between 1.5 and 2 mV synaptic weights). The bifurcation of the one-loop theory is close to the bifurcation of the mean field theory, and before that point the mean field theory and one-loop corrections both lose accuracy (Fig 18a and 18b). This makes sense: when the mean field theory fails, the only reason that the one-loop correction to the rates would be accurate is if all third- and higher-order spike train cumulants were small. Those higher-order correlations are not small near the instability.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 18. Stability of a network with exponential transfer functions.

A) Mean firing rate of the excitatory neurons. B) Integral of the auto-covariance function of the summed population spike train. C) Spectral radius of the stability matrix of mean field theory, all (A-C) vs. excitatory-excitatory synaptic weight. While the mean excitatory-excitatory weight is plotted on the horizontal axis, all other synaptic weights increase proportionally with it. Black line: tree-level theory. Red line: one-loop correction. Dots: simulation. If a simulation exhibits divergent activity, the spike train statistics are averaged over the transient time before that divergence for visualization.

https://doi.org/10.1371/journal.pcbi.1005583.g018

Next, we broke the symmetry between excitatory and inhibitory neurons in the network by giving inhibitory neurons a lower baseline drive (λ_I = −2.). This shifted the bifurcation of the mean field theory and the one-loop correction to much higher synaptic weights (Fig 19c). For intermediate synaptic weights, we saw that the one-loop correction provided a better match to simulations than the tree-level theory (Fig 19a and 19b, between 1 and 1.5 mV synaptic weights). For stronger synapses, however, the simulations diverged strongly from the tree-level and one-loop predictions (Fig 19a and 19b, around 1.5 mV synaptic weights). In principle, we could continue to calculate additional loop corrections in an attempt to control this phenomenon. The exponential has arbitrary-order derivatives, however, preventing a perturbative expansion of the nonlinear terms—suggesting a renormalization approach [39], which is beyond the scope of this article. In sum, with an exponential transfer function, we saw that for intermediate synaptic weights, the one-loop correction improved on the tree-level theory. For strong enough synaptic weights, however, both failed to predict the simulations. How soon before the mean-field bifurcation this failure occurred depended on the specific model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 19. Failure of one-loop corrections with exponential transfer functions.

A) Mean firing rate of the excitatory neurons. B) Integral of the auto-covariance function of the summed population spike train. C) Spectral radius of the stability matrix of mean field theory, all (A-C) vs. excitatory-excitatory synaptic weight. While the mean excitatory-excitatory weight is plotted on the horizontal axis, all other synaptic weights increase proportionally with it. Black line: tree-level theory. Red line: one-loop correction. Dots: simulation. If a simulation exhibits divergent activity, the spike train statistics are averaged over the transient time before that divergence for visualization.

https://doi.org/10.1371/journal.pcbi.1005583.g019

Relationship to other theoretical methods

Dynamical mean field theory is a classic tool for analyzing the rate models with Gaussian-distributed synaptic weights, which reveals a transition to chaotic rate fluctuations with strong connectivity [60]. This proceeds, briefly, by taking the limit of large networks and replacing interactions through the quenched heterogeneity of the synaptic weights by an effective Gaussian process mimicking their statistics. Recent extensions have incorporated a number of simple biological constraints, including non-zero firing rates [61–63] and certain forms of cell type-specific connectivity [61, 64–66]. In this framework, spiking is usually only described in the limit of slow synapses as additive noise in the rates, which can shift the transition to chaotic rate fluctuations to higher coupling strengths and smooth the dynamics near the transition [61, 67].

An alternative approach is to start from the bottom up: to posit an inherently stochastic dynamics of single neurons and specify a finite-size network model, and from these derive a set of equations for statistics of the activity [68, 69]. This approach provides a rigorous derivation of a finite-size rate model as the mean field theory of the underlying stochastic activity, as well as the opportunity to calculate higher-order activity statistics for the activity of a particular network [40, 70–72]. This is the approach taken here with the popular and biologically motivated class of linear-nonlinear-Poisson models. A similar approach underlies linear response theory for computing spike train covariances [73], which corresponds to the tree level of the loop expansion presented here. For integrate-and-fire neuron models receiving Poisson inputs, Fokker-Planck theory for the membrane potential distributions can be used to calculate the linear response function of an isolated neuron [74], which together with the synaptic filter and weight matrix determines the propagator Δ.

Convergence and truncation of the loop expansion

The loop expansion is organized by the dependence of lower-order activity cumulants to higher-order ones. The first-order (tree level) description of the nth activity cumulant does not depend on higher-order cumulants. One-loop corrections correspond to dependence of the order n cumulants on the tree-level n + 1-order cumulants, two-loop corrections correspond to dependence of the order n cumulant on the tree-level n + 1 and n + 2 order cumulants and so on. This coupling arises from the nonlinearity of the single-neuron transfer function ϕ (Results: Nonlinearities impose bidirectional coupling between different orders of activity: nonlinearly self-exciting process; Methods: Path integral representation). When the transfer function is linear at the mean-field rates, the tree-level theory provides an accurate description of activity (Fig 14). This corresponds to the second- and higher-order derivatives of the transfer function ϕ with respect to the total input, evaluated at the mean-field rates, being zero. When ϕ has non-zero second- or higher derivatives at the mean-field rates, higher orders of the loop expansion can be important. The magnitude of n-loop corrections depends on two things: the magnitude of the n + 1-order tree-level activity cumulants and the magnitude of the n + 1st derivative of ϕ at the mean-field rates (i.e., the strength of the coupling to that cumulant).

Recent work has shown that the the magnitude of spike train correlations depend on the motif structure of the network (Fig 17; [41, 55, 56, 75]), as well as on the correlation structure of the inputs it receives. For a particular model network of interest, the accuracy of a truncated loop expansion can be evaluated empirically by comparing against simulations. Computing higher-order loop corrections can quickly become unwieldy. The ratio of the n-loop correction to a particular spike train cumulant to the n + 1 loop correction can provide an estimate (but not a guarantee) of whether truncating at n loops is sufficient. If the loop contributions scale with a small parameter such as the inverse system size, this provides a natural justification for truncation. This can occur if the tree-level cumulants scale with 1/N, which happens if synaptic weights scale as 1/N [23, 71], or (in the thermodynamic limit) if connections are strong but sparse [76], or if strong, dense, balanced excitation and inhibition actively cancel correlations [77]. Finally, if the system lies close to a bifurcation of the mean field theory so that the eigenvalues of the propagator diverge, then the mean field theory and this expansion around it can also fail. In that case, renormalization arguments can allow the discussion of the scaling behavior of correlations [40].

Dynamics and stability in spiking networks

A new potential impact of correlations on population coding.

Many studies have examined the impact of “noise correlations” on population coding, examining the vector of neural responses. If all responses are conditionally independent given the stimulus, the distribution of responses to a particular stimulus is spherical. The discriminability of the responses to two stimuli corresponds to the area of overlap of those multi-neuron response distributions. To tree level in the loop expansion of the population responses, correlations stretch that response distribution. These correlations can either improve or lower coding performance, depending on how they relate to the stimulus-evoked responses [12, 16–18, 85]. In the presence of a nonlinear transfer function, a further potential impact of correlations is to change neurons’ mean activities (Fig 15. This corresponds to a translation of the multi-neuron response distributions (Fig 20, bottom), which could, in principle, either increase or decrease their discriminability (Fig 20, bottom).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 20. Potential impacts of correlations on coding in a presence of a nonlinearity.

A) The independent Poisson assumption for neurons gives rise to uncorrelated distributions of population activity. B) Correlations could increase or decrease the overlap of those distributions by stretching them, decreasing or increasing the decoding performance (top and bottom, respectively). C) The impact of correlations on the mean responses can shift those distributions, potentially counteracting the impact of stretching the distributions (as shown), or exaggerating it.

https://doi.org/10.1371/journal.pcbi.1005583.g020

Non-Erdős-Rényi network model

For Fig 17, we generated the excitatory-excitatory connectivity with a truncated power-law degree distribution. The marginal distributions of the number of excitatory-excitatory synaptic inputs (in-degree) or outputs (out-degree) obeyed: (58) where d is the in- or out-degree. Parameter values are contained in Table 1; C₁ and C₂ are normalization constants to make the degree distribution continuous at the cutoff L₁. The in- and out-degree distributions were then coupled by a Gaussian copula with correlation coefficient ρ to generate in- and out-degree lists. These lists generated likelihoods for each possible connection proportional to the in-degree of the postsynaptic neuron and the out-degree of the presynaptic neuron. We then sampled the excitatory-excitatory connectivity according to those likelihoods.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Model parameters (unless otherwise specified in text).

https://doi.org/10.1371/journal.pcbi.1005583.t001

Path integral representation

Here we outline the derivation of a path integral formalism for a network of processes with nonlinear input-rate transfer, following methods developed in nonequilibrium statistical mechanics [86–91]. We will begin by developing the formalism for a simple model, where a spike train is generated stochastically with intensity given by some input ν(t). We will specify the cumulant generating functional of the spike train given ν(t) below.

The general strategy is to introduce an auxiliary variable, called the “response variable,” whose dynamics will determine how information from a given configuration (e.g. the spike counts n(t)) at one time will effect future configurations of the dynamics. Introducing the response variable allows us to write the probability density functional for the process in an exponential form. The integrand of that exponential is called the “action” for the model, which can then be split into a “free” action (the piece bilinear in the configuration and response variables) and an “interacting” one (the remainder). Cumulants of the process can then be computed, in brief, by taking expectations of the configuration and response variables and the interacting action against the probability distribution given by the free action.

Let n(t) be the number of spike events recorded since some fiducial time t₀. In a time bin dt, Δn events are generated with some distribution p(Δn) and added to n(t). Let the events generated in any two time bins be conditionally independent given some inhomogeneous rate ν(t), so that p(Δn) = p(Δn|ν). So, assuming that initially n(t₀) = 0, the probability density functional of the vector of events over M time bins is: (59) where is the Laplace transform of p(Δn_i|ν_i) and is the cumulant generating functional for the auxiliary variable. In the third step we have written the distribution of p(Δn_i) as the inverse Laplace transform of the Laplace transform. The Laplace transform variable is our auxiliary response variable. In the fourth step we identified the Laplace transform of the probability density functional as the moment-generating functional, so that is the cumulant generating functional of the spike count. Note that these are complex integrals. The contour for the integration over is parallel to the imaginary axis.

Taking the continuum limit M → ∞, dt → 0 then yields the probability density functional of the spike train process : (60) where and and we suppress the conditional dependence of on ν(t). In the continuum limit the integral is a functional or path integral over realizations of . We will call the negative exponent of the integrand in Eq 60 the action: (61) We have slightly abused notation here in that a factor of 1/dt has been absorbed into . We will justify this below.

We have not yet specified the conditional distribution of the events given the input ν(t), leaving unspecified. Here, we will take the events to be conditionally Poisson [92], so that (62) (In the continuum limit, the rate ν(t) allowed us to absorb the factor of 1/dt into W. A finite size time bin would produce ν(t)dt events in bin dt.)

This representation of the probability density functional yields the joint moment-generating functional (MGF) of and : (63) and the moment-generating functional of : (64)

The above strictly applies only to the inhomogeneous Poisson process. This formalism is adapted to the self-exciting process by introducing conditional dependence of the rate ν(t) on the previous spiking history. In the discrete case, before taking the limit M → ∞, we say that the rate ν_i = ϕ[n_i−], where ϕ is some positive function and n_i− indicates all spiking activity up to but not including bin i. This requirement is equivalent to an Ito interpretation for the measure on the stochastic process . Because of this assumption, the previous derivation holds and we can write (65) where . In the continuum limit, there is an ambiguity introduced by the appearance of the time variable t in both and . This is resolved in the definition of the measure for the functional integral and affects the definition of the linear response (below). Again, this is a manifestation of the Ito assumption for our process.

The specific model used in this paper assumes a particular form for the argument of ϕ. We assume that the input is given by (66) where g(t) is a filter that defines the dynamics of the process in question and λ(t) is an inhomogeneous rate function. The result is that the action for nonlinearly self-exciting process is given by (67)

The only extension required to move from the above action to the network model is to introduce indices labelling the neurons and couplings specific for each neuron pair. Nothing of substance is altered in the above derivation and we are left with (68)