empty square/.style= /tikzfeynman/square dot, /tikz/fill=none ,

††thanks: co-senior authors††thanks: co-senior authors

# Metastability in networks of nonlinear stochastic integrate-and-fire neurons

Siddharth Paliwal  Department of Neurobiology and Behavior, Stony Brook
University, Stony Brook, NY, 11794, USA  Gabriel Koch Ocker<EMAIL_ADDRESS>Department of Mathematics and Statistics, Boston University, Boston, MA,
02215, USA  Braden A. W. Brinkman<EMAIL_ADDRESS>Department
of Neurobiology and Behavior, Stony Brook University, Stony Brook, NY, 11794,
USA

###### Abstract

Neurons in the brain continuously process the barrage of sensory inputs they
receive from the environment. A wide array of experimental work has shown that
the collective activity of neural populations encodes and processes this
constant bombardment of information. How these collective patterns of activity
depend on single neuron properties is often unclear. Single-neuron recordings
have shown that individual neural responses to inputs are nonlinear, which
prevents a straightforward extrapolation from single neuron features to
emergent collective states. In this work, we use a field theoretic formulation
of a stochastic leaky integrate-and-fire model to study the impact of
nonlinear intensity functions on macroscopic network activity. We show that
the interplay between nonlinear spike emission and membrane potential resets
can i) give rise to metastable transitions between active firing rate states,
and ii) can enhance or suppress mean firing rates and membrane potentials in
opposite directions.

††preprint: APS/123-QED

##  I  Introduction

Populations of neurons can showcase a wide array of complex collective
activity patterns that underlie sensory and information processing in the
brain. Understanding how these macroscopic patterns of activity emerge from
the properties of individual neurons has been a central question in
neuroscience. Recently, metastable activity—sharp stochastic- or input-driven
changes in the firing patterns of the network—has come under investigation as
a potential means by which populations could flexibly encode sensory
information [1, 2]. These states are frequently observed as sub-populations of
neurons switching between up (high activity) and down (low activity) states
[3, 4, 5, 6, 7], or transitions between multiple clusters of neurons [8, 9,
10]. Models of up-down transitions commonly describe the activity in these
states as bistable attractors in dynamical models of network activity [11, 12,
13, 14, 15, 16]. Similarly, models with multiple clusters can produce
multistable attractors that describe metastability between multiple states
[17, 18, 19, 20, 21].

Classically, theories for describing emergent macroscopic population activity
(e.g., the Wilson-Cowan [22, 23] or Amari-Grossberg equations [24, 25, 26])
are commonly understood as coarse grained models of the underlying populations
of neurons. These equations are derived under the implicit assumption of a
separation of timescales between the synaptic and intrinsic dynamics of
neurons [27]. Additionally, these models describe the mean population
activity, ignoring the role that stochastic fluctuations play in shaping the
emergent activity of a population of neurons. Some models introduce
stochasticity into firing rate models by introducing an external fluctuating
current input, with variance sometimes self-consistently matched to the firing
rates of the network to mimic a recurrent Poissonian spiking input. However,
this phenomenological addition of noise assumes that correlations in spiking
activity are negligible, precluding strong coordination within the network.

To properly understand the influence of recurrent spiking activity, we follow
recent modeling efforts that employ stochastic field theory representations of
the spiking network dynamics without inserting phenomenological noise to mimic
the effects of spiking variability [28, 29, 30, 31, 32]. In particular, we use
a stochastic spiking network that incorporates a hard reset of the neurons’
membrane potentials after they emit a spike [33]. This hard reset aids in
stabilizing network activity, allowing the stochastic neurons to fire
guaranteed spikes if their input is large enough, yet remain stable because
their membrane potentials are always reset after a spike.

The stability of a network and the possible patterns of collective activity it
admits is also intrinsically tied to nonlinearities in neural activity, even
at the single neuron level. For example, neurons are often characterized by an
intensity (or transfer) function that characterizes their instantaneous rate
of firing given their current membrane potential [34, 35, 28, 36, 37, 38, 39,
40]. This intensity function cannot be measured directly; instead, the firing
rate as a function of mean membrane potential has been estimated by fitting
single neuron data. The measured nonlinearities are assumed to be a noisy
version of the underlying intensity function, giving insight into the
appropriate families of intensity functions to use in theoretical and
computational models. For example, in the primary visual cortex the
nonlinearities are fit well by threshold-power law functions of the form
$\lfloor V-\theta\rfloor_{+}^{\alpha}$, where $\lfloor x\rfloor_{+}=x$ if
$x>0$ and $0$ otherwise, and $\alpha$ is an exponent [37, 38, 39, 40].
Observed exponents $\alpha$ typically fall in the range of 2 $\sim$ 5, larger
than the rectified linear units often used to model neural firing in
mathematically tractable models of stochastic spiking neurons [41, 42, 33].

The rectifying nature of the threshold power-law stands in contrast with the
exponential nonlinearity, an alternative modeling choice that has been used
extensively in fitting point process generalized linear models (GLMs) and
generalized integrate-and-fire models to spiking activity data [43, 44, 45,
46, 47, 29, 48, 49, 50, 31]. The use of exponential nonlinearities in fitting
GLMs is primarily for technical convenience, but it motivates us to also
investigate whether the collective activity generated by such networks
approximates the activity of threshold-power law networks. As shown in Fig. 1,
we find that both power law and exponential nonlinearities can fit intensity
functions well across cell types and layers of visual cortex in mouse data
from the Allen Institute for Brain Science [47], but we will show that the two
classes of nonlinearities exhibit several key differences in the collective
activity they can generate.

To show this, we use a stochastic field theory formalization [51, 52, 53, 54,
55, 56, 33] of a biologically realistic model of stochastic leaky and fire
(sLIF) neurons with nonlinear intensity functions. We study its dynamics to
compare the impact of threshold power-law and exponential intensity functions
on population activity. We first calculate the mean-field phase diagram of the
stochastic spiking network as a function of its synaptic connectivity strength
and external input. We show that superlinear threshold-power law intensity
functions, coupled with a hard reset of the membrane potential after spiking,
enable metastability between i) a quiescent and active state for subthreshold
currents and ii) two active states for super-threshold current. In the latter
case the networks can stochastically switch between these two active states.
Moreover, we show that the quiescent steady state is a result of the
rectifying nature of the intensity. Networks with exponential intensity
functions exhibit a monostable regime of continuously varying firing rate or a
metastable regime between low and high firing rates.

To capture the impact of fluctuations on population activity, we go beyond
mean-field theory to compute fluctuation corrections. The fluctuation
corrections to the mean-field theory come from the nonlinear firing intensity
and the spike reset. While the spike reset always suppresses activity, the
nonlinear firing intensity promotes activity when the nonlinearity is
superlinear and suppresses activity when it is sublinear [28]. We show that
this interplay between the nonlinear intensity and spike reset fluctuations
shapes the overall dynamics of the network.

The paper is organized as follows: in Section II.1 we introduce the model and
the corresponding mean-field equations using which we can obtain the phase
diagram for the network. In Section II.2, II.3 we obtain the mean-field phase
diagrams for the two nonlinear intensity functions and find that they can have
at most two stable steady-state solutions. Finally, we analyze the impact of
fluctuations on the mean firing rate and membrane potential of the network.
This helps us dissociate the impact of the two sources of fluctuation
corrections: the spike reset and nonlinear firing intensity.

Figure 1: a) Voltage trace showcasing multiple hard-resets after a spike is
emitted by a neuron in the model. b) Voltage trace of an example neuron from
the cell electrophysiological recording showcasing the full spike generation
and post spike hyper-polarization. The model differs from the true dynamics of
the cell since the details of the spike generation and hyper-polarization are
replaced by the hard-reset. c) Mean membrane potential vs firing rate data for
an example cell (empty red circles), with the power law (pink) and exponential
(dark blue) fits. The text shows the estimated parameters for both fits. d)
The average estimated exponent of the threshold power law fit for different
CRE lines. The inset shows the estimated gain for the fit. e) Same as d but
for the exponential fit showing the threshold and gain (inset).

##  II Results

###  II.1 Network model and mean-field equations

We model a network of $N$ stochastic leaky-integrate and fire neurons (sLIF).
The membrane dynamics are governed by the stochastic differential equation:

$\displaystyle\frac{dV_{i}(t)}{dt}$
$\displaystyle=-\frac{V_{i}(t)}{\tau}+\sum_{j=1}^{N}J_{ij}\dot{n}_{j}(t)+\mathcal{E}_{i}-\dot{n}_{i}(t^{+})V_{i}(t)$
$\displaystyle\dot{n}_{i}(t)dt$ $\displaystyle\sim{\rm
Poiss}\left[\phi(V_{i}(t))dt\right].$ (1)

where $V_{i}(t)$ is the membrane potential of neuron $i$ at time $t$, $J_{ij}$
is the connection strength from neuron $j$ to neuron $i$ (equal to $J/(pN)$
with probability $p$, the connection probability, and $0$ otherwise), and
$\mathcal{E}_{i}$ is the net external current that each neuron receives. We
non-dimensionalize the model by setting $\tau=1$. The spike trains
$\dot{n}_{i}(t)$ are conditionally Poisson, with instantaneous firing rate
$\phi(V_{i}(t))$; i.e., the instantaneous firing rate is a function of the
neuron’s membrane potential. After the emission of a spike, the last term
($-\dot{n}_{i}(t^{+})V_{i}(t)$) resets the membrane voltage to $0$; see Fig.
1a). We aim to characterize the steady-states for rectified power-law
intensity functions, $\phi(V)\propto\lfloor V-\theta\rfloor^{\alpha}_{+}$ with
exponent $\alpha>0$ and exponential intensity functions $\phi(V)\propto
e^{V-\theta}$, where $\lfloor x\rfloor_{+}=x$ if $x>0$ and $0$ otherwise.
Here, $\theta$ is the activation threshold in power-law networks and a “soft”
threshold in exponential networks. In both cases it is not the threshold at
which spikes are guaranteed to occur, but above which the probability of
spiking increases sharply.

To motivate the use of these particular nonlinear intensity functions, we fit
the average membrane potential vs firing rate to individual cells from the
openly available electrophysiology data from the Allen Institute for Brain
Science [47] (Appendix A). The dynamics of individual cells differs from the
model dynamics where the details of the spiking and hyper-polarization are
replaced by the hard-reset (Fig. 1b). The fitted exponents fall in the range
of 2 $\sim$ 4 as has been observed previously, with the mean exponent of 3.2
(Fig. 1c-e).

The stochastic dynamics in (1) gives rise to a probability distribution for
the membrane potential and firing rate for the network. The joint moment
generating functional (MGF) for membrane potentials and spike trains can be
expressed as a path integral in the Martin-Siggia-Rose-De Dominicis-Janssen
(MSRDJ) formalism, described by the action $S[\tilde{V},V,\tilde{n},\dot{n}]$:

$\displaystyle Z[\tilde{j},j,\tilde{h},h]$
$\displaystyle=\int\mathcal{D}\tilde{V}\,\mathcal{D}V\mathcal{D}\tilde{n}\,\mathcal{D}\dot{n}\,\exp\left(-S[\tilde{V},V,\tilde{n},\dot{n}]+\tilde{j}\cdot
V+j\cdot\tilde{V}+\tilde{h}\cdot\dot{n}+h\cdot\tilde{n}\right),$
$\displaystyle S\left[\tilde{V},V,\tilde{n},\dot{n}\right]$
$\displaystyle=\sum_{i=1}^{N}\left\\{\tilde{V_{i}}(t)\left[{\dot{V}_{i}(t)}+V_{i}(t)-\mathcal{E}_{i}+\dot{n}_{i}(t^{+})V_{i}(t)-\sum_{j=1}^{N}J_{ij}\dot{n}_{j}(t)\right]+\tilde{n}_{i}(t)\dot{n}_{i}(t)-\left(e^{\tilde{n}_{i}(t)}-1\right)\phi(V_{i}(t))\right\\}.$
(2)

where $\tilde{x}\cdot y=\sum_{i}\int dt~{}\tilde{x}_{i}(t)y_{i}(t)$,
$\tilde{V},\tilde{n}$ are called the response variables, which can be
interpreted as noise variables driving the membrane potential $V$ and firing
rate $\dot{n}$, and $\\{\tilde{j},j,\tilde{h},h\\}$ are “source” fields.
Moments of the joint probability distribution can be obtained by taking
derivatives of the MGF with respect to the source fields. Generating
functionals are a powerful tool for studying stochastic dynamics [54, 57, 55,
56]. In neuroscience, these functional methods have previously been applied to
firing rate networks, e.g., [58, 42, 59, 60] and coarse-grained networks of
neural activity, e.g., [61, 62, 63, 64, 65]. Recent work has applied such
methods to spiking networks with “soft resets” that reduce the membrane
potential by fixed amounts after a spike [28, 29, 30]. In this work we follow
[33] in applying this formulation to networks with hard resets.

For random synaptic connections $J_{ij}$ we can average the MGF over
realizations of the synaptic connections to derive an effective dynamical
description of the population dynamics. In the weak coupling regime in which
the weights $J_{ij}$ scale as $1/N$, one can ignore the higher cumulants of
the connectivity matrix and describe the collective dynamics of the network
using only the mean connection strength $J$. If we assume that there are $M$
clusters of neurons that are homogeneous within a cluster, the mean-field
treatment reduces the $N$-dimensional dynamics to an effective $M$-dimensional
description of the population statistics within each cluster corresponding to
the action

$\displaystyle S[\tilde{V},V,\tilde{n},\dot{n}]$
$\displaystyle=\sum_{\mu=1}^{M}\left\\{\tilde{V}_{\mu}(t)\left[\dot{V}_{\mu}(t)+V_{\mu}(t)-\mathcal{E}_{\mu}+\dot{n}_{\mu}(t^{+})V_{\mu}(t)-\sum_{\nu=1}^{M}J_{\mu\nu}\langle\dot{n}_{\nu}(t)\rangle\right]+\tilde{n}_{\mu}(t)\dot{n}_{\mu}(t)-\left(e^{\tilde{n}_{\mu}(t)}-1\right)\phi(V_{\mu}(t))\right\\}.$
(3)

where the Greek indices label the different clusters and
$\langle\dot{n}_{\nu}\rangle$ is the population averaged activity that needs
to be determined self-consistently from the mean-field equations of motion.
This decouples the neurons in the network, such that the average behavior of
an individual neuron in a cluster depends only on the mean synaptic field
$\sum_{\nu=1}^{M}J_{\mu\nu}\langle\dot{n}_{\nu}\rangle$. This description is
exact in the $N\to\infty$ limit [66]. The MGF in (3) can then be used to
obtain the deterministic mean-field approximation for the population membrane
potential and firing rate:

$\displaystyle\frac{d\bar{V}_{\mu}}{dt}$
$\displaystyle=-\bar{V}_{\mu}+\mathcal{E}_{\mu}+\sum_{\nu=1}^{M}J_{\mu\nu}\bar{n}_{\nu}-\bar{n}_{\mu}\bar{V}_{\mu},$
$\displaystyle\bar{n}_{\mu}$ $\displaystyle=\phi(\bar{V}_{\mu}).$ (4)

Unlike standard dynamic mean-field treatments of firing rate networks, Eq.
II.1 is not exact. The stochasticity of spiking remains in Eq. 3 and
influences the dynamics of the population-averaged membrane potentials
$V_{\mu}(t)$ through the spike reset $\dot{n}_{\mu}(t^{+})V_{\mu}(t)$.
Nevertheless, deterministic mean-field approximations often paint a good
qualitative picture of the dynamics of the system, and we will first
investigate the different phases of network activity within this
approximation. Then, in Section II.4, we account for the influence of spiking
fluctuations on network activity by expanding around steady-state solutions of
the mean-field equations to calculate Gaussian (one-loop) corrections.

Alternatively, we can take advantage of the hard reset to estimate the exact
firing rate of the networks using renewal theory. For self-averaging networks,
the density in (3) factorizes over clusters. Each neuron ($\mu\in M$) spikes
independently given a self-consistent mean-field input. When the network
reaches a steady-state, the self-consistent input becomes independent of time,
in which case the spiking dynamics can be treated as a renewal process due to
the hard reset of the membrane potential after spike emission. After each
spike, the membrane potential asymptotically evolves to the net input from the
reset potential (=0 here):

$V_{\mu}(t)=C_{\mu}(1-e^{-t}),$ (5)

where
$C_{\mu}=\mathcal{E}_{\mu}+\sum_{\nu=1}^{M}J_{\mu\nu}\langle\dot{n}_{\nu}\rangle$
is the net input that the neuron receives with $\langle\dot{n}_{\nu}\rangle$
to be determined self-consistently. For an arbitrary nonlinearity $\phi(V(t))$
the inter-spike interval density $p_{\mu}(s)$ is the product of the
instantaneous firing probability at time $s$ and the survival probability
until time $s$:

$p_{\mu}(s)=\phi(V_{\mu}(s))\exp{\Big{(}-\int_{0}^{s}\phi(V_{\mu}(t))dt\Big{)}}$
(6)

The mean inter-spike interval for each cluster is $\langle
s_{\mu}\rangle=\int_{0}^{\infty}ds~{}s\,p_{\mu}(s)$. The self-consistent
firing rate is the inverse of the mean inter-spike interval:
$\langle\dot{n}_{\mu}\rangle=1/\langle s_{\mu}\rangle$, and the roots of this
system of equations allow us to estimate the steady-state firing rates for the
network.

Figure 2: a) Schematic for the homogeneous network. The network has a coupling
strength $J$ and receives an external current $\mathcal{E}$. b) Mean-field
(MF) phase diagram for the homogeneous network with $\alpha=1$ in the
$\mathcal{E}{-}J$ plane, separating the bistable (B–H, Q) high firing rate and
quiescent region from the monostable high (H) firing rate and quiescent (Q)
regions. c) Same as b for $\alpha=2$ separating the bistable (B–H, Q) high
firing rate and quiescent, bistable (B–H, L) high and low firing rate regions
from the monostable high (H) firing rate and quiescent (Q) regions. d) Same as
c for $\alpha=3$. e) Raster plot for the stochastic spiking network
($\alpha=2$) at the parameter marked with cross in panel c ($J=3.2$,
$\mathcal{E}=1.05$), illustrating input driven transition between the two
active states. f) Raster plot for the stochastic spiking network ($\alpha=2$)
at the parameter marked with square in panel c ($J=3.0$, $\mathcal{E}=1.07$),
illustrating stochastic transition between the two active states. g) MF
(green), 1-loop (cyan) and renewal theory (dark blue) firing rate predictions
compared to simulations (brown) for fixed $J=3.0$ ($\alpha=2$). Inset
highlights the low firing rate state in the black square. h) MF phase diagram
that depicts that the critical point continuously increases with the exponent
of the intensity and reaches a limiting value for arbitrarily large $\alpha$.
Violet: $\alpha=1$, Blue: $\alpha=2$, Yellow: $\alpha=3$, Red:
$\alpha\rightarrow\infty$.

###  II.2 Homogeneous Networks

We’re interested in examining the steady states of the mean-field equations
(II.1). For simplicity, we start with a homogeneous excitatory network of $N$
coupled neurons (Fig. 2a). In this case, the mean-field theory is one-
dimensional and a phase diagram can be obtained by investigating the fixed
points of the mean-field dynamics. The phase diagram will be a function of the
external input ($\mathcal{E}$) and the mean coupling strength ($J$).

####  II.2.1 Threshold Power-Law Intensity

We first investigate population activity of threshold power-law networks. The
mean-field equation of motion simplifies to

$\frac{d\bar{V}}{dt}=-\bar{V}+\mathcal{E}+(J-\bar{V})\lfloor\bar{V}-1\rfloor_{+}^{\alpha},$
(7)

where we have fixed the threshold at $\theta=1$ without loss of generality.
The intensity function is zero below threshold so $V=\mathcal{E}$ is the only
possible sub-threshold solution, which does not exist if $\mathcal{E}>1$.
Above threshold, one can show that for any value of the exponent $\alpha$, the
above equation can have a maximum of 3 steady-state solutions, giving rise to
bistability in the mean-field dynamics (see Appendix B). Bistability in the
mean-field dynamics corresponds to metastability in the stochastic spiking
network dynamics. It was previously shown in Ref [33], that with a threshold
linear intensity function ($\alpha=1$), the homogeneous network can either be
monostable or metastable (Fig. 2b). It is not possible to explicitly calculate
the full phase diagram solution for arbitrary $\alpha$, but the cases
$\alpha=2,3$ are analytically tractable and we give the full solutions here.
For $\alpha>3$ we can still characterize several general properties of the
phase diagram.

We begin with the case $\alpha=2$. The mean-field equation of motion (II.1)
for the threshold-quadratic intensity function is a cubic equation (assuming
$\bar{V}$ is above threshold):

$\displaystyle\bar{V}^{3}-\bar{V}^{2}(J+2)+\bar{V}(2J+2)-(J+\mathcal{E})$
$\displaystyle=0.$ (8)

While one can solve for the roots of this equation, we are more interested in
the phase diagram of the network in the $\mathcal{E}{-}J$ plane, which
motivates us to find the boundaries between cases for which we have multiple
real solutions above threshold. We can identify the phase boundaries by
analyzing the discriminant, which determines the boundary between a single
real solution and three real solutions of Eq. (8). The discriminant is a
quadratic in $\mathcal{E}$ whose roots are functions of $J$, yielding the
mean-field boundaries between bistable and monostable phases; see Appendix B
for details. Analogous to the threshold-linear intensity function, the mean-
field dynamics are bistable for sufficiently large values of $J$. However, in
threshold-linear networks the bistability is only between a quiescent state
and an active state that exists only for subthreshold input $\mathcal{E}$. The
superlinear networks retain this regime, but can also exhibit bistability for
superthreshold input, and in this region, the bistable state has two active
steady-state solutions (Fig. 2c). This allows for stochastic switching between
two corresponding metastable active states in the full stochastic network
(Fig. 2f).

Turning to $\alpha=3$, the mean-field equation of motion (II.1) for the
threshold-cubic intensity function is a quartic polynomial when $\bar{V}$ is
above threshold:

$\displaystyle\bar{V}^{4}+b\bar{V}^{3}-c\bar{V}^{2}+d\bar{V}+f$
$\displaystyle=0$ (9)

where the coefficients $\\{b,c,d,f\\}$ are
$b=-(J+3),~{}c=3(J+1),~{}d=-3J,~{}f=(J+\mathcal{E})$. As in the previous case,
we are interested in identifying phase boundaries, which we can again identify
by investigating the discriminant of the quartic equation, this time a cubic
in $\mathcal{E}$; we give the details in Appendix B. We again find that only
two stable fixed points exist above threshold, meaning the bistable state can
have two active steady-state solutions for superthreshold input (Fig. 2d).
Furthermore, we observe that increasing the exponent $\alpha$ from $2$ to $3$
expands the region in which both the bistable states are active (Fig. 2c vs
d).

For the homogeneous network, the membrane dynamics can be treated as a renewal
process where after each spike the membrane potential evolves as:
$V(t)=C(1-e^{-t})$ with $C=\mathcal{E}+J\langle\dot{n}\rangle$. The mean
firing rate can then be estimated by numerically finding the roots of (see
Appendix C for details):

$\langle\dot{n}\rangle=\Bigg{(}\int_{\ln{\frac{C}{C-1}}}^{\infty}ds~{}s\big{(}C(1-e^{-s})-1\big{)}^{\alpha}e^{-I_{p}(s|C)}\Bigg{)}^{-1}$
(10)

where

$\displaystyle I_{p}(s|C)=$
$\displaystyle\frac{(C(1-e^{-s})-1)^{\alpha+1}}{(C-1)(\alpha+1)}$
$\displaystyle\times\;{}_{2}F_{1}\Bigg{(}1,\alpha+1,\alpha+2,\frac{(C(1-e^{-s})-1)}{(C-1)}\Bigg{)}$

and ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function. (The subscript $p$
refers to the power-law intensity function.) Fig. 2g shows the firing rate
predictions of the mean-field theory, the one-loop correction (discussed in
Section II.4), and renewal theory compared to simulations.

While analytic results are difficult to obtain for $\alpha>3$, we can show
that there are at most two stable fixed points above threshold, meaning the
mean-field dynamics can be at most bistable (Appendix B). We can also
determine the coordinates of the ‘cusp’ of the phase diagram for arbitrary
$\alpha$, finding that it traces a non-monotonic path through the
$J{-}\mathcal{E}$ plane (Fig. 2h). Our result suggests that for $\alpha\gtrsim
2$ an increase in exponent of the intensity function expands the region of
bistability where two active metastable steady-states can coexist.

In all cases the bistable region of the $(J,\mathcal{E})$ plane is subdivided
in two. Below the threshold $\mathcal{E}=1$, there exists a quiescent steady
state and a high-rate steady state. Above threshold, the quiescent state
becomes a low-rate state. Both active states are metastable in the full
stochastic network, which can stochastically transition between them (Fig.
2f). If the states are far apart a transition can be induced with the
application of an external current (Fig. 2e).

####  II.2.2 Exponential intensity functions

For an exponential intensity function, $\phi(V)=e^{V-\theta}$, we find that
the mean-field dynamics can be bistable between two states or monostable,
similar to the threshold power law networks. However, the exponential
nonlinearity is never zero for any finite input so the network does not admit
a quiescent state. Instead, in the monostable regime the firing rate changes
continuously as the parameters $J$ or $\mathcal{E}$ are tuned. We show in
Appendix D that the mean-field dynamics of the homogeneous exponential network
are bistable if:

$\displaystyle J$ $\displaystyle>\theta+2,$ (11) $\displaystyle\mathcal{E}$
$\displaystyle>J-(1-W_{0})\left[1+\exp\left(J-1-\theta+W_{0}\right)\right],$
$\displaystyle\mathcal{E}$
$\displaystyle<J-(1-W_{-1})\left[1+\exp\left(J-1-\theta+W_{-1}\right)\right],$

where $W_{0,-1}\equiv W_{0,-1}(-e^{\theta+1-J})$ is the Lambert W function.
Figure 3a shows the mean-field bistability boundaries for various values of
$\theta\in[0,5]$. We see these states in the stochastic spiking network as
well, where the network can either be monostable (Fig. 3b) or metastable (Fig.
3c).

We can also estimate the exact firing rates for the exponential intensity
using renewal theory. In this case, the self-consistent firing rates are roots
of (see Appendix C for details):

$\langle\dot{n}\rangle=\Bigg{(}\int_{0}^{\infty}ds\,s\,\exp\left[C(1-e^{-t})-\theta-
I_{e}(s|C)\right]\Bigg{)}^{-1},$ (12)

where $I_{e}(s|C)=\exp\left(C-\theta\right)[-{\rm Ei}(-Ce^{-s})+{\rm Ei}(-C)]$
and ${\rm Ei}(x)=-\int_{-x}^{\infty}dt~{}e^{-t}/t$ is the exponential
integral. Figure 3d shows the firing rate predictions for a homogeneous
network with an exponential intensity function.

Figure 3: a) MF phase diagram for the exponential intensity function
separating the bistable (B) regime from the monostable (M) regime for various
values of $\theta\in[0,5]$ b) Raster plot for the stochastic spiking network
($\theta=1$) at the parameter marked with square in panel a ($J=4.0$,
$\mathcal{E}=-0.75$), illustrating the monostable states. c) Raster plot for
the stochastic spiking network ($\theta=1$) at the parameter marked with cross
in panel a ($J=4.0$, $\mathcal{E}=-2.0$), illustrating input driven transition
between the the two active states. d) MF (green), 1-loop (cyan) and renewal
theory (dark blue) firing rate predictions compared to simulations (brown) for
fixed $J=6.0$ ($\theta=3$). Inset highlights the low firing rate state in the
main plot.

###  II.3 Excitatory-Inhibitory Networks

Neuronal networks in the brain are not homogeneous. Realistic networks follow
Dale’s law: there are both excitatory and inhibitory neurons in these
networks. To look at more biologically realistic networks, here we consider an
excitatory-inhibitory (EI) network with pulse coupling (Fig. 4a). We introduce
an additional parameter ($g$) that quantifies the relative strength of
inhibition such that the mean connectivity matrix has the following structure:

$\begin{bmatrix}J_{EE}&J_{EI}\\\ J_{IE}&J_{II}\\\
\end{bmatrix}=\begin{bmatrix}J&-gJ\\\ J&-gJ\\\ \end{bmatrix}$ (13)

If both populations receive the same input, the mean-field theory for this EI
network reduces to the one-dimensional model with the replacement
$J\rightarrow J(1-g)$ in the homogeneous network results.

In Fig. 5 we show the results for threshold power law intensities in the
$\mathcal{E}{-}g$ plane (Fig. 4b), along with a raster plot of a simulation to
confirm the full stochastic network exhibits metastability (Fig. 4c). We
compare the mean-field predictions of the firing rates with the one-loop
corrections (discussed in the next section), and the estimates from renewal
theory and the simulations (Fig. 4d).

In Fig. 6 we show the corresponding results for E-I networks with exponential
nonlinearity, this time in the $J{-}g$ plane (Fig. 5b), along with
demonstration of input-driven metastability (Fig. 5c) and the comparisons of
mean-field, one-loop, renewal theory, and simulation estimates of the firing
rates (Fig. 5d).

Figure 4: a) Schematic for the excitatory-inhibitory (EI) network (E: Orange,
I: Blue). The excitatory cluster has a coupling strength $J$, the relative
strength of inhibition is $g$, and both populations receive an external
current $\mathcal{E}$. b) MF phase diagram for threshold power law intensity
in the $\mathcal{E}{-}g$ plane for $\alpha=1$ (dashed), $\alpha=2$ (dot
dashed) and $\alpha=3$ (solid). The phase are same as labelled in Fig. 2. c)
Raster plot for the stochastic spiking network ($\alpha=2$) at $J=5.0$,
$g=0.4$, and $\mathcal{E}=1.09$ illustrating stochastic transition between the
two active states. d) MF (green), 1-loop (cyan) and renewal theory (dark blue)
firing rate predictions compared to simulations (brown) for fixed $J=5.0$,
$g=0.4$ ($\alpha=3$). Inset highlights the low firing rate state in the main
plot. Figure 5: a) MF phase diagram for the exponential intensity in the
$\mathcal{E}{-}g$ plane for various values of $\theta\in[0,5]$. The phase are
same as labelled in Fig. 3. b) Same as a but in the $J{-}g$ plane for
$\theta=1$. c) Raster plot for the stochastic spiking network ($\theta=1$) at
the parameter marked with cross in panel a ($J=8.0$, $g=0.4$, and
$\mathcal{E}=-2.5$), illustrating input driven transition between the the two
active states. d) MF (green), 1-loop (cyan) and renewal theory (dark blue)
firing rate predictions compared to simulations (brown) for fixed $J=8.0$ and
$g=0.4$ ($\theta=1$). The inset shows the low firing rate state in the main
figure.

###  II.4 Impact of fluctuations on the mean activity

While the mean-field theory provides a reasonably accurate qualitative
description of the network dynamics and phase transitions in the network, it
is quantitatively inaccurate (e.g., Figs. 2g, 3-5d) because it ignores all
fluctuations. These fluctuations impact the mean activity due to i) the
nonlinear spike reset and ii) the nonlinear intensity function. In isolation,
these have opposite effects: the spike reset suppresses activity [33], while a
concave-up intensity function promotes activity [28]. We next ask: in the
presence of both sources of nonlinearities, do fluctuations enhance firing or
suppress it?

To answer this we calculate the so-called perturbative 1-loop corrections to
the mean-field predictions [33]. In the path integral formalism these
corrections can be evaluated by perturbatively expanding the MGF under a
Gaussian approximation. This can be achieved by the use of a diagrammatic
expansion (Appendix E). These Feynman diagrams are constructed from edges and
vertices.

The edges in the Feynman diagrams correspond to the linear responses of
configuration variables ($V,\dot{n}$) to the response variables
($\tilde{V},\tilde{n}$) – also called propagators. Since we have 2 each of
configuration and response variables, they give rise to 4 propagators in the
model, given here in the frequency domain:

$\displaystyle\Delta_{n\tilde{n}}(\omega)=$
$\displaystyle\frac{1+\bar{n}+i\omega}{1+\bar{n}+\bar{V}\phi^{(1)}(\bar{V})+i\omega}\equiv\vbox{\hbox{
\feynmandiagram[layered layout, horizontal=a to b] { a -- b };}}$ (14)
$\displaystyle\Delta_{V\tilde{n}}(\omega)=$
$\displaystyle\frac{-\bar{V}}{1+\bar{n}+\bar{V}\phi^{(1)}(\bar{V})+i\omega}\equiv\vbox{\hbox{\feynmandiagram[layered
layout, horizontal=a to b] { a -- [photon] b };}}$
$\displaystyle\Delta_{n\tilde{V}}(\omega)=$
$\displaystyle\frac{\phi^{(1)}(\bar{V})}{1+\bar{n}+\bar{V}\phi^{(1)}(\bar{V})+i\omega}\equiv\vbox{\hbox{\feynmandiagram[layered
layout, horizontal=a to b] { a -- [gluon] b };}}$
$\displaystyle\Delta_{V\tilde{V}}(\omega)=$
$\displaystyle\frac{1}{1+\bar{n}+\bar{V}\phi^{(1)}(\bar{V})+i\omega}\equiv\vbox{\hbox{\feynmandiagram[layered
layout, horizontal=a to b] { a -- [scalar] b };}}$

where $\phi^{(n)}(V)$ is the $n$th derivative of the intensity and $\bar{V}$
and $\bar{n}$ are the mean-field estimates of the membrane potential and
firing rate, respectively (Eq. II.1 for a single population).

The vertices in the Feynman diagrams correspond to source and interaction
terms in the action (Appendix E). The source terms correspond to sources of
fluctuations—here, the stochastic spike emission. The interaction terms are
generated by the nonlinearities, i.e., the spike reset and the nonlinear
intensity function. The vertices that give rise to 1-loop corrections are:

$\displaystyle\vbox{\hbox{ \leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman[small] \vertex[dot] (a)
{}; \vertex[above left=of a] (f1) {$\tilde{n}$}; \vertex[below left=of a] (f2)
{$\tilde{n}$}; \diagram* { {(f1), (f2)} -- (a), };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}},\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman[small] \vertex[empty
square] (a) {}; \vertex[left=of a] (b) {$\tilde{V}$}; \vertex[above right=of
a] (f1) {$n$}; \vertex[below right=of a] (f2) {$V$}; \diagram* { (b) -- (a) --
{(f1), (f2)}, }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}},\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman[small] \vertex[empty dot]
(a) {}; \vertex[left=of a] (b) {$\tilde{n}$}; \vertex[above right=of a] (f1)
{$V$}; \vertex[below right=of a] (f2) {$V$}; \diagram* { (b) -- (a) -- {(f1),
(f2)}, }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}$ (15)

where the first vertex is a source vertex that emits 2 response variables
$\tilde{n}$ and has an amplitude $\phi^{(0)}(\bar{V})$, the second vertex
comes from the nonlinear spike reset and has an amplitude $-1$, and the last
vertex comes from the nonlinear intensity function and has an amplitude
$\phi^{(2)}(\bar{V})/2$. The second and third vertices are interaction
vertices. Using the Feynman rules derived by [33], we can evaluate the loop
corrections to $V$ and $\dot{n}$. Both quantities get corrections from the
intensity function as well as the spike reset. The perturbative 1-loop
corrections to the mean-field theory are:

$\displaystyle\langle V\rangle=\bar{V}\quad+\quad\vbox{\hbox{ \leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [scalar] b
[empty square] -- [photon, half left] c [dot] c -- [half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\quad+\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [photon] b
[empty dot] -- [photon, half left] c [dot] c -- [photon, half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}$ (16)

$\displaystyle\langle\dot{n}\rangle=\bar{n}\quad+\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [gluon] b
[empty square] -- [photon, half left] c [dot] c -- [half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\quad+\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- b [empty
dot] -- [photon, half left] c [dot] c -- [photon, half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}$ (17)

which can be evaluated using the residue theorem (see [30, 56, 67] for an
introduction to evaluating frequency integrals):

$\displaystyle\langle V\rangle$
$\displaystyle=\bar{V}-\frac{\bar{V}^{2}\phi^{(0)}\phi^{(1)}}{2(1+\bar{n}+\bar{V}\phi^{(1)})^{2}}-\frac{\bar{V}^{3}\phi^{(0)}\phi^{(2)}}{4(1+\bar{n}+\bar{V}\phi^{(1)})^{2}}$
(18) $\displaystyle\langle\dot{n}\rangle$
$\displaystyle=\bar{n}-\frac{\bar{V}^{2}\phi^{(0)}\phi^{(1)2}}{2(1+\bar{n}+\bar{V}\phi^{(1)})^{2}}+\frac{\bar{V}^{2}(1+\bar{n})\phi^{(0)}\phi^{(2)}}{4(1+\bar{n}+\bar{V}\phi^{(1)})^{2}}.$
(19)

As noted earlier, the corrections due to fluctuations in (18), (19) come from
two terms. The first correction comes due to the spike reset, and always
suppresses both firing rate and membrane potential. The second correction is
due to the nonlinear intensity function. It can either suppress or promote the
firing rate/membrane potential depending on the concavity of $\phi$.
Furthermore, the second correction appears with opposite sign in the two
equations. This is because the linear response of the rate to spike
fluctuations is is positive, while the linear response of the voltage to spike
fluctuations is negative—another manifestation of the spike reset. That term
could suppress firing rates while paradoxically increasing the mean membrane
potential, or vice versa. This raises the question: under what conditions do
fluctuations promote or suppress the mean membrane potential and firing rate
of the network?

####  II.4.1 Fluctuation corrections for threshold power Law intensities

The threshold power law intensity function, $\phi(V)=\lfloor
V-1\rfloor^{\alpha}$, can be sublinear ($\alpha<1$) or superlinear
($\alpha>1$). We find that fluctuations can promote the membrane potential
only when the intensity function is sublinear and promote the firing rate only
when the intensity function is superlinear. The boundary between enhancement
versus suppression occurs when the two fluctuation corrections in (18) and
(19) cancel. With $\alpha>1$ there are two mutually exclusive possibilities:
the corrections to the mean membrane potential vanish or the corrections to
the firing rates vanish. These occur when:

$\begin{cases}\bar{V}=\frac{2\theta}{\alpha+1}&;~{}\delta V=0\\\
\bar{V}=\Big{(}\frac{\alpha-1}{\alpha+1}\Big{)}^{1/\alpha}+\theta&;~{}\delta\dot{n}=0.\end{cases}$
(20)

These curves define the surfaces where the one-loop fluctuation corrections
vanish (Fig. 6a). These can be viewed in the $\mathcal{E}{-}J$ plane by
substituting the above solutions for $\bar{V}$ in the mean-field equations of
motion (II.1):

$\begin{cases}\mathcal{E}+J\Big{(}\frac{\alpha+\theta}{\alpha+1}\Big{)}^{\alpha}-\Big{[}\Big{(}\frac{\alpha+\theta}{\alpha+1}\Big{)}^{\alpha}+1\Big{]}\Big{(}\frac{2\theta}{\alpha+1}\Big{)}=0&;~{}\delta
V=0\\\
\mathcal{E}+J\Big{(}\frac{\alpha-1}{\alpha+1}\Big{)}-\Big{[}\Big{(}\frac{\alpha-1}{\alpha+1}\Big{)}^{1/\alpha}+\theta\Big{]}\Big{(}\frac{2\alpha}{\alpha+1}\Big{)}=0&;~{}\delta\dot{n}=0.\end{cases}$
(21)

For the monostable regimes, these are lines in the $\mathcal{E}{-}J$ plane
that separate the regions where the membrane potential/firing rate are
promoted/suppressed. For the sublinear intensity functions, increasing
$\alpha$ increases the intercept of the lines (Fig. 6b), while for superlinear
intensity functions, increasing $\alpha$ increases the magnitude of the slope
as well as intercept of the lines (Fig. 6c). The equations suggest that for a
stronger coupling and a larger external input, fluctuations always suppress
the firing rate and the membrane potential. This is a result of the nonlinear
spike reset that prevents the dynamics from exploding despite the unbounded
nonlinearity. Fig. 6d shows the difference between the firing rate of the
stochastic spiking network and the mean-field prediction of the firing rate.

Figure 6: a) Phase diagram in the $V{-}\alpha$ plane showing the regions where
fluctuations promote/suppress the mean firing rate and membrane potential. The
curves for $\alpha<1$ and $\alpha>1$ meet at the same point as
$\alpha\rightarrow 1$, where the corrections due to the intensity vanish. b
Phase diagram in the $\mathcal{E}{-}J$ plane depicting the lines where the
fluctuations vanish for sublinear intensity functions ($\alpha\in[0,0.9]$).
Fluctuations promote the membrane potential to the left of each line and
suppress the membrane potential to the right of each line. c) Same as b but
for superlinear intensity functions ($\alpha\in[1.5,5.0]$). Fluctuations
promote the firing rate to the left of each line and suppress the firing rate
to the right of each line. d) The difference between simulation and mean-field
predictions of the firing rate for a homogeneous network with a superlinear
intensity ($\alpha=2$). The solid black line represents the theoretical
predictions for the points where fluctuations vanish. e) Same as c but in the
$J{-}g$ plane for an EI network. Here, fluctuations promote the firing rate to
the right of each line and suppress the firing rate to the left of each line.
f) Same as d but in the $J{-}g$ plane for fixed $\mathcal{E}=1.5$. Figure 7:
a) Phase diagram in the $J{-}g$ plane for fixed $\mathcal{E}=1.5$ depicting
the lines where the fluctuations vanish for exponential intensity functions
($\theta\in[0,5]$). Fluctuations promote the firing rate to the left of each
line and suppress the firing rate to the right of each line. b) The difference
between simulation and mean-field predictions of the firing rate for a
homogeneous network with exponential intensity ($\theta=1$). The solid black
line represents the theoretical predictions for the points where fluctuations
vanish. c) Same as a but in the $J{-}g$ plane for for an EI network. Here,
fluctuations promote the firing rate to the right of each line and suppress
the firing rate to the left of each line. d) Same as b but in the $J{-}g$
plane for fixed $\mathcal{E}=-0.5$.

For the EI network, the same boundaries can be estimated with the replacement
$J\rightarrow J(1-g)$ in (21). The boundaries for superlinear intensity
functions and the comparisons with simulation are shown in Figure 6e,f. The
equations suggest that–similar to the homogeneous case–fluctuations always
suppress the firing rates and membrane potentials for the parameter regimes
where the network activity becomes very large.

####  II.4.2 Fluctuation corrections for exponential intensities

The exponential intensity functions receive the same fluctuations corrections
as in (18) and (19). For an exponential, however, $\phi^{(2)}(V)>0$ for all
$V>-\infty$. Hence, there are no parameter regimes where the membrane
potential can be promoted by the fluctuations, since both correction terms to
$\bar{V}$ come with a negative sign. As before, we assume the intensity
function has the form $\phi(V)=e^{V-\theta}$ where $\theta$ is the soft
threshold of the intensity function. The nullcline where the fluctuations
corrections vanish for an exponential intensity function is:

$\mathcal{E}+J=2\theta$ (22)

which again is a straight line in the $\mathcal{E}{-}J$ plane (Fig. 7a). For
the EI network, the same boundaries can be estimated with the replacement of
$J\rightarrow J(1-g)$ in (22) (Fig. 7b). Fig. 7c, d show the difference
between the firing rates of the stochastic spiking network and the mean-field
prediction for the homogeneous population ($\mathcal{E}{-}J$ plane) and the EI
populations ($J{-}g$ plane) for fixed $\mathcal{E}$, respectively.

##  III Discussion

We analyzed how single neuron intensity functions impact the overall network
activity in a network of stochastic leaky-integrate and fire neurons. Here, we
specifically focused on the threshold-power law intensity functions, which
have been fit well to experimentally recorded activity in the visual cortex
[37, 38, 39, 40], and the exponential intensity functions that have been
extensively used in fitting point process GLMs to neuronal recordings [43, 44,
45, 46, 47, 29, 48, 49, 50, 31]. We demonstrated that networks with
superlinear threshold-power law intensity functions can have a maximum of 2
stable steady states in the weak coupling regime ($J_{ij}\sim 1/N$). Similar
to threshold linear intensity functions [33], the homogeneous and EI networks
can be monostable (active or quiescent), or metastable between a high firing
rate and quiescent state. However, for superlinear intensities a new
metastable region that admits active high- and low-firing rates emerges [68,
11], reminiscent of the up-down transitions seen in numerous cortical areas
[3, 4, 5, 6, 7]. We find that increasing the exponent of the intensity
function continuously expands the region where these active states co-exist.
(Fig. 2). This high-low metastability also appears in networks with
exponential intensity functions (Fig. 3), owing to the fact that the firing
intensity is never zero—the “subthreshold” firing probability is small but
still finite. Depending on the network parameters, the transitions between the
two active states can be driven stochastically, or with the application of an
external current.

While the mean-field theory provides a decent qualitative description of the
network dynamics, it neglects all fluctuations. Using a diagrammatic
expansion, one can add the effects of nonlinear fluctuations on the network
activity. In the sLIF model, the corrections to the mean-field come from the
nonlinear intensity function and the hard reset of the membrane potential
after a spike [33]. The effect of the nonlinear intensity depends on its
concavity. For the threshold power law intensity, when the intensity function
is sublinear ($\alpha<1$) fluctuations always suppress the firing rate and
promote the membrane potential. When the intensity function is superlinear
($\alpha>1$)—thought to be the more relevant case experimentally—fluctuations
always suppress the firing rate and promote the membrane potential. The two
corrections to both variables (firing rate and membrane potential) can vanish
simultaneously only when the intensity is threshold linear $\alpha=1$. For the
exponential intensity, the concavity is always positive, so that the nonlinear
intensity always suppresses the membrane potential and always promotes the
firing rate.

Stochastic transitions in multistable firing states have been ubiquitously
observed in experimental recordings, and are often associated with behavioral
readouts in the tasks [3, 16, 12, 4, 5, 18, 19, 20, 9, 10]. Understanding how
networks parameters can shape such stochastic transitions in multistable
networks of spiking neurons has been a long standing question in neuroscience
[17, 8]. The field theoretic formulation provides a powerful method to
investigate these transitions. Previous work has estimated transition
probabilities in networks of binary units (loosely interpreted as active or
inactive neurons) [21, 62, 69]. Adaptations of this formalism can be
potentially be applied to these networks of spiking neurons to delineate the
underlying mechanisms that drive these stochastic transitions and elucidate
how network parameters shape the transition rates between different states. We
leave this as a direction for future work.

Another avenue for extension of our work is the investigation of networks with
“strong” synaptic weights that scale as $1/\sqrt{N}$ rather than $1/N$. In
strongly coupled enetworks the variance of the inputs to a neuron is
$\mathcal{O}(1)$, and typically balanced in such a way that the mean input
cancels out. Firing rate networks with strong synaptic connections have been
studied extensively [42, 58, 70, 60, 71], but extending these methods to
networks of spiking neurons is non-trivial due to the feedback between the
spiking and membrane potential evolution.

Extending the analysis to a finite-network size and strong synaptic
connections would yield valuable insights about the mechanisms that drive
stochastic transitions in networks of spiking neurons.

###### Acknowledgements.

SP and BAWB thank the National Institute of Mental Health and National
Institute for Neurological Disorders and Stroke grant UF-1NS115779-01 and
Stony Brook University for financial support for this work.

##  Appendix A Fits to data

We use publicly available cell data from the Allen Institute for Brain Science
to fit the nonlinear intensity functions $\phi(V)$ [47]. We use data from 6
different cell types from mouse visual cortex; corresponding Cre-line, area,
cortical layer, and dendritic type are given in Table 1. In the ‘long square’
stimulation protocol, the cells were stimulated at various levels of input
current amplitudes to get their steady-state response. We use the data from
this stimulation protocol to fit the mean firing rate during the stimulation
period as a function of the mean (time-averaged) membrane potential during the
inter-spike intervals (excluding the hyperpolarization period), minimizing the
squared distance between the data and fit. Within the mean-field approximation
we use in this work, a scatterplot of these means is expected to follow
$\phi(V)$. The general form of the intensity we consider for the fits are:

$\displaystyle\begin{cases}\phi(V)=\lfloor\gamma(V-\theta)\rfloor_{+}^{\alpha}&\rm{Power~{}Law}\\\
\phi(V)=\exp{[\gamma(V-\theta)]}&\rm{Exponential}\end{cases}$ (23)

where $\gamma$ is the gain of the intensity. The details for the fit are as
follows:

  1. 1.

For each level of stimulation, we extracted the mean membrane potential of the
cell, excluding the membrane potential values between the time of reaching the
threshold (extracted from the spike feature data) and the time of reaching the
trough (trough time).

  2. 2.

We estimated the firing rate by counting the number of spikes during the
stimulation divided by the total stimulation time.

  3. 3.

For the power law intensity, we fixed the threshold at the maximum value of
the mean membrane potential where the firing rate was zero. We then fit the
gain ($\gamma$) and exponent ($\alpha$).

  4. 4.

For the exponential intensity, we treated both gain ($\gamma$) and threshold
($\theta$) as free parameters in our fit.

Cell Type # | Cre Line | Brain Area | Layer |  Dendritic Type
---|---|---|---|---
1 | Ctgf-T2A-dgCre | Lateral visual | 6a/b | Spiny
2 | Scnn1a-Tg2-Cre | Primary visual | 5 | Spiny
3 | Scnn1a-Tg2-Cre | Primary visual | 4 | Spiny
4 | Vip-IRES-Cre | Primary visual | 6a | All
5 | Ctgf-T2A-dgCre | Primary visual | 6b | All
6 | Sim1-Cre-KJ18 | Lateral visual | 5 | Spiny
Table 1:  Cre Line, brain area, cortical layer, and dendritic type for each
‘Cell Type $\\#$’ label in 1. Dendritic type ‘All’ includes spiny, aspiny and
sparsely spiny types. All cells are from mouse brains.

##  Appendix B Steady-states for Threshold Power Law Intensity Function

Let the right-hand-side of the mean-field equation of motion, assuming a
superthreshold mean membrane potential, be

$f(\bar{V})=-\bar{V}+\mathcal{E}+(J-\bar{V})(\bar{V}-1)^{\alpha}$ (24)

where we have set the threshold to $\theta=1$. The values of $\bar{V}$ where
$f(\bar{V})$ vanishes correspond to the steady states of the mean-field
dynamics. The solutions cannot generally be obtained in closed form, but we
can show that at most two solutions exist that satisfy $\bar{V}\geq 1$. We do
so by investigating the points of inflection of $f(\bar{V})$, looking at the
zeros of the second derivative:

$f^{\prime\prime}(\bar{V})=-(\bar{V}-1)^{\alpha-2}[\bar{V}(-\alpha(\alpha+1))+\alpha(\alpha-1)J+2\alpha]$
(25)

From this we see that there are only two points of inflection: $\bar{V}_{{\rm
inf},1}=1$ or $\bar{V}_{{\rm inf},2}=\frac{(\alpha-1)J+2}{\alpha+1}$. Thus,
there is only one point of inflection above threshold for a general $\alpha$.
We also see that $\lim_{V\to\infty}f(V)\rightarrow-\infty$. If
$\mathcal{E}>1$, the value of the function at $\bar{V}_{{\rm inf},1}$,
$f(\bar{V}_{{\rm inf},1})$ is positive, which implies the function can have a
maximum of 3 roots above threshold. If $\mathcal{E}<1$ then $f(\bar{V}_{{\rm
inf},1})<0$, and the function can have a maximum of 2 roots above threshold.
In both cases ones of those roots is an unstable fixed point of the mean-field
dynamics, and will not be observed in practice. Accordingly, the mean-field
steady-state equation can have a maximum of two stable steady states.

We can also use this to determine the location of the critical point where the
two roots meet for a general $\alpha$. The critical point is the point where
the first three derivatives of the steady-state equation vanish simultaneously
($f(\bar{V})=f^{\prime}(\bar{V})=f^{\prime\prime}(\bar{V})=0$). This gives

$\displaystyle J$
$\displaystyle=1+\Big{(}\frac{\alpha-1}{\alpha+1}\Big{)}^{\frac{1-\alpha}{\alpha}},$
$\displaystyle\mathcal{E}$
$\displaystyle=1+\Big{(}\frac{\alpha-1}{\alpha+1}\Big{)}^{\frac{1+\alpha}{\alpha}}.$
(26)

From this we see that $\lim_{\alpha\rightarrow
1}[J,~{}\mathcal{E}]=[2.0,~{}1.0]$, and
$\lim_{\alpha\rightarrow\infty}[J,~{}\mathcal{E}]=[2.0,~{}2.0]$. The critical
point continuously increases for $\alpha>1$, thus, expanding the region above
threshold where the active-active states can co-exist.

For specific values of $\alpha$ we can obtain exact solutions for the mean-
field phase boundaries. In the case $\alpha=2$ the roots of the fixed point
equation

$\bar{V}^{3}-\bar{V}^{2}(J+2)+\bar{V}(2J+2)-(J+\mathcal{E})=0$

can be shown to have $1$ or $3$ real-valued roots. The transition from $1$ to
$3$ can be determined by investigating the discriminant, a quadratic equation
in $\mathcal{E}$:

$\displaystyle\mathbb{D}_{3}$
$\displaystyle=-27\mathcal{E}^{2}+A(J)\mathcal{E}+B(J)$ (27)

where $A(J)=[32(J+2)(J+1)-4(J+2)^{3}-54J]$ and
$B(J)=[36J(J+1)(J+2)-4J(J+2)^{3}+4(J+2)^{2}(J+1)^{2}-32(J+1)^{3}-27J^{2}$. The
roots of this quadratic equation give us the mean-field boundaries that
separate bistable and monostable phases.

For $\alpha=3$ the root equation is quartic in $\bar{V}$,

$\displaystyle\bar{V}^{4}+b\bar{V}^{3}-c\bar{V}^{2}+d\bar{V}+f$
$\displaystyle=0$ (28)

where the coefficients $\\{b,c,d,f\\}$ are
$b=-(J+3),~{}c=3(J+1),~{}d=-3J,~{}f=(J+\mathcal{E})$. The discriminant is
cubic in $\mathcal{E}$:

$\mathbb{D}_{4}=-256\mathcal{E}^{3}+\mathcal{E}^{2}A(J)-\mathcal{E}B(J)+C(J)$
(29)

where
$A(J)=(768J-192bd-128c^{2}+144b^{2}c-27b^{4}),B(J)=768J^{2}+2J[-192bd-128c^{2}+144b^{2}c-27b^{4}]+[144cd^{2}-6b^{2}d^{2}-80bc^{2}d+16c^{4}+18b^{3}cd-4b^{2}c^{3}],C(J)=256J^{3}+J^{2}[-192bd-128c^{2}+144b^{2}c-27b^{4}]+J[144cd^{2}-6b^{2}d^{2}-80bc^{2}d+16c^{4}+18b^{3}cd-4b^{2}c^{3}]+[18bcd^{3}-27d^{4}-4c^{3}d^{2}-4b^{3}d^{3}+b^{2}c^{2}d^{2}]$.
The roots of this cubic equation, along with the condition that $\bar{V}$ must
be above threshold and the following inequalities:

$\displaystyle 8c-3b^{2}$ $\displaystyle<0,$ $\displaystyle
64f-16c^{2}+16b^{2}c-16bd-3b^{4}$ $\displaystyle<0,$ (30)

define the boundaries separating the bistable and the monostable phases.

##  Appendix C Exact firing rates from renewal theory

After averaging over the synaptic weights $J_{ij}$ and neglecting variance in
the weights, which is negligible when $N$ is large, we obtain a population-
averaged stochastic differential equation that follows from the action Eq.
(3):

$\displaystyle\dot{V}_{\mu}$
$\displaystyle=-V_{\mu}(t)+\mathcal{E}_{\mu}-\dot{n}_{\mu}(t^{+})V_{\mu}(t)+\sum_{\nu=1}^{N}J_{\mu\nu}\left\langle\dot{n}_{\nu}(t)\right\rangle,$
(31) $\displaystyle\dot{n}_{\mu}(t)dt$ $\displaystyle\sim{\rm
Poiss}\left[\phi(V_{\mu}(t))dt\right],$ (32)

where $\langle\dot{n}_{\mu}(t)\rangle$ is the population-averaged input to
each cluster. In the large $N$ limit this concentrates to a constant value
that must be determined self-consistently. In this limit the clusters are
formally decoupled and interact only through the self-consistent input.

Due to the hard reset of the membrane potential after a spike, each cluster’s
spike train can be treated as a renewal process. Between spikes the membrane
dynamics of each cluster evolves independently according to:

$V_{\mu}(t)=C_{\mu}(1-e^{-t})$ (33)

where
$C_{\mu}=\mathcal{E}_{\mu}+\sum_{\nu=1}^{M}J_{\mu\nu}\langle\dot{n}_{\nu}\rangle$
is the net input that drives the membrane potential of the cluster $\mu$. The
instantaneous firing intensity, from the model definition is simply
$\phi(V_{\mu}(t))$. The inter-spike interval density is then the product of
the instantaneous firing intensity and the survival probability independently
for each cluster:

$p_{\mu}(s)=\phi(V_{\mu}(s))\exp{\Bigg{(}-\int_{0}^{s}\phi(V_{\mu}(t))dt\Bigg{)}}.$
(34)

The mean inter-spike interval can then be evaluated as $\langle
s_{\mu}\rangle=\int_{0}^{\infty}ds~{}sp_{\mu}(s)$, and the mean firing rate is
the inverse of the mean inter-spike interval
$\langle\dot{n}_{\mu}\rangle=1/\langle s_{\mu}\rangle$. For the nonlinearities
we consider here, it is not possible to get a closed form expression for the
density. However, we can take a semi-analytic approach to estimate the exact
firing rates for the two nonlinearities using renewal theory. In the following
we drop the $\mu$ subscripts in the intermediate results for simplicity, as
different clusters interact only through the coefficients $C$.

###  C.1 Threshold Power Law Intensity

The inter-spike interval density for the threshold power law intensity
function is:

$p(s)=\begin{cases}0&;~{}s\leq\ln{\frac{C}{C-\theta}}\\\
\big{(}C(1-e^{-s})-\theta\big{)}^{\alpha}\\\
\times\exp{\Big{[}-\int_{0}^{s}\big{(}C(1-e^{-t})-\theta\big{)}^{\alpha}dt\Big{]}}&;~{}s>\ln{\frac{C}{C-\theta}}.\end{cases}$
(35)

While it is not possible to get a closed form expression for this integral,
the integral in the exponential of Eqn. (35) can be evaluated in closed form:

$I_{p}(s|C)\equiv\int_{0}^{s}\big{(}C(1-e^{-t})-\theta\big{)}^{\alpha}dt$

To get a closed form expression for this integral, we make a change of
variables to $y=C(1-e^{-t})-\theta$, after which the integral can be
identified as proportional to a representation of the hypergeometric function:

$\displaystyle I_{p}(s|C)$
$\displaystyle=\frac{(C(1-e^{-s})-\theta)^{\alpha+1}}{(C-\theta)(\alpha+1)}$
$\displaystyle\times\;{}_{2}F_{1}\Big{(}1,\alpha+1,\alpha+2,\frac{C(1-e^{-s})-\theta}{C-\theta}\Big{)}.$
(36)

The mean firing rate can then be estimated by numerically finding the roots of
the following system of equations, restoring the cluster subscripts and
recalling that
$C_{\mu}=\mathcal{E}_{\mu}+\sum_{\nu=1}^{M}J_{\mu\nu}\langle\dot{n}_{\nu}\rangle$:

$\langle\dot{n}_{\mu}\rangle=\Bigg{(}\int_{\ln{\frac{C_{\mu}}{C_{\mu}-\theta}}}^{\infty}ds~{}s\big{(}C_{\mu}(1-e^{-s})-\theta\big{)}^{\alpha}e^{-I_{p}(s|C_{\mu})}\Bigg{)}^{-1}.$
(37)

###  C.2 Exponential Intensity

The inter-spike interval density for the exponential intensity function is:

$p(s)=\begin{cases}0&;~{}s\leq 0\\\ e^{(C(1-e^{-s})-\theta)}\times\\\
\exp{\Big{(}-\int_{0}^{s}e^{(C(1-e^{-t})-\theta)}dt\Big{)}}&;~{}s>0.\end{cases}$
(38)

Again, the integral in the exponential of Eqn. (38) can be evaluated in closed
form by making a change of variables $y=Ce^{-t}$, identifying the result in
terms of the “Exponential Integral” special function:

$I_{e}(s|C)=e^{(C-\theta)}[-{\rm Ei}(-Ce^{-s})+{\rm Ei}(-C)].$ (39)

where ${\rm Ei}(x)=-\int_{-x}^{\infty}dt~{}e^{-t}/t$ is the exponential
integral. The mean firing rate can then be estimated by numerically finding
the roots of the following system of equations (after restoring the explicit
subscipts):

$\langle\dot{n}_{\mu}\rangle=\Bigg{(}\int_{0}^{\infty}ds~{}se^{(C_{\mu}(1-e^{-t})-\theta)}e^{-I_{e}(s|C_{\mu})}\Bigg{)}^{-1}.$
(40)

##  Appendix D Steady-states for exponential intensity functions

The steady-state equation for the exponential intensity function is:

$f(\bar{V})=-\bar{V}+\mathcal{E}+(J-\bar{V})e^{\bar{V}-\theta}$ (41)

We can find the extremum of the function by demanding that its derivative
vanishes:

$\displaystyle f^{\prime}(V)$
$\displaystyle=-1-e^{\bar{V}-\theta}+(J-V)e^{\bar{V}-\theta}$ (42)
$\displaystyle\bar{V}_{0,-1}$
$\displaystyle=(J-1)+W_{0,-1}\left(-e^{\theta+1-J}\right)$ (43)

where, $W_{0,-1}$ are the two real branches of the Lambert W function. Thus,
the derivative can vanish at either 2 points, or does not vanish at all. Since
the argument is the negative of an exponential, the derivative can vanish at
exactly 2 points if $-1/e<-e^{\theta+1-J}<0$; otherwise the function is
strictly decreasing. Therefore, the necessary condition for the derivative to
vanish at 2 points is:

$J>\theta+2$ (44)

Since the function can have either 2 extrema or no extrema at all, the mean-
field equations can have 3 real roots if the values of the function $f(V)$ at
$\bar{V}_{0,-1}$ have opposite signs, and will have only 1 real root if the
values of the function $f(V)$ at $\bar{V}_{0,-1}$ have the same sign or if
$f^{\prime}(\bar{V})$ never vanishes ($J<\theta+2$). Substituting the solution
(43) in (41):

$f(\bar{V}_{0,-1})=(\mathcal{E}-J)+\big{(}1-W_{0,-1}\big{)}\big{(}1+e^{J-1-\theta+W_{0,-1}}\big{)}$
(45)

where we have suppressed the argument of Lambert W function. Thus, there exist
3 steady-states in the network if $f(\bar{V}_{0})>0$ and $f(\bar{V}_{-1})<0$,
or:

$\displaystyle\mathcal{E}$ $\displaystyle>J-(1-W_{0})(1+e^{J-1-\theta+W_{0}})$
$\displaystyle\mathcal{E}$
$\displaystyle<J-(1-W_{-1})(1+e^{J-1-\theta+W_{-1}}).$

The mean-field dynamics are bistable between the two curves above, and
monostable otherwise.

##  Appendix E Feynman rules for fluctuation corrections

If the fluctuations are weak, we can expand the MGF in (3) in a perturbative
series that adds corrections to the mean-field predictions. This can be
formally achieved by splitting the action into two parts: i) the free part
($S_{F}$) with terms that are up to bi-linear in the fields
$[\tilde{V},V,\tilde{n},\dot{n}]$; ii) the interaction part ($S_{I}$)
containing the terms with higher powers of the corresponding fields. Splitting
the fields in a background mean and fluctuations ($V\rightarrow\bar{V}+\delta
V,\dot{n}\rightarrow\bar{n}+\delta\dot{n}$), we can find the free and
interacting part of the action in the MGF for the homogeneous case in (3):

$\displaystyle S$ $\displaystyle=S_{F}+S_{I},$ $\displaystyle S_{F}$
$\displaystyle=\tilde{V}[\partial_{t}{\bar{V}}+\bar{V}-\mathcal{E}+J\langle\dot{n}\rangle-\bar{n}\bar{V}]+\tilde{n}[\bar{n}-\phi^{(0)}(\bar{V})],$
$\displaystyle\tilde{V}[\partial_{t}+1+\bar{n}]\delta
V+\tilde{V}[\bar{V}]\delta\dot{n}+\tilde{n}\delta\dot{n}-\tilde{n}[\phi^{(1)}(\bar{V})]\delta
V,$ $\displaystyle S_{I}$ $\displaystyle=\tilde{V}\,\delta\dot{n}\,\delta
V-\sum\limits_{\begin{subarray}{c}p=q=0\\\ \setminus\\{p=1,q=0\\}\\\
\setminus\\{p=1,q=1\\}\end{subarray}}^{\infty}\frac{\tilde{n}^{p}}{p!}~{}\frac{\phi^{(q)}(\bar{V})}{q!}(\delta
V)^{q}.$ (46)

The MGF can now be expanded around the solution to the mean-field theory in
the first line of the free action above:

$\displaystyle\frac{d\bar{V}}{dt}$
$\displaystyle=-\bar{V}+\mathcal{E}+J\langle\dot{n}\rangle-\bar{n}\bar{V},$
$\displaystyle\bar{n}$ $\displaystyle=\phi^{(0)}(\bar{V}).$

The perturbative expansion uses the fact that the case of a Gaussian action
($S_{F}$) can be solved exactly. The interacting action ($S_{I}$) is then
expanded in a functional Taylor series around this Gaussian solution to
systematically add the effect of non-Gaussian features of the distribution. An
arbitrary moment can be evaluated using:

$\displaystyle\langle V\rangle$ $\displaystyle=\bar{V}+\langle\delta
V\rangle,$ $\displaystyle\langle\dot{n}\rangle$
$\displaystyle=\bar{n}+\langle\delta\dot{n}\rangle,$
$\displaystyle\Big{\langle}\prod_{i=1}^{a}\delta\dot{n}(t_{i})\prod_{k=1}^{b}\delta
V(t_{k})\Big{\rangle}$
$\displaystyle=\prod_{i=1}^{a}\prod_{k=1}^{b}\frac{1}{Z[0]}\frac{\delta}{\delta\tilde{h}(t_{i})}\frac{\delta}{\delta\tilde{j}(t_{k})}Z[\tilde{j},j,\tilde{h},h]$
(47)

where $Z[0]$ is the MGF evaluated at vanishing source terms, equal to $1$
because of normalization of the probability density. At the lowest order in
the perturbation, the expectation value of the fluctuations is zero
$(\langle\delta V\rangle=\langle\delta\dot{n}\rangle=0)$, so that $(\langle
V\rangle=\bar{V},\langle\dot{n}\rangle=\bar{n})$. This is the mean-field
solution. The loop expansion gives these fluctuations a non-zero expectation
which provides fluctuation corrections to the mean-field solution. The MGF can
now be expanded around the free part of the action to calculate expectations
beyond the lowest order:

$\displaystyle Z[\tilde{j},j,\tilde{h},h]$
$\displaystyle=\int\mathcal{D}\tilde{V}\mathcal{D}\delta
V\mathcal{D}\tilde{n}\mathcal{D}\delta\dot{n}~{}e^{-S_{F}-S_{I}+\tilde{j}\cdot
V+j\cdot\tilde{V}+\tilde{h}\cdot\dot{n}+h\cdot\tilde{n}},$ $\displaystyle
Z[\tilde{j},j,\tilde{h},h]$
$\displaystyle=\int\mathcal{D}\tilde{V}\mathcal{D}\delta
V\mathcal{D}\tilde{n}\mathcal{D}\delta\dot{n}~{}e^{-S_{F}}\Bigg{[}1$
$\displaystyle+\Big{(}-S_{I}+\tilde{j}\cdot
V+j\cdot\tilde{V}+\tilde{h}\cdot\dot{n}+h\cdot\tilde{n}\Big{)}$
$\displaystyle+\frac{1}{2!}\Big{(}-S_{I}+\tilde{j}\cdot
V+j\cdot\tilde{V}+\tilde{h}\cdot\dot{n}+h\cdot\tilde{n}\Big{)}^{2}$
$\displaystyle+\frac{1}{3!}\Big{(}-S_{I}+\tilde{j}\cdot
V+j\cdot\tilde{V}+\tilde{h}\cdot\dot{n}+h\cdot\tilde{n}\Big{)}^{3}+\dots\Bigg{]}$
(48)

This expansion can be organized diagrammatically using Feynman rules to
estimate various moments of a given distribution in terms of the propagators
evaluated at the mean-field solution. Using Wick’s theorem, any expectation
can be broken into a sum over product of propagators $\langle
x(t)\tilde{x}(t^{\prime})\rangle\rightarrow[\Delta_{n\tilde{n}},\Delta_{V\tilde{n}},\Delta_{n\tilde{V}},\Delta_{V\tilde{V}}]$
(where $x=V$ or $\dot{n}$ and $\tilde{x}=\tilde{V}$ or $\tilde{n}$). The Ito
condition $\langle x(t)\tilde{x}(t)\rangle=0$ ensures that the propagators are
purely causal in the time domain. Thus, the only terms that survive in the
above expansion are the terms that have an even number of response
$[\tilde{V},\tilde{n}]$ and field $[V,\dot{n}]$ variables, and where each
variable in these pairs has a different time argument (comes from different
term in the expansion).

Vertex | Factor | In-degree $(\delta V,\delta\dot{n})$ |  Out-degree $(\tilde{V},\tilde{n})$
---|---|---|---
| $\phi^{(0)}(\bar{V})$ | $(0,0)$ |  $(0,\geq 2)$
| -1 | $(1,1)$ |  $(1,0)$
| $\frac{\phi^{(q)}(\bar{V})}{q!}$ | $(\geq 1,0)$ |  $(0,\geq q)$
Table 2:  The vertices in Feynman diagrams correspond to the terms in the interacting part of the action $(S_{I})$. Each vertex can have incoming field variables (in-degree) or outgoing response variables (out-degree). Propagator | Edge |  Factor $(t,s)$
---|---|---
$\Delta_{n\tilde{n}}(t,s)$ |  [layered layout, horizontal=a to b] a – b ; |  $\delta(t-s)$ $-\bar{V}\phi^{(1)}e^{-(1+\bar{n}+\bar{V}\phi^{(1)})(t-s)}\Theta(t-s)$
$\Delta_{n\tilde{V}}(t,s)$ |  [layered layout, horizontal=a to b] a – [gluon] b ; |  $\phi^{(1)}e^{-(1+\bar{n}+\bar{V}\phi^{(1)})(t-s)}\Theta(t-s)$
$\Delta_{V\tilde{V}}(t,s)$ |  [layered layout, horizontal=a to b] a – [scalar] b ; |  $e^{-(1+\bar{n}+\bar{V}\phi^{(1)})(t-s)}\Theta(t-s)$
$\Delta_{V\tilde{n}}(t,s)$ |  [layered layout, horizontal=a to b] a – [photon] b ; |  $-\bar{V}e^{-(1+\bar{n}+\bar{V}\phi^{(1)})(t-s)}\Theta(t-s)$
Table 3:  The edges in the Feynman diagrams correspond to the propagators
evaluated from the bilinear part of the free action $(S_{F})$. Each propagator
measures the linear response of a field variable to perturbations in the other
variable, and are causal due to the Ito condition. Here $\Theta(t-s)$ is the
Heaviside step function.

The diagrams are constructed out of vertices and edges. The vertices (Table 2)
are the terms in the interacting part of the action $(S_{I})$, and the edges
(Table 3) are constructed out of the propagators evaluated at the solution to
the mean-field theory. The diagrammatic rules to evaluate an arbitrary
expectation $\langle\delta\dot{n}(t_{1})\dots\delta\dot{n}(t_{a})\delta
V(t_{a+1})\dots\delta V(t_{a+b})\rangle$ are as follows:

  1. 1.

Place $a+b$ external vertices corresponding to the expectation value that
needs to be evaluated.

  2. 2.

Use the vertices and edges in Tables 2, 3 to construct all possible connected
graphs such that the external vertices have only one incoming propagator.

  3. 3.

The analytic expression for the diagram can be obtained by multiplying all the
factors of edges and vertices together. In the frequency domain, each edge is
assigned its own frequency variable and the expression can be evaluated by
integrating over all internal frequencies. In the time domain, each vertex
carries its own time argument, and the propagators come with the time
arguments of the vertices they connect. The expression can be evaluated by
again integrating over all internal time points.

  4. 4.

The resulting expectation is the sum of all such diagrams that can be
constructed in (2).

For the action in (46), the vertices and the resulting 1-loop diagrams are
(diagrams generated with [72]):

$\displaystyle\langle\delta V\rangle=\quad\vbox{\hbox{ \leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [scalar] b
[empty square] -- [photon, half left] c [dot] c -- [half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\quad+\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [photon] b
[empty dot] -- [photon, half left] c [dot] c -- [photon, half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}$ (49)

$\displaystyle\langle\delta\dot{n}\rangle=\quad\vbox{\hbox{ \leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- [gluon] b
[empty square] -- [photon, half left] c [dot] c -- [half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\quad+\quad\vbox{\hbox{
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a);
\vertex[left=of a] (b); \vertex[right=of a] (c); \diagram* { a -- b [empty
dot] -- [photon, half left] c [dot] c -- [photon, half left] b };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}$ (50)

The time arguments corresponding to each vertex are: (  \- $t_{1}$), ( ,  \-
$t_{2}$). The above diagrams in terms of the propagators and edges are:

$\displaystyle\langle\delta V\rangle$
$\displaystyle=-\phi^{(0)}\int_{-\infty}^{t}dt_{1}dt_{2}\Delta_{V\tilde{V}}(t,t_{1})\Delta_{V\tilde{n}}(t_{1},t_{2})\Delta_{n\tilde{n}}(t_{1},t_{2})$
$\displaystyle+\frac{\phi^{(0)}\phi^{(2)}}{2}\int_{-\infty}^{t}dt_{1}dt_{2}\Delta_{V\tilde{n}}(t,t_{1})\Delta_{V\tilde{n}}(t_{1},t_{2})\Delta_{V\tilde{n}}(t_{1},t_{2}),$
(51) $\displaystyle\langle\delta\dot{n}\rangle$
$\displaystyle=-\phi^{(0)}\int_{-\infty}^{t}dt_{1}dt_{2}\Delta_{n\tilde{V}}(t,t_{1})\Delta_{V\tilde{n}}(t_{1},t_{2})\Delta_{n\tilde{n}}(t_{1},t_{2})$
$\displaystyle+\frac{\phi^{(0)}\phi^{(2)}}{2}\int_{-\infty}^{t}dt_{1}dt_{2}\Delta_{n\tilde{n}}(t,t_{1})\Delta_{V\tilde{n}}(t_{1},t_{2})\Delta_{V\tilde{n}}(t_{1},t_{2}).$
(52)

These expression can be evaluated either in the time domain or can be
transformed into the frequency domain and evaluated using the residue theorem.
For a detailed introduction to expansions of the moment generating functions
see [67, 56], path integrals specifically applied to stochastic differential
equations see [54], calculations involving loop integrals in the frequency
domain see [30], perturbative expansion for this specific model and its
connections to the effective action see [33].

## References

  * Brinkman  _et al._ [2022] B. A. W. Brinkman, H. Yan, A. Maffei, I. M. Park, A. Fontanini, J. Wang, and G. La Camera, Metastable dynamics of neural circuits and networks, Applied Physics Reviews 9, 011313 (2022), https://pubs.aip.org/aip/apr/article-pdf/doi/10.1063/5.0062603/16495968/011313_1_online.pdf .
  * Rossi  _et al._ [2023] K. L. Rossi, R. C. Budzinski, E. S. Medeiros, B. R. R. Boaretto, L. Muller, and U. Feudel, A unified framework of metastability in neuroscience (2023), arXiv:2305.05328 [q-bio.NC] .
  * Kajikawa  _et al._ [2022] K. Kajikawa, B. K. Hulse, A. G. Siapas, and E. V. Lubenov, Up-down states and ripples differentially modulate membrane potential dynamics across dg, ca3, and ca1 in awake mice, eLife 11, e69596 (2022).
  * Anderson  _et al._ [2000] J. Anderson, I. Lampl, I. Reichova, M. Carandini, and D. Ferster, Stimulus dependence of two-state fluctuations of membrane potential in cat visual cortex, Nat Neurosci 3, 617–621 (2000).
  * Lampl  _et al._ [1999] I. Lampl, I. Reichova, and D. Ferster, Synchronous membrane potential fluctuations in neurons of the cat visual cortex, Neuron 22, 361 (1999).
  * Steriade  _et al._ [2001] M. Steriade, I. Timofeev, and F. Grenier, Natural waking and sleep states: A view from inside neocortical neurons, Journal of Neurophysiology 85, 1969 (2001).
  * Haider  _et al._ [2007] B. Haider, A. Duque, A. R. Hasenstaub, Y. Yu, and D. A. McCormick, Enhancement of visual responsiveness by spontaneous local network activity in vivo, Journal of Neurophysiology 97, 4186 (2007).
  * Schaub  _et al._ [2015] M. T. Schaub, Y. N. Billeh, C. A. Anastassiou, C. Koch, and M. Barahona, Emergence of slow-switching assemblies in structured neuronal networks, PLOS Computational Biology 11, e1004196 (2015).
  * Miller and Katz [2010] P. Miller and D. B. Katz, Stochastic transitions between neural states in taste processing and decision-making, Journal of Neuroscience 30, 2559 (2010), https://www.jneurosci.org/content/30/7/2559.full.pdf .
  * Recanatesi  _et al._ [2022] S. Recanatesi, U. Pereira-Obilinovic, M. Murakami, Z. Mainen, and L. Mazzucato, Metastable attractors explain the variable timing of stable behavioral action sequences, Neuron 110, 139 (2022).
  * Tartaglia and Brunel [2017] E. M. Tartaglia and N. Brunel, Bistability and up/down state alternations in inhibition-dominated randomly connected networks of lif neurons, Sci Rep 7, 11916 (2017).
  * Mejias  _et al._ [2010] J. F. Mejias, H. J. Kappen, and J. J. Torres, Irregular dynamics in up and down cortical states, PLOS ONE 5, e13651 (2010).
  * Mattia and Sanchez-Vives [2012] M. Mattia and M. Sanchez-Vives, Exploring the spectrum of dynamical regimes and timescales in spontaneous cortical activity, Cogn Neurodyn 6, 239–250 (2012).
  * Jercog  _et al._ [2017] D. Jercog, A. Roxin, P. Barthó, A. Luczak, A. Compte, and J. de la Rocha, Up-down cortical dynamics reflect state transitions in a bistable network, eLife 6, e22425 (2017).
  * Kraynyukova and Tchumatchenko [2018] N. Kraynyukova and T. Tchumatchenko, Stabilized supralinear network can give rise to bistable, oscillatory, and persistent activity, Proceedings of the National Academy of Sciences of the United States of America 115, 3464 (2018), https://www.pnas.org/content/pnas/115/13/3464.full.pdf.
  * Holcman and Tsodyks [2006] D. Holcman and M. Tsodyks, The emergence of up and down states in cortical networks, PLoS computational biology 2, e23 (2006).
  * Litwin-Kumar and Doiron [2012] A. Litwin-Kumar and B. Doiron, Slow dynamics and high variability in balanced cortical networks with clustered connections, Nat Neurosci 15, 1498–1505 (2012).
  * Mazzucato  _et al._ [2015] L. Mazzucato, A. Fontanini, and G. L. Camera, Dynamics of multistable states during ongoing and evoked cortical activity, Journal of Neuroscience 35, 8214 (2015), https://www.jneurosci.org/content/35/21/8214.full.pdf .
  * Mazzucato  _et al._ [2019] L. Mazzucato, G. La Camera, and A. Fontanini, Expectation-induced modulation of metastable activity underlies faster coding of sensory stimuli, Nat Neurosci 22, 787–796 (2019).
  * La Camera  _et al._ [2019] G. La Camera, A. Fontanini, and L. Mazzucato, Cortical computations via metastable activity, Current Opinion in Neurobiology 58, 37 (2019), computational Neuroscience.
  * Bressloff [2010a] P. C. Bressloff, Metastable states and quasicycles in a stochastic wilson-cowan model of neuronal population dynamics, Phys. Rev. E 82, 051903 (2010a).
  * Wilson and Cowan [1972] H. Wilson and J. Cowan, Excitatory and inhibitory interactions in localized populations of model neurons, Biophys J. 12, 1–24 (1972).
  * Wilson and Cowan [1973] H. Wilson and J. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Kybernetik 13, 55–80 (1973).
  * Cohen and Grossberg [1983] M. A. Cohen and S. Grossberg, Absolute stability of global pattern formation and parallel memory storage by competitive neural networks, IEEE Transactions on Systems, Man, and Cybernetics SMC-13, 815 (1983).
  * Amari [1977] S. Amari, Dynamics of pattern formation in lateral-inhibition type neural fields, Biol. Cybern. 27, 77–87 (1977).
  * Potthast [2013] R. Potthast, Amari model, in  _Encyclopedia of Computational Neuroscience_, edited by D. Jaeger and R. Jung (Springer New York, New York, NY, 2013) pp. 1–6.
  * Pinto  _et al._ [1996] D. J. Pinto, J. C. Brumberg, D. J. Simons, G. B. Ermentrout, and R. Traub, A quantitative population model of whisker barrels: Re-examining the wilson-cowan equations, J Comput Neurosci 3, 247–264 (1996).
  * Ocker  _et al._ [2017] G. K. Ocker, K. Josić, E. Shea-Brown, and M. A. Buice, Linking structure and activity in nonlinear spiking networks, PLOS Computational Biology 13, e1005583 (2017).
  * Brinkman  _et al._ [2018] B. A. W. Brinkman, F. Rieke, E. Shea-Brown, and M. A. Buice, Predicting how and when hidden neurons skew measured synaptic interactions, PLOS Computational Biology 14, e1006490 (2018).
  * Kordovan and Rotter [2020] M. Kordovan and S. Rotter, Spike train cumulants for linear-nonlinear poisson cascade models (2020), arXiv:2001.05057 [q-bio.NC] .
  * Liang and Brinkman [2023] T. Liang and B. A. W. Brinkman, Statistically inferred neuronal connections in subsampled neural networks strongly correlate with spike train covariance, bioRxiv 10.1101/2023.02.01.526673 (2023).
  * Brinkman [2023] B. A. Brinkman, Non-perturbative renormalization group analysis of nonlinear spiking networks, arXiv preprint arXiv:2301.09600 (2023).
  * Ocker [2023] G. K. Ocker, Republished: Dynamics of stochastic integrate-and-fire networks, Phys. Rev. X 13, 041047 (2023).
  * Brunel  _et al._ [2014] N. Brunel, V. Hakim, and M. J. Richardson, Single neuron dynamics and computation, Current Opinion in Neurobiology 25, 149 (2014), theoretical and computational neuroscience.
  * van Meegen and van Albada [2021] A. van Meegen and S. J. van Albada, Microscopic theory of intrinsic timescales in spiking neural networks, Phys. Rev. Res. 3, 043077 (2021).
  * Ekelmans  _et al._ [2023] P. Ekelmans, N. Kraynyukova, and T. Tchumatchenko, Targeting operational regimes of interest in recurrent neural networks, PLOS Computational Biology 19, e1011097 (2023).
  * Miller and Troyer [2002] K. D. Miller and T. W. Troyer, Neural noise can explain expansive, power-law nonlinearities in neural response functions, Journal of Neurophysiology 87, 653 (2002).
  * Persi  _et al._ [2011] E. Persi, D. Hansel, L. Nowak, P. Barone, and C. van Vreeswijk, Power-law input-output transfer functions explain the contrast-response and tuning properties of neurons in visual cortex, PLOS Computational Biology 7, e1001078 (2011).
  * Nicholas J. Priebe [2008] D. F. Nicholas J. Priebe, Inhibition, spike threshold, and stimulus selectivity in primary visual cortex, Neuron 57, 482 (2008).
  * Hansel and van Vreeswijk [2002] D. Hansel and C. van Vreeswijk, How noise contributes to contrast invariance of orientation tuning in cat visual cortex, Journal of Neuroscience 22, 5118 (2002), https://www.jneurosci.org/content/22/12/5118.full.pdf .
  * Mastrogiuseppe and Ostojic [2017] F. Mastrogiuseppe and S. Ostojic, Intrinsically-generated fluctuating activity in excitatory-inhibitory networks, PLoS computational biology 13, e1005498 (2017).
  * Kadmon and Sompolinsky [2015] J. Kadmon and H. Sompolinsky, Transition to chaos in random neuronal networks, Phys. Rev. X 5, 041030 (2015).
  * Badel  _et al._ [2007] L. Badel, S. Lefort, R. Brette, C. C. H. Petersen, W. Gerstner, and M. J. E. Richardson, Dynamic i-v curves are reliable predictors of naturalistic pyramidal-neuron voltage traces, Journal of Neurophysiology 99, 656 (2007).
  * Pillow  _et al._ [2008] J. W. Pillow, J. Shlens, L. Paninski, A. Sher, A. M. Litke, E. J. Chichilnisky, and E. P. Simoncelli, Spatio-temporal correlations and visual signalling in a complete neuronal population, Nature 454, 995–999 (2008).
  * Mensi  _et al._ [2011] S. Mensi, R. Naud, C. Pozzorini, M. Avermann, C. C. H. Petersen, and W. Gerstner, Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechanisms, Journal of Neurophysiology 107, 1756 (2011).
  * Pozzorini  _et al._ [2015] C. Pozzorini, S. Mensi, O. Hagens, R. Naud, C. Koch, and W. Gerstner, Automated high-throughput characterization of single neurons by means of simplified spiking models, PLOS Computational Biology 11, e1004275 (2015).
  * Teeter  _et al._ [2018] C. Teeter, R. Iyer, V. Menon, N. Gouwens, D. Feng, J. Berg, A. Szafer, N. Cain, H. Zeng, M. Hawrylycz, C. Koch, and S. Mihalas, Generalized leaky integrate-and-fire models classify multiple neuron types, Nature 9, 709 (2018).
  * Simoncelli  _et al._ [2004] E. P. Simoncelli, L. Paninski, J. Pillow, O. Schwartz, _et al._ , Characterization of neural responses with stochastic stimuli, The cognitive neurosciences 3, 1 (2004).
  * Paninski [2004] L. Paninski, Maximum likelihood estimation of cascade point-process neural encoding models, Network: Computation in Neural Systems 15, 243 (2004).
  * Pillow  _et al._ [2005] J. W. Pillow, L. Paninski, V. J. Uzzell, E. P. Simoncelli, and E. J. Chichilnisky, Prediction and decoding of retinal ganglion cell responses with a probabilistic spiking model, Journal of Neuroscience 25, 11003 (2005), https://www.jneurosci.org/content/25/47/11003.full.pdf .
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8, 423 (1973).
  * Dominicis [1976] C. D. Dominicis, Techniques de renormalisation de la th´eorie des champs et dynamique des ph´enom‘enes critiques, Z Physik B 37, 247 (1976).
  * Janssen [1976] H. Janssen, On a lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties, Z Physik B 23, 377–380 (1976).
  * Chow and Buice [2015] C. Chow and M. Buice, Path integral methods for stochastic differential equations, J. Math. Neurosc. 5, 10.1186/s13408-015-0018-5 (2015).
  * Hertz  _et al._ [2016] J. A. Hertz, Y. Roudi, and P. Sollich, Path integral methods for the dynamics of stochastic and disordered systems, Journal of Physics A: Mathematical and Theoretical 50, 033001 (2016).
  * Helias and Dahmen [2020] M. Helias and D. Dahmen,  _Statistical field theory for neural networks_ , Vol. 970 (Springer, 2020).
  * Buice and Chow [2013a] M. A. Buice and C. C. Chow, Beyond mean field theory: statistical field theory for neural networks, Journal of Statistical Mechanics: Theory and Experiment 2013, P03003 (2013a).
  * Sompolinsky  _et al._ [1988] H. Sompolinsky, A. Crisanti, and H. J. Sommers, Chaos in random neural networks, Phys. Rev. Lett. 61, 259 (1988).
  * van Meegen  _et al._ [2021] A. van Meegen, T. Kühn, and M. Helias, Large-deviation approach to random recurrent neuronal networks: Parameter inference and fluctuation-induced transitions, Phys. Rev. Lett. 127, 158302 (2021).
  * Keup  _et al._ [2021] C. Keup, T. Kühn, D. Dahmen, and M. Helias, Transient chaotic dimensionality expansion by recurrent networks, Phys. Rev. X 11, 021064 (2021).
  * Buice and Chow [2013b] M. A. Buice and C. C. Chow, Dynamic finite size effects in spiking neural networks, PLOS Computational Biology 9, e1002872 (2013b).
  * Bressloff and Newby [2014] P. C. Bressloff and J. M. Newby, Path integrals and large deviations in stochastic hybrid systems, Phys. Rev. E 89, 042701 (2014).
  * Bressloff [2015] P. Bressloff, Path-integral methods for analyzing the effects of fluctuations in stochastic hybrid neural networks, J. Math. Neurosc. 5, 10.1186/s13408-014-0016-z (2015).
  * Buice  _et al._ [2010] M. A. Buice, J. D. Cowan, and C. C. Chow, Systematic Fluctuation Expansion for Neural Network Activity Equations, Neural Computation 22, 377 (2010), https://direct.mit.edu/neco/article-pdf/22/2/377/831286/neco.2009.02-09-960.pdf .
  * Stapmanns  _et al._ [2020] J. Stapmanns, T. Kühn, D. Dahmen, T. Luu, C. Honerkamp, and M. Helias, Self-consistent formulations for stochastic nonlinear neuronal dynamics, Phys. Rev. E 101, 042124 (2020).
  * Robert and Touboul [2016] P. Robert and J. Touboul, On the dynamics of random neuronal networks, Journal of Statistical Physics 165, 545 (2016).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum Field Theory and Critical Phenomena_ (Oxford University Press, 2002).
  * Brunel [2000] N. Brunel, Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons, Journal of computational neuroscience 8, 183 (2000).
  * Bressloff [2010b] P. C. Bressloff, Stochastic neural field theory and the system-size expansion, SIAM Journal on Applied Mathematics 70, 1488 (2010b), https://doi.org/10.1137/090756971 .
  * Schuecker  _et al._ [2018] J. Schuecker, S. Goedeke, and M. Helias, Optimal sequence memory in driven random networks, Phys. Rev. X 8, 041029 (2018).
  * Muscinelli  _et al._ [2019] S. P. Muscinelli, W. Gerstner, and T. Schwalger, How single neuron properties shape chaotic dynamics and signal transmission in random neural networks, PLOS Computational Biology 15, e1007122 (2019).
  * Ellis [2017] J. P. Ellis, Tikz-feynman: Feynman diagrams with tikz, Computer Physics Communications 210, 103 (2017).