Applications of information geometry to spiking neural network behavior
The space of possible behaviors complex biological systems may exhibit is
unimaginably vast, and these systems often appear to be stochastic, whether due
to variable noisy environmental inputs or intrinsically generated chaos. The
brain is a prominent example of a biological system with complex behaviors. The
number of possible patterns of spikes emitted by a local brain circuit is
combinatorially large, though the brain may not make use of all of them.
Understanding which of these possible patterns are actually used by the brain,
and how those sets of patterns change as properties of neural circuitry change
is a major goal in neuroscience. Recently, tools from information geometry have
been used to study embeddings of probabilistic models onto a hierarchy of model
manifolds that encode how model behaviors change as a function of their
parameters, giving a quantitative notion of "distances" between model
behaviors. We apply this method to a network model of excitatory and inhibitory
neural populations to understand how the competition between membrane and
synaptic response timescales shapes the network's information geometry. The
hyperbolic embedding allows us to identify the statistical parameters to which
the model behavior is most sensitive, and demonstrate how the ranking of these
coordinates changes with the balance of excitation and inhibition in the
network.
Jacob T. Crosser, Braden A. W. Brinkman
# Applications of information geometry to spiking neural network behavior

###### Abstract

The space of possible behaviors complex biological systems may exhibit is unimaginably vast, and these systems often appear to be stochastic, whether due to variable noisy environmental inputs or intrinsically generated chaos. The brain is a prominent example of a biological system with complex behaviors. The number of possible patterns of spikes emitted by a local brain circuit is combinatorially large, though the brain may not make use of all of them. Understanding which of these possible patterns are actually used by the brain, and how those sets of patterns change as properties of neural circuitry change is a major goal in neuroscience. Recently, tools from information geometry have been used to study embeddings of probabilistic models onto a hierarchy of model manifolds that encode how model behaviors change as a function of their parameters, giving a quantitative notion of "distances" between model behaviors. We apply this method to a network model of excitatory and inhibitory neural populations to understand how the competition between membrane and synaptic response timescales shapes the network's information geometry. The hyperbolic embedding allows us to identify the statistical parameters to which the model behavior is most sensitive, and demonstrate how the ranking of these coordinates changes with the balance of excitation and inhibition in the network.

## I Introduction

A major obstacle to understanding the computational underpinnings of the brain is the high dimensionality of its inputs--environmental stimuli such as light and sound--and its outputs--the activity of neurons and the organismal behaviors they enact [1]. The behavioral space of a neural circuit with \(N\) neurons is unmanageably large: the number of possible spike train patterns such a network can in principle produce over a trial of time length \(T\) divided into time bins of size \(\Delta t\) is of the order \(\sim\)\(2^{NT/\Delta t}\), assuming at most one spike per time bin. As \(\Delta t\to 0\), this output space becomes infinite-dimensional. However, the behavior of a neural population does not occupy this entire space, as activity is correlated across time and neurons, and the actual behavior of any given neural circuit constitutes just a subset of all possible observations. Perhaps surprisingly, analysis of experimental data has repeatedly found that under many conditions collective neural activity is low dimensional, often comprising less than \(\sim\)\(10^{2}\) dimensions of this infinite space [2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13].

Theoretical and computational work in neuroscience has largely focused on investigating the role that synaptic connections between neurons play in shaping the possible activity patterns of a network [14; 15; 16; 17; 18; 19; 20], which can be represented by manifolds (hyper-surfaces) in the behavioral output space of a neural circuit. These manifolds are complicated by the fact that many distinct neural circuits give rise to essentially identical patterns of activity [21; 22; 23; 24], meaning many different configurations are mapped to nearby points on these manifolds. Understanding how network activity and function changes as network properties or states change is a fundamental problem in neuroscience, and learning how to manipulate this activity most efficiently could lead to new and more effective treatments of neurological disorders.

Taming the possible behavioral repertoires of neural circuits by brute force simulations of network activity is computationally expensive and impractical for circuits larger than a few neurons. The tools of information geometry offer a possible means of representing network activity in an abstract way, but one that is easier to apply to larger networks and begin to understand how to most effectively move a network through its parameter space to achieve desired output behaviors [25; 26; 27; 28; 29; 30; 31; 32]. Note that this use of information geometry is a means of understanding the structure of complex models themselves, in contrast to applications of information theory in neuroscience as a modeling tool for understanding sensory coding [33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45].

Many models in complex biology generate a hierarchy of "hyperribbons" in their behavioral output space. These hyperribbons are manifolds with a few long directions of the manifold, representing "stiff" directions that separate disparate activity states, and many thin directions, which represent "sloppy" directions that describe networks with very similar behavior [46; 47; 48; 49; 50; 51; 52; 53]. These model manifolds come equipped with a natural metric that measures a sense of difference in behavior that is like a distance. We can determine the combinations of parameters that predict the bulk of the behavioral space of the network by identifying these model manifolds. This opens a path for better understanding of how to manipulate network properties to tune a circuit between different regimes of behavior.

In this work, we apply tools of information geometry to models of neural circuitry, and investigate how the bal
ance of single-neuron properties and the properties of the synaptic connections between neurons shape the hierarchy of possible behaviors of the networks. Specifically, we study how changing the membrane and synaptic time constants of the networks shape the manifold hierarchy. We also investigate how adjusting the balance of excitation and inhibition in the network change the rankings of the different hierarchical modes of the behavioral space. Previous work applying ideas from information geometry to neuroscience have primarily to study abstracted representations of spiking networks [25; 29], networks of rate models [30; 31; 28; 32], or neural field and pool models [26; 27]. By contrast, the work presented in this paper studies a class of leaky integrate-and-fire neurons--a commonly used modeling framework--with explicit consideration of some biophysical properties of individual neurons to make closer contact with the biological reality of neural systems. This is done by leveraging the specific properties of recently developed tools in information geometry [52; 53]. We organize the paper as follows: in Sec. II we introduce the class of stochastic spiking models we will be working with and the reduction to a population-based formalism. Then, in Sec. III, we give a self-contained explanation of the "isKL" embedding method introduced by [53], and how it applies to our population model. We detail the results of the application of the isKL method in Sec. IV, and finally discuss the interpretation and significance of our results and methodology in Sec. V.

## II Models

### Nonlinear Hawkes process

To model the spiking dynamics of individual neurons, we consider a nonlinear Hawkes process [54; 17]

\[\frac{dV_{i}}{dt}= -\tau_{m}^{-1}(V_{i}-\varepsilon_{i})+I_{i}\] \[+\tau_{s}^{-1}\left(\mu_{\text{ext}}-J_{\text{self}}\dot{n}_{i}( t)+\sum_{j=1}^{n}w_{ij}\dot{n}_{j}(t)\right) \tag{1a}\] \[\dot{n}_{i}(t)dt\sim\text{Poiss}[\phi(V_{i}(t))dt], \tag{1b}\]

where \(V_{i}\) is the membrane potential of neuron \(i\), \(\varepsilon_{i}\) is the leak reversal potential, \(w_{ij}\) is the strength of a synaptic connection from neuron \(j\) to neuron \(i\), and \(-J_{\text{self}}\) is an inhibitory self-coupling to implement post-spike refractory dynamics. The two currents \(\mu_{\text{ext}}\) and \(I_{i}\) represent an average current received from an external network and an experimentally injected current that differs by neuron, respectively. The process \(\dot{n}_{i}(t)\) is the spike train of neuron \(i\), and \(\phi(V_{i}(t))dt\) is the instantaneous firing rate nonlinearity that determines a Poisson event rate conditioned on the membrane potential of a given neuron. For the specific models studied here, \(\phi(x)=\frac{1}{2}(x+\sqrt{x^{2}+1/2})\). Finally, \(\tau_{m}\) and \(\tau_{s}\) are modulated parameters corresponding membrane and synaptic timescales, respectively. Eqn. 1a of this model assigns leaky integration dynamics to the membrane potential of each individual neuron, while Eqn. 1b assigns conditionally Poisson spiking dynamics to each neuron. Taken together, this model can be thought of as a soft-threshold leaky integrate-and-fire system.

Foreshadowing the coming analysis, we note that analytically calculating the statistical properties of the models in Eqn. 1 is generally intractable, and to make headway we will implement a Gaussian-process approximation of the network dynamics around the mean-field activity.

We can obtain a mean-field approximation of the steady-state solution for the membrane potential dynamics in Eqn. 1a by marginalizing out the spiking dynamics and assuming the distribution is sharply peaked around the most probable path of \(V_{i}(t)\). Assuming the network achieves a steady state at long times, this procedure gives us a set of transcendental equations that can be solved numerically:

\[V_{i}^{\text{mf}} =\varepsilon_{I}+\tau_{m}I_{i}\] \[+\frac{\tau_{m}}{\tau_{s}}\left(\mu_{\text{ext}}-J_{\text{self}} \phi(V_{i}^{\text{mf}})+\sum_{j}w_{ij}\phi(V_{j}^{\text{mf}})\right), \tag{2}\]

where the \(V_{i}^{mf}\), the solutions of these equations, are the mean-field predictions of the steady-state values of membrane potentials, with \(\phi(V_{i}^{mf})\) the corresponding mean-field prediction of the firing rates. We find the solutions to these transcendental equations using a forward-Euler integration scheme.

Following the prescription of Ref. [55; 56], the time-dependent distribution of model behaviors described in Eqn. 1 can be written in the form of a path integral,

\[P[\mathbf{V}(t),\mathbf{\dot{n}}(t)]=\int\mathfrak{D}[\tilde{\mathbf{V}}, \mathbf{\ddot{n}}]e^{-S[\tilde{\mathbf{V}},\mathbf{V},\tilde{\mathbf{n}}, \dot{\mathbf{n}}]}, \tag{3}\]

with an action \(S\) given by
s The Gaussian process approximation of the membrane dynamics in Eqn. 1a is obtained by marginalizing out the spiking dynamics from the action in Eqn. 4 and taking a saddle point approximation of the action around the mean-field solution in Eqn. 2 (see Appendices A & B for details). The resulting action corresponds to the Gaussian stochastic process given by [17]

\[d\mathbf{V}=\mathbf{A}\left(\mathbf{V}^{mf}-\mathbf{V}\right)dt+\mathbf{ \Sigma}d\mathbf{W}_{t} \tag{5a}\] \[A_{ij}=\delta_{ij}\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self} }\phi^{\prime}(V_{j}^{\text{mf}})\right)-\tau_{s}^{-1}w_{ij}\phi^{\prime}(V_{ j}^{\text{mf}})\] (5b) \[\left(\mathbf{\Sigma}\mathbf{\Sigma}^{T}\right)_{ij}=\tau_{s}^{- 2}\sum_{k}\Biggl{[} \left(-\delta_{ik}J_{\text{self}}+w_{ik}\right)\] \[\times\left(-\delta_{jk}J_{\text{self}}+w_{jk}\right)\phi(V_{k}^ {\text{mf}})\Biggr{]} \tag{5c}\]

where \(d\mathbf{W}_{t}\) is a standard Wiener process and we use the Ito convention. We note that Eqn. 5a is an Ornstein-Uhlenbeck (OU) process, albeit one in which the drift and diffusion matrices are dependent on the mean-field values of the membrane potential.

In principal, the network modeled in Eqn. 1 and approximated in Eqn. 5 could be of arbitrary size. To make our information geometric analysis tractable, however, we will reduce the model to a three-population model, comprising excitatory and inhibitory populations, and a single neuron targeted with an injected current; this is depicted diagrammatically in Fig. 1. We start by considering the connectivity matrix to be random with each entry being a Bernoulli variable with probability \(p\) being scaled by a connection type-dependent value \(w_{IJ}\). To produce the more tractable reduced model, we take a population-averaging approach to the approximated process in Eqn. 5. We now use an uppercase subscript to denote a population averaged variable. For example

\[V_{I}\equiv\frac{1}{N_{I}}\sum_{i\in I}V_{i}(t),\]

where we will use uppercase indices \(I,J,K\in\{0,1,2\}\) to denote the different populations, with \(I=0\) the single test neuron, \(I=1\) the excitatory population, and \(I=2\) the inhibitory population. The dynamics of the population-averaged membrane potentials under the Gaussian approximation now follow a lower-dimensional version of Eqn. 5a with drift and diffusion matrices given by

\[A_{IJ}= \delta_{IJ}(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{\prime }(V_{I}^{\text{mf}}))\] \[-\tau_{s}^{-1}pw_{IJ}N_{J}\phi^{\prime}(V_{J}^{\text{mf}}) \tag{6a}\]

Figure 1: **Network model** (A) A graphical representation of the network architecture being studied. (B) An example raster plot generated from an extended network of spiking neurons modeled by Eqn.1
\[\left(\mathbf{\Sigma\Sigma}^{T}\right)_{IJ}= \tau_{s}^{-2}\sum_{K=0,1,2}\Bigg{[}\left(-\delta_{IK}\frac{J_{\rm self }}{N_{K}}+pw_{IK}\right)\] \[\times\left(-\delta_{JK}\frac{J_{\rm self}}{N_{K}}+pw_{JK}\right)N _{K}\phi(V_{K}^{\rm mf})\Bigg{]} \tag{6b}\] \[\approx\tau_{s}^{-2}\sum_{K=1,2}p^{2}w_{IK}w_{JK}N_{K}\phi(V_{K}^ {\rm mf})\]

The approximation in the last line above comes from the fact that \(N_{1},\ N_{2}\gg 0\). The population-averaged mean-field equations are now

\[V_{I}^{\rm mf}= \varepsilon_{I}+\tau_{m}I_{I}+\frac{\tau_{m}}{\tau_{s}}\mu_{\rm ext}\] \[+\frac{\tau_{m}}{\tau_{s}}\left(-J_{\rm self}\phi(V_{I}^{\rm mf}) +\sum_{J=0,1,2}pw_{IJ}N_{J}\phi(V_{J}^{\rm mf})\right) \tag{6c}\]

We formally derive the Gaussian-process approximation of the full-network (Eqn. 5) and the population-averaged approximation (Eqn. 6) in Appendix A. We also note that the statistics of the model in Eqn. 6 are equivalent to those derived by first taking a population average of the membrane potential dynamics and then applying the Gaussian approximation framework. This second derivation is provided in Appendix B.

The spiking model we study centers around a balanced network, specifically a network that is not finely tuned. This notion of fine-tuning arises from a standard derivation of balance equations for the model system (see Appendix D). In short, we can look at the average external input \(\kappa_{I}\) into population \(I\). For our model, we can approximate \(\kappa_{I}\) to leading order as:

\[\tau_{s}^{-1}\kappa_{I}\approx\sqrt{N} \left(\frac{1}{\sqrt{N}}\left(I_{I}+\tau_{s}^{-1}\mu_{\rm ext} \right)\right.\] \[\left.+\tau_{s}^{-1}\left\{pw_{I1}\frac{N_{1}}{\sqrt{N}}\phi(V_{ 1})+pw_{I2}\frac{N_{2}}{\sqrt{N}}\phi(V_{2})\right\}\right)\]

where \(V_{I}\) is the population-averaged membrane potential for population \(I\). For the model to be in a balanced state, the variance of the synaptic input should be \(\mathcal{O}(N^{0})\) which in turn implies the synaptic weights should scale as \(w_{IJ}\sim 1/\sqrt{N}\). Additionally, we assume that \(N_{I}\propto N\) and \(\langle I_{I}\rangle,\ \mu_{\rm ext}\propto\sqrt{N}\). The balanced state of the model also requires that all \(\kappa_{I}\) be \(\mathcal{O}(1)\). For this to be true as \(N\to\infty\), the terms in the parentheses must vanish. This provides gives us a linear system that uniquely defines (\(\phi(V_{1})\), \(\phi(V_{2})\)):

\[-\begin{bmatrix}I_{1}+\tau_{s}^{-1}\mu_{\rm ext}\\ I_{2}+\tau_{s}^{-1}\mu_{\rm ext}\end{bmatrix}=\frac{1}{\tau_{s}}\begin{bmatrix} pw_{11}N_{1}&pw_{12}N_{2}\\ pw_{21}N_{1}&pw_{22}N_{2}\end{bmatrix}\begin{bmatrix}\phi(V_{1})\\ \phi(V_{2})\end{bmatrix} \tag{7}\]

From this set of equations, we derive two cases. First, if the matrix on the right-hand side of Eqn. 7 is singular and neither of the columns of the matrix are trivially the zero-vector, the columns must be scalar multiples of each other. We refer to this as a "fine-tuned" spiking model. If the left-hand side of Eqn. 7 is also a multiple of the columns, the system admits an infinite set of solutions (\(\phi(V_{1})\), \(\phi(V_{2})\)). Otherwise, it admits no solution. Such a network is thus finely-tuned to specific inputs. In contrast to this, we have "un-tuned" spiking models. In this case, the matrix on the right-hand side of Eqn. 7 is invertible and the system admits a unique solution (\(\phi(V_{1})\), \(\phi(V_{2})\)). This in effect applies constraints on the values of \(\{w_{IJ}\}\), which we refer to as the balance equations Eqs. D2 & D3 (see Appendix D for a derivation). Moving forward, we consider only spiking models derived from a balanced, un-tuned network. We also introduce a linear non-spiking model that will serve as a baseline comparison.

### Linear non-spiking model

Although the Gaussian process approximation of the spiking network will have a Gaussian steady-state distribution of the membrane potentials, the parameters of this distribution vary nonlinearly with the self-consistent mean-field solutions. To demonstrate that the behaviors we observe are consequences of the mean-field treatment of the spiking network, and not just the behavior of Gaussian processes, we also construct a simpler model of networked, linear non-spiking (or "graded potential") neurons. We assume the neurons are injected with large numbers of synaptic input that sum together to be approximately Gaussian, with non-zero mean \(\mu_{\rm ext}\), creating a stochastic system with dynamics described by:

\[\frac{dV_{i}}{dt}= -\tau_{m}^{-1}(V_{i}-\varepsilon_{I})+I_{i}+\tau_{s}^{-1}\mu_{\rm ext }-\tau_{s}^{-1}J_{\rm self}\phi(V_{i})\] \[\qquad\qquad+\tau_{s}^{-1}\sum_{j}w_{ij}\phi(V_{j})+\xi_{i}(t) \tag{8}\]

Here, the transfer function \(\phi(\cdot)\) is simply the identity function (i.e. \(\phi(x)=x\)). The processes \(\xi_{i}(t)\) are zero-mean Gaussian noise synaptic input from neurons external to the network being examined, and thus they scale with \(\tau_{s}^{-1}\). We define the covariance of the noise processes \(\{\xi_{i}(t)\}\) as follows.

\[\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\tau_{s}^{-2}\delta_{ij}\mu_{\rm ext }\delta(t-t^{\prime})\]

After population-averaging, the non-spiking model becomes another OU process:

\[d\mathbf{V}= \mathbf{A}\left(\mathbf{A}^{-1}\left(\tau_{s}^{-1}\mu_{\rm ext}+ \tau_{m}^{-1}\varepsilon_{I}+\mathbf{I}\right)-\mathbf{V}\right)dt+\mathbf{\Sigma}d \mathbf{W}_{t}\] \[= \mathbf{A}\left(\mu-\mathbf{V}\right)dt+\mathbf{\Sigma}d\mathbf{W}_{t}. \tag{9}\]

The drift and diffusion matrices are defined as follows

\[A_{IJ}= \delta_{IJ}\tau_{m}^{-1}+\tau_{s}^{-1}w_{IJ}^{*}\] \[w_{IJ}^{*} =-\delta_{IJ}J_{\rm self}+pw_{IJ}^{\rm mod}N_{J}\] \[\left(\Sigma\Sigma^{T}\right)_{IJ} =\tau_{s}^{-2}\delta_{IJ}\frac{\mu_{\rm ext}}{N_{I}}\]

Here and in the following sections, \(\mathbf{w}^{*}\) denotes the effective connectivity matrix for the linear non-spiking models. The values of \(\mathbf{w}^{\rm mod}\) are modulated depending on
the desired excitation-inhibition conditions, which will be discussed in Sec. II.4.

The linear form of the population-averaged non-spiking model permits more analytic study than the corresponding spiking models. Ornstein-Uhlenbeck processes like those in Eqns. 6 and 9 admit a Gaussian steady-state distribution if all eigenvalues of the drift matrix are positive [57]. From the form of the drift matrix for the linear model (Eqn. 9) there is a correspondence between eigenvalues of the drift matrix \(\mathbf{A}\) and the connectivity matrix \(\mathbf{w}^{*}\).

\[\lambda_{i,\mathbf{A}}=\tau_{m}^{-1}-\tau_{s}^{-1}\lambda_{i,\mathbf{w}^{*}}.\]

From the stationarity condition on the eigenvalues of \(\mathbf{A}\) and this correspondence between eigenvalues of \(\mathbf{A}\) and \(\mathbf{w}^{*}\), we can derive a stability boundary for the \((\tau_{m}^{-1},\tau_{s}^{-1})\) inverse timescale-space

\[\tau_{m}^{-1}>\tau_{s}^{-1}\lambda_{\mathbf{w}^{*}},\ \forall\ \lambda_{ \mathbf{w}^{*}}. \tag{10}\]

The loss of stability observed in OU processes often corresponds to a non-stationary regime in which the random variables may grow without bound. As the firing rate nonlinearity \(\phi(x)\) used in the spiking model is quasi-linear in the \(x>0\) regime, we expect the stability of the spiking models to be similar to the non-spiking models when they have the same E/I-dependent connectivity given by \(w_{IJ}\).

### Stationary distributions

As mentioned above, the stationary distributions admitted by Ornstein-Uhlenbeck processes are Gaussian when they exist [57]. Consider a general \(N\)-dimensional OU process:

\[d\mathbf{X}=\mathbf{A}\left(\mu-\mathbf{X}\right)dt+\mathbf{\Sigma}d\mathbf{W }_{t}\]

The stationary distribution, when it exists, is described by the multivariate normal probability density [57]

\[p\left(\mathbf{X}\right)=\frac{1}{(2\pi)^{N/2}\sqrt{\det(\mathbf{C})}}e^{- \frac{1}{2}(\mathbf{X}-\mu)^{T}\mathbf{C}^{-1}(\mathbf{X}-\mu)},\]

where the stationary covariance \(\mathbf{C}\) is given by the solution to the matrix equation [57]

\[\mathbf{\Sigma}\mathbf{\Sigma}^{T}=\mathbf{A}\mathbf{C}+\mathbf{C}\mathbf{A}^ {T}.\]

In practice, the stationary covariance matrix can by found by linearizing the matrix equation and solving the resulting linear system numerically.

### Network architectures

Now, we turn back to our network models. We consider a population of excitatory and inhibitory neurons in which a single excitatory target neuron is injected with an external driving current. The full network contains \(N=1000\) sparsely connected neurons. We condensed the full network model into representative 3-neuron network by population-averaging, as depicted in Fig.1A and described in Eqn. 6, representing the excitatory target neuron, the excitatory population, and the inhibitory population. Table 1 contains descriptions and numerical values for the parameters used in the present study.

In addition, we would like to adjust the relative recurrent excitation and inhibition in the networks. To accomplish this, the base connection weights given in Table 1 are scaled by a ratio \(r>0\) depending on the desired activity regime:

\[\mathbf{w}_{Xe}^{\mathrm{mod}}=r_{e}(r)w_{Xe,\mathrm{base}}=\begin{cases}rw_ {Xe,\mathrm{base}}&\text{if }r\geq 1\\ w_{Xe,\mathrm{base}}&\text{otherwise}\end{cases} \tag{11a}\] \[\mathbf{w}_{Xi}^{\mathrm{mod}}=r_{i}(r)w_{Xi,\mathrm{base}}= \begin{cases}w_{Xi,\mathrm{base}}&\text{if }r\geq 1\\ \frac{1}{r}w_{Xi,\mathrm{base}}&\text{otherwise}\end{cases}. \tag{11b}\]

The ratio serves to boost the recurrent excitatory weights in the excitatory regime (\(r>1\)) and the recurrent inhibitory weights in the inhibitory regime (\(r<1\)), through the functions \(r_{e}(r)\) and \(r_{i}(r)\), respectively.

The connection matrices \(\mathbf{w}\) and \(\mathbf{w}^{*}\) of the population-averaged spiking and non-spiking models, respectively, are now constructed from the full-network parameters and scaling of excitation and inhibition. All matrices \(\mathbf{w}\) and \(\mathbf{w}^{*}\) use the same indexing with \(I=0\) denoting the target neuron "population," \(I=1\) denoting the remaining excitatory neurons, and \(I=2\) denoting all inhibitory neurons. The connection strengths used in the linear non-spiking model in Eqn. 9 is then given by
\[\mathbf{w}^{*}= -\begin{bmatrix}J_{\mathrm{self}}&0&0\\ 0&J_{\mathrm{self}}&0\\ 0&0&J_{\mathrm{self}}\end{bmatrix}+\frac{1}{\sqrt{pN}}\begin{bmatrix}pw_{ee}&p \left(N_{e}-1\right)w_{ee}&pN_{i}w_{ei}\\ pw_{ee}&p\left(N_{e}-1\right)w_{ee}&pN_{i}w_{ei}\\ pw_{ie}&p\left(N_{e}-1\right)w_{ie}&pN_{i}w_{ii}\end{bmatrix}\] \[=-\begin{bmatrix}J_{\mathrm{self}}&0&0\\ 0&J_{\mathrm{self}}&0\\ 0&0&J_{\mathrm{self}}\end{bmatrix}+\frac{1}{\sqrt{pN}}\begin{bmatrix}pr_{e}(r )w_{ee,\mathrm{base}}&pr_{e}(r)\left(N_{e}-1\right)w_{ee,\mathrm{base}}&pr_{i} (r)N_{i}w_{ei,\mathrm{base}}\\ pr_{e}(r)w_{ee,\mathrm{base}}&pr_{e}(r)\left(N_{e}-1\right)w_{ee,\mathrm{base }}&pr_{i}(r)N_{i}w_{ei,\mathrm{base}}\\ pr_{e}(r)w_{ie,\mathrm{base}}&pr_{e}(r)\left(N_{e}-1\right)w_{ie,\mathrm{ base}}&pr_{i}(r)N_{i}w_{ii,\mathrm{base}}\end{bmatrix}. \tag{12}\]

The \(1/\sqrt{pN}\) scaling of the connection weights arises from the balance conditions mentioned at the end of Sec. II.1 and derived in Appendix D. The connection matrices used by the spiking models described generally in Eqn. 6 is given by

\[\mathbf{w} =\frac{1}{\sqrt{pN}}\mathbf{w}^{\mathrm{mod}}\] \[=\frac{1}{\sqrt{pN}}\begin{bmatrix}r_{e}(r)w_{ee,\mathrm{base}}&r _{e}(r)w_{ee,\mathrm{base}}&r_{i}(r)w_{ei,\mathrm{base}}\\ r_{e}(r)w_{ee,\mathrm{base}}&r_{e}(r)w_{ee,\mathrm{base}}&r_{i}(r)w_{ei, \mathrm{base}}\\ r_{e}(r)w_{ie,\mathrm{base}}&r_{e}(r)w_{ie,\mathrm{base}}&r_{i}(r)w_{ii, \mathrm{base}}\end{bmatrix}. \tag{13}\]

Finally, we would like a measure of the balance of excitation and inhibition ("E/I") within a class of models. As each model type corresponds to many particular models with different values of the inverse timescales (\(\tau_{m}^{-1},\tau_{s}^{-1}\)), we require a proxy measure for the E/I ratio to describe the whole class. In line with the method for adjusting the relative strength of recurrent excitation and inhibition introduced above, we assign a ratio of connection weights into the bulk excitatory population for a given model and a given modulation \(r\). For the non-spiking models, we give the log-ratio \(R\) of these weights

\[R=\log_{10}\left|\frac{\mathbf{w}_{2,1}^{*}+\mathbf{w}_{2,2}^{*}}{\mathbf{w}_{ 2,3}^{*}}\right|\]

To make an accurate comparison to the non-spiking models, the E/I values for the spiking models are reported using this same measure (for a given value of modulation parameter \(r\)).

We note here that the same balanced-network calculations (see Appendix D) that gave rise to the definitions of "fine-tuned" and "un-tuned" formally define a notion of balance. A balanced spiking network based on the model architecture used here must satisfy constraints on the weights of \(\mathbf{w}\) (Eqn. 13), either Eqns. D2 or Eqns. D3. The base connection weights for the unadjusted network--i.e. \(r=1\) in Eqns. 11--were chosen to meet these balance criteria, and the functions \(r_{e}(r)\) and \(r_{i}(r)\) serve to tilt the excitation-inhibition balance with respect to this measure.

\begin{table}
\begin{tabular}{|l||l|l|} \hline Parameter & Description & Value \\ \hline \hline \(N\) & Total number of neurons & \(1000\) \\ \hline \(N_{e}\) & Number of excitatory neurons & \(0.8N\) \\ \hline \(N_{i}\) & Number of inhibitory neurons & \(0.2N\) \\ \hline \(p\) & Probability of a directional synaptic connection \(w_{ij}\) between any two neurons & \(0.1\) \\ \hline \(-J_{\mathrm{self}}\) & The self connection for a neuron of type \(E/I\) designed to capture post-spike re- & \(-5\) \\  & factory dynamics & \\ \hline \(\varepsilon_{I}\) & Leak reversal potential for neuron \(i\) & \(0\) \\ \hline \(I_{I}\) & Injected current impinging on neuron population \(I\) & \(\begin{cases}0.02\text{ if target}\\ 0\text{ otherwise}\end{cases}\) \\ \hline \(\mu_{\mathrm{ext}}\) & Mean input from network-external neurons & \(0.1\) \\ \hline \(w_{ee,\mathrm{base}}\) & The total expected synaptic input weight from exc. neurons to exc. neurons & \(285\) \\ \hline \(w_{ie,\mathrm{base}}\) & The total expected synaptic input weight from exc. neurons onto inh. neurons & \(300\) \\ \hline \(w_{ei,\mathrm{base}}\) & The total expected synaptic input weight from inh. neurons to exc. neurons & \(-902.5\) \\ \hline \(w_{ii,\mathrm{base}}\) & The total expected synaptic input weight from inh. neurons to inh. neurons & \(-950\) \\ \hline \(\phi(x)\) & Firing rate transfer function & \(\begin{cases}x\text{ if Non}-\text{spiking}\\ \frac{x+\sqrt{x^{2}+\frac{1}{2}}}{2}\text{ if Spiking}\end{cases}\) \\ \hline \(\tau_{m}\) & Membrane timescale & variable \\ \hline \(\tau_{s}\) & Synaptic timescale & variable \\ \hline \end{tabular}
\end{table}
Table 1: **Model parameters** Descriptions and numerical values for the parameters for the non-spiking and spiking model types.
### Timescale sampling

To embed and visualize the model manifolds of interest, we must sample points on the manifold characterized by different values of the two modulated parameters. We do this by sampling a portion of the inverse-timescale parameter space that satisfies the stability condition given by Eq. (10) and where both inverse-timescales are positive. We apply a curvilinear grid to this region, uniformly sampling the radial and angular components. The radial distance components \(d\) of the grid are taken over a fixed range:

\[d\in[0.0025,\,0.03]\text{ ms}^{-1}\]

To apply both the stability boundary and positivity constraints, the lower bound of the angular component \(\alpha\) of the sample grid is set to a fixed value while the upper bound is set either by the stability boundary described by Eqn. 10 or to a fixed value, whichever is more stringent:

\[\tan(\alpha)\in\left[0.1,\,\min\left(\frac{1}{\max\{\lambda_{\mathbf{w}^{*}} \}},500\right)\right]\]

The conditions for this maximal sampling are summarized in Table 2. It is important to note that the stability boundary is determined by the eigenvalues of the connectivity matrix \(\mathbf{w}^{*}\), and thus the stability boundary and the sampling region are affected by the induced E/I balance is adjusted through its affect on \(\mathbf{w}^{*}\). The spiking models use the connection matrix from the equivalent non-spiking model to set the sampling range. The maximal sampling scheme is depicted diagrammatically in Fig. 2.

After the maximal sampling of parameter space for each model type for each E/I condition, sample points from the inverse-timescale space are subject to further exclusionary criteria. For both the spiking- and non-spiking-type models, sample points are excluded if they cause either the drift matrix \(\mathbf{A}\) or the covariance matrix \(\mathbf{C}\) to become singular. The singularities in these matrices have been observed to occur numerically near the theoretical stability boundary (Eqn. 10). In addition, sample points for the spiking models are excluded if the Euler integration used to find the mean-field solutions to Eqn. 6c does not converge. The integration is determined to be numerically non-convergent if the rate of change of the system either exceeds a predetermined value during integration or does not hit a convergence threshold before reaching the maximum number of steps.

All model manifolds studied here were generated from between 210,000 and 211,000 sampled parameter pairs.

## III Iskl embedding

In this section, we recapitulate the methods developed by Teoh and colleagues [53]. This framework revolves around using the symmetric Kullback-Liebler divergence \(D_{sKL}\) as a measure of separation for different probabilistic models of the same form but with different parameters:

\[D_{sKL}(\theta,\theta^{\prime})= D_{KL}(\theta:\theta^{\prime})+D_{KL}(\theta^{\prime}:\theta)\] \[=\mathbb{E}_{\theta}\left[\ln\frac{p(x|\theta)}{p(x|\theta^{ \prime})}\right]-\mathbb{E}_{\theta^{\prime}}\left[\ln\frac{p(x|\theta)}{p(x| \theta^{\prime})}\right].\]

Teoh _et al._ apply this measure to exponential family models, which have the general form

\[p(x|\theta)=\exp\left[\sum_{i=1}^{n}t_{i}(x)\eta_{i}\left(\theta\right)+k(x)- A\left(\eta\left(\theta\right)\right)\right],\]

where \(\{\eta_{i}\left(\theta\right)\}\) are the \(n\) natural parameters of the model and \(\{t_{i}(x)\}\) are the corresponding sufficient statistics.

Figure 2: **Maximal sampling of the inverse timescale plane:** We sample pairs of inverse-timescale values from the depicted region of the \(\tau_{m}^{-1}\)-\(\tau_{s}^{-1}\) plane. Sampled points are distributed evenly on a curvilinear grid between a predefined lower bound and the stability boundary for the specific connection matrix being used (A). If sampling to the stability boundary would produce samples with negative timescales, the space is instead sampled up to a predefined value in the angular direction (B).
The \(D_{sKL}\) for exponential family models can be analytically decomposed into a _finite_ number of component functions

\[D_{sKL}[\theta,\theta^{\prime}]=\sum_{i=1}^{n}\left\{\left[\mathcal{T}_{i}^{+}( \theta)-\mathcal{T}_{i}^{+}(\theta^{\prime})\right]^{2}-\left[\mathcal{T}_{i}^ {-}(\theta)-\mathcal{T}_{i}^{-}(\theta^{\prime})\right]^{2}\right\}.\]

These component functions form a set of \(n\) space-like (\(\mathcal{T}_{i}^{+}\)) and \(n\) time-like (\(\mathcal{T}_{i}^{-}\)) coordinates by which the model manifold may be embedded in a Minkowski-like behavioral space [53]. These coordinate functions are given in terms of the natural parameters and sufficient statistics [53] by

\[\mathcal{T}_{i}^{\pm}=\frac{1}{2}\left[\eta_{i}\left(\theta\right)\pm\langle t _{i}(x)\rangle_{\theta}\right].\]

Alternatively, we may use an isometric embedding given by shifting and rotating the manifold [53]:

\[T_{i}^{\pm}(\theta)=\frac{1}{2}\left\{\lambda_{i}\left[\eta_{i}(\theta)- \overline{\eta_{i}}\right]\pm\frac{1}{\lambda_{i}}\left[\langle t_{i}\rangle_ {\theta}-\overline{\langle t_{i}\rangle}\right]\right\}. \tag{14}\]

We use \(T^{\pm}\) to distinguish the isometric embedding coordinates from the unscaled coordinates \(\mathcal{T}^{\pm}\). Here, an over-bar denotes a mean over sampled parameters and \(\lambda_{i}=\left[\operatorname{var}\left(\langle t_{i}\rangle\right)/ \operatorname{var}\left(\eta_{i}\right)\right]^{1/4}\). These coordinates can be understood as an alternative definition of the exponential family. We can straightforwardly express the log-likelihood function for an exponential family in terms of the isKL coordinates:

\[\ln p(x|\theta)= \ln k(\mathbf{x})+\sum_{i}\mathcal{T}_{i}^{+}(\theta)t_{i}( \mathbf{x})\] \[\qquad+\sum_{i}\mathcal{T}_{i}^{-}(\theta)t_{i}(\mathbf{x})-A(\theta)\] \[=g(\mathbf{x})+\sum_{i}T_{i}^{+}(\theta)t_{i}(\mathbf{x})\] \[\qquad+\sum_{i}T_{i}^{-}(\theta)t_{i}(\mathbf{x})-A(\theta),\]

where \(g(\mathbf{x})=\ln k(\mathbf{x})+\sum_{i}\overline{\eta}_{i}t_{i}(\mathbf{x})\). The authors [53] show that the coordinates \(\{T_{i}\}\) can also be understood in relation to the data visualization procedure multidimensional scaling (MDS). In standard MDS the data points are recorded data, whereas here each "data point" corresponds to the full distribution of an exponential family evaluated at a specific set of parameters. The double mean centered matrix of MDS can be constructed in this context from the pairwise separation matrix measured by the symmetric KL-divergence \(\mathbf{D}_{c}=-\mathbf{P}\mathbf{D}_{sKL}\mathbf{P}\) with \(\mathbf{P}_{ij}=1/n-\delta_{ij}\). In the continuous sampling limit, the eigenvalue problem for MDS is formulated as an integral equation:

\[\int D_{c}(\tilde{\theta},\theta)v(\theta)d\rho(\theta)=\Lambda v(\tilde{ \theta}), \tag{15}\]

where \(d\rho(\theta)=\rho(\theta)d\theta\) is the measure of the distribution of parameters \(\theta\). Teoh and colleagues show [53] that the coordinates \(T_{i}^{\pm}\) are solutions of this eigenvalue problem with corresponding eigenvalues

\[\Lambda_{i}^{\pm}=\frac{1}{2}\left[\operatorname{Cov}(\eta_{i},\langle t_{i} \rangle)\pm\sqrt{\operatorname{var}(\eta_{i})\operatorname{var}(\langle t_{i }\rangle)}\right] \tag{16}\]

This procedure produces an embedding with only a finite and relatively small number of non-zero modes, contrasting sharply with the infinite or data-proportional embedding produced by other methods for continuous or discrete parameter sampling, respectively [53].

We complement this perspective by viewing this embedding procedure as an eigenmode expansion of the conditional probability \(p(x|\theta)\) around the marginalized distribution \(p(x)\) for a given prior on the parameters \(\theta\):

\[p(x|\theta)=p(x)+\sum_{i}c_{i}^{+}(x)T_{i}^{+}(\theta)+\sum_{i}c_{i}^{-}(x)T_{ i}^{-}(\theta) \tag{17}\]

By defining the inner product of functions on \(\Theta\) as

\[\langle f(\theta),g(\theta)\rangle=\int d\rho(\theta)f(\theta)g(\theta),\]

the modes \(\sqrt{\mu(\theta)}v(\theta)\) of Eq. (15) can be shown to be orthogonal as long as the corresponding eigenvalues are distinct. Thus, the coordinate functions \(T_{i}^{\pm}(\theta)\) are orthogonal with respect to the weight \(\rho(\theta)\). Taking advantage of this orthogonality of the coordinate functions, it follows

\[\int d\rho(\theta)p(x|\theta)T_{j}^{\pm}(\theta)\] \[=p(x)\int d\rho(\theta)T_{j}^{\pm}(\theta)+\sum_{i,\pm}\int d\rho (\theta)c_{j}^{\pm}(x)T_{j}^{\pm}(\theta)T_{i}^{\pm}(\theta)\] \[=c_{j}^{\pm}(x)\int d\rho(\theta)\left(T_{j}^{\pm}(\theta) \right)^{2}\]

The first term on the right-hand side vanishes because the mean of each coordinate function is zero by construction, while only the \(i=j\) term from the sum survives due to the orthogonality. Thus, we may calculate the coefficient functions \(c_{i}^{\pm}(x)\) as

\[c_{i}^{\pm}(x)=\langle T_{i}^{\pm},T_{i}^{\pm}\rangle^{-1}\int d\rho(\theta)p (x|\theta)T_{i}^{\pm}(\theta).\]

\begin{table}
\begin{tabular}{|l||l|l||} \hline  & Radial distance \(d\) & Angle \(\alpha\) \\ \hline \hline Minimum & & \\ Value & 0.0025 ms\({}^{-1}\) & 0.1 \\ \hline Maximum & & \\ Value & 0.03 ms\({}^{-1}\) & \(\min\left(\frac{1}{\max\{\lambda_{\mathbf{y}^{-1}}\}},500\right)\) \\ \hline Number of & & \\ Samples & 301 & 701 \\ \hline \end{tabular}
\end{table}
Table 2: **Maximal Sampling Parameters** Descriptions and numerical values for the parameters that are constant across the non-spiking, spiking with fine-tuning, and spiking without fine-tuning model types.
In this work we focus on applying these embedding methods to the stationary distributions of the various network models, both multivariate normal within our mean-field approximation. For a \(M\)-dimensional multivariate normal distribution with a set of means \(\{\mu_{i}\}\) and covariance values \(\{C_{ij}\}\), the \(M(M+3)/2\) distinct natural parameters and sufficient statistics are given by

\[\eta=\begin{bmatrix}\sum_{i}C_{1i}^{-1}\mu_{i}\\ \vdots\\ \sum_{i}C_{1i}^{-1}\mu_{i}\\ -\frac{1}{2}C_{11}^{-1}\\ \vdots\\ -\frac{1}{2}C_{M1}^{-1}\\ -\frac{1}{2}C_{22}^{-1}\\ \vdots\\ -\frac{1}{2}C_{MM}^{-1}\end{bmatrix},\quad\langle t_{i}\rangle_{\theta}= \begin{bmatrix}\langle x_{1}\rangle\\ \vdots\\ \langle x_{M}\rangle\\ \langle x_{1}^{2}\rangle\\ \vdots\\ \langle x_{M}x_{1}\rangle\\ \vdots\\ \langle x_{2}^{2}\rangle\end{bmatrix} \tag{18}\]

Before we present the embedding and analysis of the

Figure 3: **Visualizations of the isKL embedding coordinate functions:** (A) The two coordinate functions for the 1-dimensional exponential model. (B-E) Coordinate functions for the 1-dimensional Gaussian model. (B) \(T_{1}^{+}\). (C) \(T_{1}^{-}\). (D) \(T_{2}^{+}\). (E) \(T_{2}^{-}\). (F) The model manifold for the 1-dimensional exponential model colored by \(\log\upsilon\). (G) The model manifold for the 1-dimensional Gaussian model projected onto just \(T_{1}^{\pm}\), colored by \(\mu\).
models from Section II, we provide two simpler models as illustrative examples that are related to the Poissonian and Gaussian characteristics of our model.

### Example: 1-dimensional exponential model

Let \(X\) be exponentially distributed with rate \(\upsilon\), i.e. \(X\sim\text{Exp}(\upsilon)\). In the exponential family formalism, we have

\[\eta=-\upsilon,\ \ \langle t_{i}\rangle_{\theta}=\langle x\rangle_{\upsilon}= \upsilon^{-1},\ \ k(x)=1,\ \ A(\upsilon)=-\ln\upsilon.\]

The isKL embedding coordinates for this model are one-dimensional functions given by

\[T^{\pm}(\upsilon)=\frac{1}{2}\left\{\lambda\left[\overline{\upsilon}-\upsilon \right]\pm\frac{1}{\lambda}\left[\upsilon^{-1}-\overline{\upsilon^{-1}} \right]\right\}.\]

These embedding functions are shown in Fig. 3A using a parameter distribution \(\rho(\upsilon)=\left(8\upsilon\ln 10\right)^{-1}\) with support \(\upsilon\in[10^{-5},10^{5}]\) for illustration.

We may also explicitly calculate the coefficients \(c^{\pm}(x)\) for this example,

\[c^{\pm}(x)= \langle T^{\pm},T^{\pm}\rangle^{-1}\Bigg{[}\left(\frac{\lambda}{ 2}\overline{\upsilon}\mp\frac{1}{2\lambda}\overline{\upsilon^{-1}}\right) \left(-\frac{dZ}{dx}\right)\] \[-\frac{\lambda}{2}\frac{d^{2}Z}{dx^{2}}\pm\frac{1}{2\lambda}Z \Bigg{]},\]

where \(Z(x)\equiv\int e^{-vx}d\rho(\upsilon)\) is the moment-generating function of the distribution \(\rho(\upsilon)\) with source \(-x\). The full model manifold is depicted in Fig. 3F, where points are colored by the logarithm of the rate parameter \(\upsilon\). Here, we see the manifold is neatly broken into two branches corresponding to a low event-rate (\(\log\upsilon<0\)) and a high event-rate (\(\log\upsilon>0\)).

### Example: 1-dimensional Gaussian model

Let \(X\) be normally distributed as \(X\sim\mathcal{N}(\mu,\sigma)\). In the exponential family formalism, we have

\[\eta=\begin{bmatrix}\mu/\sigma^{2}\\ -\sigma^{-2}\end{bmatrix},\ \ \langle t_{i}\rangle_{\theta}=\begin{bmatrix} \langle x\rangle\\ \langle x^{2}\rangle\end{bmatrix}=\begin{bmatrix}\mu\\ \sigma^{2}+\mu^{2}\end{bmatrix}\] \[k(x)=\frac{1}{\sqrt{2\pi}},\ \ \ A(\mu,\sigma)=\frac{\mu^{2}}{2 \sigma^{2}}+\ln\sigma.\]

The isKL embedding coordinates are then two-dimensional functions

\[T^{\pm}_{1}(\mu,\sigma)=\frac{1}{2}\left\{\lambda_{1}\left[\frac{\mu}{\sigma^{ 2}}-\overline{\left(\frac{\mu}{\sigma^{2}}\right)}\right]\pm\frac{1}{\lambda_ {1}}\left[\mu-\overline{\mu}\right]\right\}\]

The Gaussian model coordinate functions are depicted in Fig. 3B-E using a parameter distribution

\[\rho(\mu,\sigma)=\begin{cases}1/800\text{ if }-20\leq\mu\leq 20,\ 0<\sigma\leq 20 \\ 0\ \text{ otherwise}\end{cases}.\]

A projection of the model manifold onto the space-like and time-like coordinates corresponding to first moment of the model is depicted in Fig. 3G. The points on this projection are colored by the mean parameter \(\mu\). Here, we see a degree of rotational symmetry in the manifold projection, separated into negative mean values on the left and positive mean values on the right of the \(T^{+}_{1}\) center-line. Also note that there are apparent breaks in this manifold projection. These breaks do not reflect a true discontinuity in the structure of the manifold, but instead reflect the density with which the \((\mu,\sigma)\)-space is sampled. We will see manifold breaks related to the sampling density in our results for the network models.

## IV Results

Before proceeding with results, it is helpful to briefly summarize the goal of this paper and the workflow constructed in prior sections. We wish to study the population-averaged behavior of stochastic spiking models as we vary synaptic and membrane timescales, repeating this across a range of relative excitation and inhibition. To do this, we approximate the full spiking network dynamics as a population-averaged multivariate Gaussian process (Eqn. 6). We choose a sample of inverse timescales as discussed in Sec. II.5, and in particular constrain the sampling based on the stability condition for the corresponding non-spiking model (Eqn. 10). Within this sampled regime of timescales, the Gaussian process approximations should be mean-reverting and thus reach a stationary Gaussian distribution. We numerically solve for the vector-mean and the covariance matrix of the stationary Gaussian distribution at each sampled timescale-point. Finally, we embed this manifold of stationary Gaussian distributions into a behavioral space using the isKL methods introduced in Sec. III. With the workflow summarized, we may proceed.

### Gaussian process approximations are stable

A key step in the analysis workflow is to find the stationary distribution for the approximated processes at each sampled timescale-point. For the stationary distribution of an Ornstein-Uhlenbeck process to exist, all of the eigenvalues \(\lambda_{\mathbf{A}}\) of the drift matrix \(\mathbf{A}\) must have a positive real component. Basing the upper sampling boundary on the theoretical stability boundary of the related linear model, as well as the check for singularities in the drift matrices \(\mathbf{A}\) and covariance matrices \(\mathbf{C}\), should ensure this requirement is met. We confirm this
by explicitly examining the eigenvalues of the sampled models.

Each individual model--as specified by the model type, E/I log-ratio \(R\), and a pair of inverse timescales (\(\tau_{m}^{-1},\tau_{s}^{-1}\))--has three drift-matrix eigenvalues. We pool together the eigenvalues of all particular models on a given model manifold as specified just by the model type and E/I log-ratio \(R\). The resulting eigenvalue distributions for a subset of E/I conditions \(R\) are given in Fig. 4. Points in the eigenvalue distribution are colored by the log-distance of the particular model (specified by \((\tau_{m}^{-1},\tau_{s}^{-1})\)) from the origin in the inverse-timescale parameter space, shown in the insert. Each individual eigenvalue distribution is accompanied by marginal histograms where appropriate.

The subset of manifolds shown in Fig. 4 highlights a portion of the inhibition-dominated regime for which the models were observed to have complex eigenvalues. These complex distributions seen in the spiking-type model manifolds have a relatively small range in the imaginary direction and the imaginary components tend to pool near the origin in the along the real-axis. Additionally, most of the density for these complex distributions are along the real-axis itself, indicating that models with complex eigenvalues are relatively rare within their corresponding manifolds. The manifolds for the remaining E/I conditions have eigenvalue distributions qualitatively very similar to those at the extremes: purely real and positive eigenvalues spanning roughly the same range and skewed toward the origin. A key takeaway from Fig. 4 is that all of the eigenvalues have strictly-positive real components, confirming that the sampled models for each manifold are stable and therefore appropriate for embedding analysis.

### Behavior of full spiking network models

As discussed in Sec. III, the isKL embedding methods take the model manifold from the parameter space and position it an a hyperbolic using the symmetric Kullback-Liebler divergence. The KL-divergence functions similarly to a distance between models based on their (sufficient) statistics which determine the behavior of a particular model from the manifold. The isKL method thus embeds the model manifold in a behavioral space. This connection to the underlying behavior of the sampled models can be obscured when looking only at the results of the embedding analysis. As such, we will take some time here to discuss the behavior of the full-network model described by Eqn. 1.

Full-network spiking models were generated from the appropriate parameters in Table 1 with a sparse, random connection matrix as described in Sec. II.1. A subset of model manifolds were chosen from across the range of E/I conditions \(R\), and individual models from these manifolds were taken from along an arc of radius \(\sim 0.01\) (see e.g. Fig. 5, left column). The membrane and spiking dynamics described by Eqn. 1 with a specific choice of timescales were simulated using a basic forward-Euler integration scheme using a time step \(dt=0.1\) ms. Most models were simulated for \(20,000\) ms. The models in bottom 3 rows of the right-most column simulated for increased durations--\(100,000\) ms, \(150,000\) ms, and \(150,000\) ms respectively--to assure convergence to a stationary behavioral regime.

The long-term dynamics of the population-averaged membrane potentials and the membrane potentials of six excitatory and six inhibitory neurons are shown in Fig. 5. The mean-field values of the membrane potentials predicted by Eqn. 6c are represented by the red and blue dotted lines. From this figure, we see that the population-averaged membrane potentials in the full spiking network do indeed reach stationary values as predicted by the drift-eigenvalues \(\lambda_{A}\) for the approximated spiking models shown in Sec. IV.1. Additionally, we see that the mean-field membrane potential values correspond fairly well to the stationary population-averaged potentials (Fig. 5, columns 1 and 2), but this breaks down near the upper bound of the arc (Fig. 5, column 3). This upper boundary of the arc corresponds to the stability boundary in first two manifolds (Fig. 5, rows 1 and 2) and a bifurcation boundary in the remaining manifolds. The breakdown of the mean-field approximations at these limits thus lines up with the colloquial understanding of their accuracy. Knowing that the population-averaged potentials of the full spiking networks converge to a stationary condition, we next look at how these stationary solutions of these models differ.

The membrane potential dynamics in the stationary regime of these same models are shown in Fig. 6. Here, we show data from the last 1000 ms of simulation and ignore the mean-field predictions. We see that the population-averaged membrane potentials visibly fluctuate around a average variable for most of the simulated models. These fluctuations in the population averages become less noticeable as the overall magnitude of the averages and standard deviations increase, e.g. along column 3 of Fig. 6. A similar trend is seen in the membrane dynamics for individual neurons in the network. Fluctuations in individual potentials are very large relative to their mean values for the first two models along the arc (Fig. 6, columns 1 and 2). The population-variance in membrane potentials for these models is thus highly dependent on the fluctuations of individual membrane potentials. Towards the upper end of the arc (Fig. 6, column 3), the magnitude of the membrane potentials increases and the fluctuations of individual potentials are less pronounced. For these models, the population-variance of the membrane potentials is much more dependent on the spread of individuals around the population-mean as opposed to the fluctuations of those individuals.

We finish this section by examining the actual spiking dynamics for the example models in Fig. 7. Here we show raster plots for each model during the last 1000 ms of simulation. First, we see that the target neuron (black
spikes) fires more frequently than other excitatory neurons in the network, which is to be expected as it receives extra current input. For each model manifold (different rows in Fig. 7), we observe an overall increase in the rate of spiking in the network as we move along the arc from point 1 to point 3. This aligns with the change in overall magnitude of the membrane potentials seen in Figs. 5 & 6. This observation also aligns with an intuitive understanding of the timescales: along this arc, the relative rate of synaptic input becomes much faster than the relaxation dynamics. This trend is taken to the extreme in the models at the top of the arc for each manifold (Fig. 7, column 3) where we see unrealistically high spiking rates in the last three rows. With this, we've built an intuitive understanding of the differences in behavior across the model manifolds and across E/I conditions. We now move on to the embedding analysis for these models.

### Network embedding is hierarchical

It has been previously reported that biological models exhibit a hierarchy of sensitivities to different parameter combinations relative to some cost function on model behavior [51]. A similar hierarchical structure has been observed in the widths of model manifolds and the corresponding eigenvalues induced by a particular embedding, and a correlation between the widths and eigenvalues has also been noted [47, 52]. The current modeling and embedding differs from these prior cases in that we

Figure 4: **Eigenvalues of drift matrices** Eigenvalue distributions for the drift matrices \(\mathbf{A}\) of all sampled network models for the non-spiking model manifolds (top) and the spiking model (bottom) as the excitation-inhibition ratio \(R\) is adjusted. Points are colored by the radial distance \(d\) of model from the origin in inverse-timescale space, as illustrated by the inset. The marginal histograms for the real (\(\mathcal{RE}\), top) and imaginary (\(\mathcal{LM}\), right) are given for each distribution. Plots without a histogram in the imaginary dimension indicates a marginal \(\delta\)-distribution. Emphasis is placed on a portion of the inhibitory regime (\(-0.19<R<-0.09\), middle columns) in which some eigenvalue distributions display imaginary components.
are embedding probabilistic models in behavior space. Considering also the limited dimensionality of the current embedding, it is unclear if this hierarchical property should manifest in the current system. We show below that the manifolds for models of the types in Eqns. 6 & 9 are indeed hierarchical under the i\(\mathrm{k}\)KL embedding framework, with coordinate eigenvalues spanning several orders of magnitude for each E/I condition.

We used the i\(\mathrm{s}\)KL methods (Sec.III) to embed the stationary distributions for both the spiking and non-spiking model types across 25 E/I conditions ranging from the excitation-dominated to the inhibition-dominated, and approximately centered at \(R=0\). The root absolute eigenvalues for the embedding coordinates \(\{\Lambda_{i}^{\pm}\}\) are plot

Figure 5: **Membrane potential dynamics of spiking network models** The dynamic population-averaged membrane potentials (solid lines) are plotted against the predicted mean-field values (dotted lines) for three different parameter pairs across the examined range of E/I ratios \(R\). The membrane dynamics of six excitatory and six inhibitory neurons (dashed lines) are also shown for each condition. The ribbons around the population-averaged potentials are the the standard deviation of the membrane potentials within the corresponding population at each time point. The sampled parameter distribution from the embedding calculations, along with the chosen points for spiking simulation, are given in the right-most column. Full-network spiking simulations were run until an apparent stationary state was reached.
ted against the observed manifold width along the same coordinate for the non-spiking models (Fig. 8A) and the spiking models (Fig. 8C). Here, the manifold width is taken to be the simple range across a given coordinate. We see that the widths and eigenvalues are indeed correlated across E/I conditions for both model types, following with previous observations [47, 52].This suggests these two measures may be used interchangeably in further analysis. The coordinate eigenvalues of the non-spiking model (Fig. 8B) span at least three orders of magnitude for any given E/I condition tested, and up to nearly fifteen orders of magnitude at the most ex

Figure 6: **Stationary membrane potential dynamics of spiking network models** The dynamic population-averaged membrane potentials (solid lines) are plotted against the predicted mean-field values (dotted lines) for three different parameter pairs across the examined range of E/I ratios \(R\). The membrane dynamics of six excitatory and six inhibitory neurons (dashed lines) are also shown for each condition. The ribbons around the population-averaged potentials are the the standard deviation of the membrane potentials within the corresponding population at each time point. The sampled parameter distribution from the embedding calculations, along with the chosen points for spiking simulation, are given in the right-most column. Full-network spiking simulations were run until an apparent stationary state was reached, and the membrane dynamics for the the last 1000 ms of simulation time are plotted.
treme. The coordinate eigenvalues for the spiking models (Fig. 8D) span roughly two to three orders of magnitude on the extreme ends of the E/I spectrum and upwards of four in at some points in the center, with eigenvalues peaking towards the center as you approach from either extreme. Taken together, both model types studied here exhibit a hierarchical structure in line with prior observations of other systems, albeit with a more limited degree of separation in the case of the spiking-type models.

Before proceeding, we make some comparative observations between the two model categories. The scale and range of eigenvalues for the non-spiking model significantly larger than those for the spiking model type when in the excitation-dominated regime. Additionally, the non-spiking models show a sharp jump in eigenvalues when moving from the inhibition-dominated regime to the excitation-dominated one. This jump in eigenvalues may indicate a sort of bifurcation in the overall manifold. The eigenvalues \(\{\Lambda_{i}^{\pm}\}\) directly reflect the covariance and--anecdotally more importantly--the variance in the corresponding sufficient statistics and natural parameters. A jump in the magnitude of the eigenvalues thus

Figure 7: **Stationary spiking behavior of full network models** Example raster plots for individual time-scale pairs for spiking models across the sampled range of E/I ratios \(R\). The sampled parameter distribution from the embedding calculations, along with the chosen points for spiking simulation, are given in the right-most column. Full-network spiking simulations were run until an apparent stationary state was reached, and the spikes from the last 1000 ms are plotted. In each model, neuron index 350 was designated as the target neuron and its spikes are shown in black.
indicates a sudden increase in the variability of model behavior, and this could correspond to sampling near the stability boundary (Eqn. 10) in the case of the transition seen the non-spiking models. A similar transition may be happening at the peaks in the eigenvalue distributions of the spiking-type models, however these are much less pronounced than the one seen in the non-spiking models and the increase does not persist through the excitation-dominated regime as in the non-spiking models. We note that the firing rate non-linearity for the spiking model-type (Table 1) is quasi-linear when \(x\gg 1\). A naive prediction would be a similar behavior between model types when the membrane potentials become more positive as in the excitation-dominated regime. However, this is not reflected in the observed distributions of \(\{\Lambda_{\rm i}^{\pm}\}\).

### Projection hierarchies

Having established that the isKL embedding of the spiking models and the non-spiking models exhibit a hierarchical structure, we next want to interrogate this structure across our model manifolds. We do this by examining projections of the manifolds onto lower-dimensional spaces along the largest widths and smallest widths. We will focus on 2-dimensional projections.

Fig. 9 shows the largest manifold projections in behavioral space for the non-spiking models and spiking models across a subset of E/I conditions. Points on these manifolds are colored by the mean value of the membrane potential for the test neuron \(\langle V_{0}\rangle\). It is visually clear that the manifold projections are shrinking from top left to bottom right for each condition, reflecting the hierarchical structure of the manifold. A large fraction of projections--for example, Fig. 9 column 3, row 2--have apparent gaps in their structure. These are similar to the gaps seen in the projection of the example Gaussian distribution shown in Fig. 3G, and are tied to the sampling density used across the inverse-timescale space near key boundaries (data not shown).

Many of the projections across model types and E/I conditions appear very linear or piece-wise linear, for example Fig. 9 column 2 rows 1-4. This thinness at the largest scales would suggest a relatively simple relationship between the largest coordinates and that the model manifold is relatively flat. The difficulty of overcoming under-sampling of the parameter space complicates this interpretation slightly. The gaps in the projections seen in the excitation-dominated regime are clear evidence of some under-sampling, but interpolating the data across gaps suggests that the projections in these conditions may still be piecewise- or quasi-linear. These stick-like projections both model-types in the excitation-dominated regime. The projections are also seen to change shape qualitatively as the E/I conditions change. In the non-spiking models, we see the appearance of spoon-shaped projections as we move into the inhibition dominated regime. In contrast we see knife-like projections in the spiking models, albeit only under the most inhibitory of E/I conditions. This qualitative change in the manifold projections seen across the two model types could be caused either by warping of the manifold along each coordinate as the conditions change or by changes in ranking of each coordinate. This point will be revisited in Sections IV.5 & IV.7.

Additionally, we note that many of the projections separate points on the manifold in alignment with \(\langle V_{0}\rangle\), as was seen in the example embeddings shown in Fig. 3F & G. This is particularly clear, for example, in the inhibition-dominated regime of the two model types (Fig. 9 row 5). The separation of the manifold into sections based on behavioral regimes depends on more than just \(\langle V_{0}\rangle\), however. For example, we see no such trend in Fig. 9 column 2 row 4. In this case the individual models within the manifold are more significantly separated by (a combination of) behavioral parameters that, in a sense, have a dependence on the timescale parameters that is orthogonal to the way \(\langle V_{0}\rangle\) depends on them. Exceptions aside, this noted \(\langle V_{0}\rangle\)-aligned separation serves as a clear demonstration of the behavioral clustering induced by the isKL method.

We can also examine the smallest projections of the model manifolds for the two model types, which correspond to the least important modes of the expansion in Eqn. 17. Fig. 10 shows the smallest projections of the model manifolds for all three model types across E/I conditions. As was the case for the largest coordinate projections, there is evidence of under-sampling in the smallest projections also. This particularly evident in the non-spiking model manifold in the excitation-dominated regime (Fig. 10, column 1, rows 1 and 2). In contrast to the largest projections, the stick-like projections comprise the minority of the small-projection shapes. The model manifolds appear instead to be highly curved at the fine-grained level. Following the observation from the large-scale projections, we see the smallest manifold projections separate points in line with \(\langle V_{0}\rangle\). In particular, the counterexample mentioned before ( Fig. 9 column 2 row 4) now also shows a degree of alignment with changes in \(\langle V_{0}\rangle\), now in Fig. 10 column 2 row 4. This highlights the fact that model separation along different directions on the manifold can be more or less tied to a particular behavioral or statistical parameter.

Before moving on, we should highlight the relationship between the projection coordinates and model behavior. Coordinates with larger eigenvalues contribute more to separation between \(p(x|\theta)\) and \(p(x|\theta^{\prime})\) as measured by the \(D_{sKL}\). More specifically, a relatively large absolute-eigenvalue \(|\Lambda_{i}^{\pm}|\) indicates that the corresponding natural parameter is a relatively better way to separate individual models on the manifold by their behavioral predictions, or alternatively that a larger part of the variance of behavioral predictions across the manifold are explained by the associated natural parameter. Given a tractable mapping between the underlying model parameters and the natural parameters, the magnitude of co
ordinate eigenvalues can also give a sense of the relative importance of different parameter combinations aligning with a given coordinate direction. Lastly, the high degree of correlation between coordinate eigenvalues and the manifold widths along those coordinates means that the relative size of a particular projection gives a visual representation of the importance of the corresponding parameter combination.

### Coordinate rankings

We saw in Sec. IV.4 that the projections hierarchies of the model manifolds changed across the examined range of E/I conditions \(R\). One possible explanation for this is that the rankings of coordinates by manifold width change with \(R\). This potential aspect of the changing projections can be interrogated by tracking the rank of each coordinate across the range of \(R\). This will additionally provide insight into what aspects of the statistical model have the greatest (or least) impact on the overall model behavior for both model-types. Fig. 11 depicts the ranking for each coordinate by both the length of the manifold along that coordinate (top row) and eigenvalue-magnitude (middle row) for the spiking and non-spiking model types for a subset of E/I conditions. As each coordinate corresponds directly to a single sufficient statistic, we color- and shape-code the rank of each coordinate according to this correspondence. The bottom row of Fig. 11 reproduces the eigenvalue distributions shown in Fig. 8, except each point is now color- and shape-coded according to the sufficient statistic instead of the E/I measure \(R\).

We see in Fig. 11 that the ranking of coordinates changes across E/I conditions for both the non-spiking and the spiking models. We also note that the rankings by observed width (top row) and by eigenvalue-magnitude (middle row) agree fairly well across the range of \(R\). This agreement between the two sets of rankings makes sense when considering the high degree of correlation between the coordinate eigenvalues and manifold widths shown in Fig. 8. It is interesting that the

Figure 8: **Hierarchies of manifold widths** Top row: The correlation between the coordinate eigenvalues \(\Lambda_{\hat{x}}^{\pm}\) and the width \(W\) across the manifold in that coordinate direction as the E/I ratio is varied for the non-spiking models (A) and spiking models (C). A simple linear regression is applied to the log-widths and log-absolute eigenvalues for visualization. The distribution of the scale of coordinate eigenvalues as the E/I ratio is varied for the non-spiking models (B) and the spiking models (D).
eigenvalue distributions--particularly those for the spiking models--seem to separate into clusters of coordinates that do not intersect for much of the range of \(R\). For example, the tan-orange-steel blue (eight point star-five point star-hexagon) cluster at the top of the eigenvalue spectrum corresponds to the second moments involving the inhibitory and excitatory populations. This cluster remains consistently above the grey-magenta (downward triangle-pentagon) and pink-red (diamond-square) clusters--which correspond to the second moments involving the target neuron and first moments for the bulk populations, respectively--across \(R\) for the spiking models.

Knowing that the eigenvalues and thus the eigenvalue-magnitude rankings form these clusters across the E/I spectrum, it is natural to examine the sufficient statistics that correspond to coordinates in these clusters. We will focus on the spiking-type models that exhibit these

Figure 9: **Largest manifold projections** Projections of the manifolds for the two model types onto the largest four coordinates as determined by the observed manifold width. The hierarchy of projections is shown as the excitation-inhibition ratio is changed from excitation dominant regime (top) to an inhibition-dominant regime (bottom). Manifolds are colored by \(\langle V_{0}\rangle\)—the mean membrane potential of the test neuron—and each projection is scaled by the largest observed width. These projections are the “stiffest,” contributing the most to the behavior of the distribution of activities.
clusters. The tan-orange-steel blue (eight point star-five point star-hexagon) cluster noted before dominates over other clusters in the spiking models, and these coordinates correspond to the second-order moments of the membrane potentials of the excitatory and inhibitory populations (Fig. 11 columns 2 and 3). This indicates that they are the most important statistics for distinguishing models across the manifolds. These second moments are also important for the non-spiking model manifolds, but they only sit at the top of the hierarchy in the excitation dominated regime \(R>0\). The grey-magenta (downward triangle-pentagon) cluster in the spiking models corresponds to the mixed second moments involving the target neuron \(V_{0}\) and either \(V_{1}\) or \(V_{0}\). This cluster is above the green (triangles) cluster in the mid-range of \(R\) and just below it in the extreme E/I conditions, and this green (triangles) cluster is the pure second moment \(\langle V_{0}\rangle\). The green cluster is gener

Figure 10: **Smallest manifold projections** Projections of the model manifold for the two model types onto the smallest four coordinates as determined by the observed manifold width. The hierarchy of projections is shown as the excitation-inhibition ratio is changed from excitation dominant regime (top) to an inhibition-dominant regime (bottom). Manifolds are colored by \(\langle V_{0}\rangle\), and each projection is scaled by the observed width along the coordinate ranked \((M-3)=15\). These projections are the “sloppiest,” contributing the least to the behavior of the distribution of activities.
ally above the pink-red (diamonds-squares) and the blue (circles) clusters, except for a brief crossing of the pink-red and the green clusters around \(R\approx-0.18\). These last two clusters correspond to the mean values of all three membrane potentials. Taken together, these observations say that for the spiking-type models the fluctuations are more important for distinguishing between individual models on a given manifold and--for both the first and second moments--the statistics that involve the test neuron are generally less important than those that do not. These observations hold for the non-spiking model manifolds in the excitation-dominated regime, but not in the inhibition-dominated regime (Fig. 11, column 1).

Let us discuss the coordinate rankings at a more granular level of detail. While both model types have the second moments at the top of their respective hierarchies in the excitation-dominated regime, it is interesting to note how they differ here. The spiking-type models have the \(\langle V_{2}^{2}\rangle^{+}\)-related coordinates at the top while the non-spiking model is topped by the \(\langle V_{1}^{2}\rangle^{+}\)-related coordinates. The suggests that the degree of fluctuations in the inhibitory population are the most varied for the spiking models in this regime, but the excitatory population fluctuations take that title in the excitation-dominated non-spiking models. The last fine-grained detail we highlight here is the increased importance of the \(V_{0}\)-moments in distinguishing the behavior of the inhibition-dominated non-spiking models relative to their importance in the inhibition-dominated spiking-type models.

In addition to visualizing the relative importance of certain sufficient statistics for distinguishing between particular models across a given model manifold, we get another piece of information visualized for free through the eigenvalue-magnitude ranking plots (Fig. 11, middle row). Recall from Eqn. 16 that the eigenvalues \(\Lambda_{i}^{\pm}\) are given by the covariance of the sufficient statistic and natural parameter across the manifold (\(\text{Cov}(\eta_{i},\langle t_{i}\rangle)\)) and the geometric means of their individual variances (\(\sqrt{\text{var}(\eta_{i})\text{var}(\langle t_{i}\rangle)}\)). As we know the \(\Lambda_{i}^{-}\) eigenvalues are negative and of the same order of magnitude as the corresponding \(\Lambda_{i}^{+}\) (Fig. 11, row 3), we know that geometric mean of those variances greatly outweighs their covariance. Furthermore, the relative ranking of \(\Lambda_{i}^{+}\) and \(\Lambda_{i}^{-}\) imply the sign of the covariance \(\text{Cov}(\eta_{i},\langle t_{i}\rangle)\): if \(\Lambda_{i}^{+}>\Lambda_{i}^{-}\) then \(\text{Cov}(\eta_{i},\langle t_{i}\rangle)>0\) and _vice versa_. For example, by looking at the eigenvalue-magnitude ranking of the steel blue (hexagon) coordinates in the spiking-type models (Fig. 11, row 2 columns 2 and 3) we see that \(\text{Cov}(-\frac{1}{2}C_{12}^{-1},\langle V_{1}V_{2}\rangle)<0\) across all E/I conditions examined here. While this could very easily be determined by looking at these covariances themselves--and they must be calculated in order to determine \(\Lambda_{i}^{\pm}\)--it is convenient to be able to glean this from a plot already produced for another purpose.

We will briefly summarize. Fig. 11 shows that the coordinate rankings do change across the sampled E/I range, thus explaining the changing projection hierarchies in Figs. 9,10 at least in part. We found that the coordinates form clusters in the eigenvalue distribution that behave in a correlated manner across the E/I spectrum and with which they share relations to similar types sufficient statistics. In particular, the cluster of coordinates corresponding to the fluctuations \(\langle V_{I}V_{J}\rangle^{\pm}\) for \(I,J\in\{1,2\}\) have the most impact on the activity of the spiking-type models, as well as in the excitation-dominated non-spiking models. We made observations of which types of fluctuations were most important to model distinction across the manifold for different model-types and different E/I conditions. Finally, we highlighted a secondary visual interpretation of the eigenvalue-magnitude ranking plots relating to the base statistical model.

### Transforming of base parameters

We highlighted in Sec. IV.5 that the statistical parameters from the stationary Gaussian distribution of membrane potentials--discussed in terms of the sufficient statistics--have a hierarchical impact on the possible behaviors of the spiking model that changes across the E/I spectrum \(R\). Further, we identified clusters of parameters that tended change in similar ways with \(R\). In the case of the spiking models, the fluctuations \(\langle V_{I}V_{J}\rangle^{\pm}\) for \(I,J\in\{1,2\}\) were the most impactful while the mean coordinates \(\langle V_{0}\rangle^{\pm}\) had a relatively small impact. While important, these observations do not directly address the role of the inverse timescales (\(\tau_{m}^{-1},\tau_{s}^{-1}\)) on model behavior. Unfortunately, the mapping from the timescale parameters to statistical parameters is intractable, owing primarily to the transcendental system of mean-field equations (Eqn. 6c). Closed forms for the stationary distribution parameters of the linear non-spiking models can be found, but these expressions are ratios of high-degree polynomial functions of the timescales and do not directly reflect the mapping in the spiking model context. To begin untangling the impact of the timescale parameters on the range of model behaviors, we must thus rely on a qualitative understanding of the relationship between the timescales and e.g. the sufficient statistics.

In Fig. 12, we plot the logarithm of several sufficient statistics as a function of position in inverse-timescale space for the spiking network with \(R=-6.06\times 10^{-3}\). We include \(\langle V_{1}^{2}\rangle\) (Fig. 12,B) and \(\langle V_{2}^{2}\rangle\) (Fig. 12,C) from the upper cluster as well \(\langle V_{0}\rangle\) (Fig. 12,A) to serve as a representative set from across the hierarchies in Fig. 11. Note that the second moments have different units than those of \(\langle V_{0}\rangle\), which should be kept in mind when comparing the color scales. That said, the \(D_{sKL}\) between two members of the same exponential family can be rewritten as [53]

\[D_{sKL}(\theta,\theta^{\prime})= \sum_{i}\left(\eta_{i}(\theta)-\eta_{i}(\theta^{\prime})\right) \left(\langle t_{i}(\theta)\rangle-\langle t_{i}(\theta^{\prime})\rangle \right).\]

Paired with the sufficient statistics of a multivariate normal distribution (Eqn. 18), we see that the \(D_{sKL}\) is in
some sense weighing the first and second moments directly against each other. This in mind, the variability in the second moments is \(\sim 2\) orders of magnitude larger than that for the mean of the test neuron, in line with their relative ranking in Fig. 11. Further, we note the similar dependence of all three statistical parameters on the inverse timescales, increasing in magnitude radially from \(\tau_{m}^{-1}\)-axis to the \(\tau_{s}^{-1}\)-axis as well as exhibiting a "cold spot" triangle on the right-most corner of the sampled wedge. The trends between the means and covariances differ most significantly along the \(\tau_{s}^{-1}\) boundary. Here, the magnitude of the mean increases towards the inverse-timescale origin (i.e. very long timescales) while the second moments increase moving away from the origin (i.e. very short timescales).

The presence of the cold spot in the mappings to the statistical parameters--particularly the sharpness of the transition seen for \(\langle V_{0}\rangle^{\pm}\)--reinforce the intuition that translating changes in the statistical parameters back onto the timescale parameters is non-trivial. That said, the shared general trend in the mappings suggest a possible avenue for model reduction if some loss of expressivity is permitted. Reducing the 2-dimensional sampled space to an arc around the origin and through the cold spot

Figure 11: **Coordinate rankings** The ranking of each manifold coordinate \(T_{i}^{\pm}\) as the E/I balance is changed in both model types. Coordinates are ranked from most important (top of each plot) to least important (bottom of each plot) based on the observed width of the manifold along said coordinate (top row) or the magnitude of the corresponding eigenvalue \(|\Lambda_{i}^{\pm}|\) (middle row). The log-magnitude of the eigenvalue for each coordinate is given in the bottom plot as the E/I balance \(R\) is changed as in Fig. 8. The legend renames each coordinate \(T_{i}^{\pm}\) to the corresponding sufficient statistic \(\langle t_{i}\rangle^{\pm}\) for ready interpretation.
could be used to capture the concomitant increases in the magnitude of the first and second moments, capturing the majority of their respective variability. Alternatively, radial sampling along the \(\tau_{s}^{-1}\) boundary could be used to study the apparent trade-off in magnitude of the means and covariances. This idea of model reduction is intimately tied to notions of model dimensionality, which we will return to in Sec. IV.8 and Sec. V.

### Manifold projections change smoothly with E/I balance

We now return to the question raised at the beginning of Sec. IV.5: What causes the projection hierarchies to change across E/I conditions? While the changing coordinate ranking observed across E/I conditions for both models can explain the changing manifold projections, it does not rule out the possibility that the projections along a given coordinate are themselves changing. To address this possibility, we project the model manifolds for both of the model-types onto the same pair of coordinates across the E/I spectrum in Fig. 13. We chose to project the manifold onto the space-like \(\langle V_{0}^{2}\rangle^{+}\) and \(\langle V_{0}\rangle^{+}\) coordinates as the statistical behavior of the test neuron may be of particular interest in some scenarios. For the sake of visualization, each projection along the E/I spectrum is scaled by the larger of the two manifold projections at each condition. The overall scale of the projection is given by the axis scale.

We can see in the projections of the non-spiking model manifolds (Fig. 13, upper diagonal) that there is a squashing and stretching of the manifold relative to the overall change in scaling as the E/I conditions are changed. Additionally, these transformations appear to act smoothly on the manifold projections until the manifold flattens going from the inhibition-dominated regime to the excitation-dominated one in the range \(-6.06\times 10^{-3}\leq R\leq 0.115\). This flattening reflects a radical increase in the manifold scale along the \(\langle V_{0}^{2}\rangle^{+}\)-coordinate relative to the \(\langle V_{0}\rangle^{+}\)-coordinate as all of the eigenvalues are seen to jump (Fig. 11, column 1 row 3). This interpretation is corroborated by the change in overall scale of the axes--from \(\sim\mathcal{O}(10^{3})\) for \(R<0.115\) to \(\sim\mathcal{O}(10^{17})\) for \(R\geq 0.115\) (Fig. 13, upper diagonal)--and the correlation between eigenvalue-magnitude and manifold width discussed in Sec. IV.3 (Fig. 8). The projections of the spiking manifolds (Fig. 13, lower diagonal) are also seen to transform smoothly with \(R\) with the fork-shaped projections (e.g. lower diagonal, \(R=-0.307\)) collapsing into the spoon projections (e.g. lower diagonal, \(R=-0.103\)) seen in the non-spiking model around \(R\approx-0.1\). The projections for the spiking-type manifolds do change along \(R\) in line with the changes in their respective eigenvalue distributions (Fig. 11, columns 2 and 3, row 3), but these changes are more subtle than in the non-spiking model. The eigenvalue distributions for the spiking-type models drift downwards as you move from \(R<-0.18\) to \(R>0.10\), and this is mirrored in the manifold projects by a slight decrease in projection scale moving in the same direction.

We have shown here that the manifold projections for both model-types do indeed change across the sampled E/I range, which plays a subsequent role in the changing of projection hierarchies across E/I conditions noted

Figure 12: **Mapping from inverse timescales to statistical parameters** The relationship between the inverse timescales \((\tau_{m}^{-1},\tau_{s}^{-1})\) and a select set of statistical parameters from the corresponding stationary Gaussian distribution for the spiking network with \(R=-6.06\times 10^{-3}\). Inverse timescale spaces are colored by the log-value of one of the sufficient statistics: A) \(\log_{10}\langle V_{0}\rangle\); B) \(\log_{10}\langle V_{1}^{2}\rangle\); C) \(\log_{10}\langle V_{2}^{2}\rangle\)
in Sec. IV.4. The individual projections were shown to undergo potentially significant rescaling across values of \(R\) that alter it visually, as noted in the non-spiking manifold. Additionally, the manifold can exhibit a warping, as in the fork-spoon-fork transition noted in the spiking-type models.

### Manifold dimensionality

As mentioned at the outset (Sec. I), a key issue when analyzing the behavior of a collection of large spiking networks is the dimensionality of the behavioral output space. The true behavioral space of the model--as expressed through the spiking activity--grows with both an increasing network size \(N\) and a decreasing time bin size \(\Delta t\). A goal of the current work is to understand the behavioral output of these models in a lower dimensional framing. Thus, we will briefly interrogate the dimensionality of our spiking and non-spiking models before proceeding to the final discussion.

The behavioral dimensionality of models following Eqn. 1 is \(NT/\Delta t\) for a discrete-time trial of length \(T\). The dimensionality of the behavioral space remains the same when moving to the Gaussian process approximation of the model in Eqn. 5. Population-averaging of

Figure 13: **Coordinate evolution as E/I balance is tuned** Projections onto the \(\langle V_{0}^{2}\rangle^{+}\)-\(\langle V_{0}\rangle^{+}\) plane of the non-spiking model manifolds (top) and the spiking model manifolds (bottom) as the excitation-inhibition ratio is adjusted. The axes of each 2-D projection are scaled to the larger width for the given model and E/I condition for the convenience of visualization. Note that the projections change in size by an order of magnitude as the E/I balance is adjusted. The projections are colored by \(\langle V_{0}\rangle\) as in prior figures.
these approximated dynamics into Eqn. 6 decreases the dimensionality to \(N_{\rm pop}T/\Delta t\), where \(N_{\rm pop}<N\) is the number of populations being considered. By simplifying our analysis to studying the _stationary distribution_ of population behaviors, the behavioral dimensionality drops to \(N_{\rm pop}(N_{\rm pop}+3)/2\) corresponding to the maximal number of independent sufficient statistics (see Sec. III). Finally, the isKL methods embed this manifold of behaviors into an \(N_{\rm pop}(N_{\rm pop}+3)\)-dimensional space which determines the upper limit of dimensionality that may be measured from embedded data (i.e. sampled models).

Having considered how the dimensionality changes across the steps for our analysis, two key questions remain. If we can only see the results of the embedding, how do we gauge the dimensionality of the manifold being analyzed? If we instead have an understanding of the maximal dimensionality of the system, is there any effective reduction in dimensionality that we can measure? To address these questions, we start with a measure of effective manifold dimensionality commonly used in principal component analysis (PCA) known as the participation ratio (PR):

\[PR=\frac{\left(\sum_{i}\Lambda_{i}\right)^{2}}{\sum_{i}\Lambda_{i}^{2}}.\]

As the isKL embedding methods are intimately tied to multidimensional scaling (MDS)--an extension of PCA--the PR should serve as a useful base for measuring the effective dimensionality of our embedded model manifolds. This is complicated slightly by the presence of negative eigenvalues \(\{\Lambda_{i}^{-}\}\) that arise in MDS, so we use an altered PR as our measure of effective dimensionality:

\[PR=\frac{\left(\sum_{i,\pm}|\Lambda_{i}^{\pm}|\right)^{2}}{\sum_{i,\pm}\left( \Lambda_{i}^{\pm}\right)^{2}}. \tag{19}\]

The effective dimensionality of our two model-types across the examined E/I spectrum is shown in Fig. 14. We see that the spiking-type models begin with a relatively high \(PR\approx 8\) in the inhibition-dominated regime before dropping to \(PR\approx 6\) in the middle regime and then rising slightly again in the excitation-dominated regime. By contrast, the non-spiking model has \(PR\approx 3\) in the excitation-dominated regime. The \(PR\) for the non-spiking model then peak at \(PR\approx 6\) around \(R\approx-0.1\) before decaying back down further into the inhibition-dominated regime. The non-spiking models thus have a lower effective dimensionality than the spiking-type models.

How do we contextualize the measured \(PR\) for these models? First, we note that the maximal possible measured dimensionality for both the spiking and non-spiking model-types is \(N_{\rm pop}(N_{\rm pop}+3)=18\), and the _statistical model_ dimensionality is \(N_{\rm pop}(N_{\rm pop}+3)/2=9\). This indicates that the approximate models show a dimensionality reduction compared to both the model dimensionality and the maximal embedding dimensionality. This seems trivial until one examines the \(PR\) measure for the example embeddings given in Sec. III. The simple Poisson example has just 1 parameter \(\Lambda\) and subsequently a maximum embedding dimension of 2. Despite the intrinsic parameter density of 1, its isKL embedding (Fig. 3F) gives an effective dimensionality is much closer to its maximal embedding dimension and gives \(PR\approx 1.982\). In a similar vein, the example Gaussian model has two parameters \((\mu,\sigma)\), yet it has an embedding dimension of \(PR\approx 3.929\) which is nearly its maximum embedding dimensionality of 4 (see Fig. 3G for one of the manifolds projections). The measured \(PR\) thus does not seem to reflect the dimensionality of the intrinsic manifold structure, but instead the number of embedding dimensions within the isKL framework required to capture most of the variability in model behaviors. This will be discussed further in Sec. V.

## V Discussion

The central motivation of this paper is to tease apart the impact of cellular and synaptic model parameters--internal timescales and relative synaptic strengths, respectively--on the complex and high-dimensional behavioral space of spiking network models. Taking inspiration from prior applications of information geometry to neural systems [25; 26; 27; 28; 29; 30; 31; 32], we approached this Herculean task by leveraging recently developed methods for studying the information geometry of complex biology models [52; 53] and applying them to spiking network models with more biological features than those considered previously. We began by defining our spiking model [17; 54] and then simplifying it through population-averaging, using a path-integral formalism to approximate the membrane dynamics as a Gaussian process [17], and then calculating the stationary distribution for that approxima

Figure 14: **Dimensionality of model manifolds** The effective dimensionality of non-spiking and spiking models across the spectrum of E/I conditions is given. Effective dimensionality is given by the altered participation ratio given by Eqn. 19.
tion [57]. The stationary distributions for these were then analyzed using the information geometric framework introduced by Teoh and colleagues [53]. This workflow is the core of the work presented.

Before diving into the results of the geometric embedding analysis, we briefly examined the behaviors of full spiking networks across various E/I conditions and for a few different timescale points. We showed that the behavior of the actual networks change distinctly across the variables at both the level of spiking and of population-averaged membrane dynamics. Importantly, the spiking models reach a stationary behavior in the long-time limit. This agreed qualitatively with the mean field predictions and supported the analysis of the stationary distributions from the reduced model.

The information-geometric analysis demonstrated that the approximated models are hierarchical. Manifold widths and coordinate eigenvalues spanned several orders of magnitudes, pointing to a "hyperribon" structures with "stiff" and "sloppy" coordinate directions. The distribution of these coordinate eigenvalues changed across E/I conditions and with it the hierarchy of 2-dimensional manifold projections. These changes in the manifold projections arose from a smooth warping of projections onto specific coordinate pairs as well as a reordering of the coordinate rankings. Identifying each coordinates with their corresponding sufficient statistic highlighted a clustered structure in the eigenvalue distribution of the spiking models across E/I conditions. From this clustered structure, it is possible to pick out the most and least important sufficient statistics for distinguishing between models on a given manifold--these are the parameter combinations that underlie the stiff and sloppy coordinate directions, respectively. In particular, the stiffest directions on the spiking-model manifolds corresponded to the second moments of the excitatory and inhibitory population membrane potentials while the sloppiest directions were those corresponding to the first moments. This suggests that bulk fluctuations are key for determining the behavior of a specific network. It is unfortunately difficult to tie this understanding of stiff and sloppy statistical parameters to the timescale parameters in a manner that is satisfactorily analytical, owing primarily to the transcendental mean-field equations (Eqn. 6c) that must be solved numerically. That said, the implication of the sloppy and stiff coordinate observations is that an adjustment of the membrane and synaptic timescales tends to have a larger effect on the large-population fluctuations than it does on the means.

At the end of our isKL analysis, we began a discussion regarding the dimensionality of the models, their behavior, and their manifolds. The largest reductions in the size of the model being discussed occur when moving to population-averaged models and when focusing on the stationary distribution. The combined effect decreases the dimensionality of the behavioral space being studied from \(NT/\Delta t\) to \(N_{\rm pop}(N_{\rm pop}+3)/2\), in which we essentially shift from a study of particular spike patterns to a study of probability distributions. From here, the isKL methods embed the distribution in \(N_{\rm pop}(N_{\rm pop}+3)\) dimensions. This sets the upper limit of dimensionality at the end of our analysis, that upper limit being 18 for the particular architectures studied here. Using an altered participation ratio to measure the effective dimensionality of our embedded spiking models gave us a range of \(6\lessapprox PR\lessapprox 8\)--depending on the E/I measure \(R\)--less than the maximum possible dimension. It was illustrated elegantly through the example embedding of the 1-dimensional Poisson model that the participation ratio measures the number of dimensions needed to hold a sufficiently representative version of the model manifold rather than the intrinsic dimensionality of the manifold. In fact, the participation ratio for both toy models indicated that they basically "filled" their respective embedding spaces. Taken together, these show that approximated spiking models are definitely undergoing a degree of dimensionality reduction as they are not filling the embedding space like the toy models did. The participation ratio can then be interpreted as giving a sense of how "pointed" a change in parameters is. If a 6-dimensional space is needed to represent most features of the manifold, this likely implies that the modulated parameters are mostly affecting 3 natural parameters. However, this is a measure of the effect of base parameter (i.e. \(\tau_{m}\) and \(\tau_{s}\)) on behavior, and not necessarily reflect a minimal structure in the base parameter space need for nearly-full expressivity of the model.

From the copious stick-like projections seen in the hierarchies (Figs. 9,10), we may intuit that the embedded manifolds are of an even smaller dimension than is represented by the participation ratio. We can take this a step further by understanding the entire embedding process as a transformation of a manifold originally in the parameter space, implying that it should intrinsically be, at most, 2-dimensional. We also noted in Sec. IV.6 that there are _ad hoc_ ways of reducing the parameter space to a 1-dimensional curve while seemingly preserving much of the variability in statistical parameters. If the goal is to find a reduced number of base parameters to approximately cover the manifold in a more principled way, this would likely require estimating the intrinsic dimensionality of the model manifolds with more sophisticated tools than those discussed here. This would provide a number of parameters--or parameter combinations--needed to understand and express the model. Thus, a combination of both an intrinsic measure and the participation ratio provides a complimentary understanding of model manifold dimensionality through the lenses of necessary base parameters and range of impact, respectively.

Lastly, the properties of these spiking model manifolds is likely to be impacted by the conditions under which it is being studied, more specifically any particular task in which it is being implemented. The models studied here are functionally in a spontaneous regime with a tonic drive that is minor in the scale of the network. The structure of a given task is known collapse high-dimensional
spontaneous activity into a lower-dimensional behavioral space [4; 10], which might be seen directly in information-geometric interrogations such as the one performed in this paper. Furthermore, this may well affect which statistical parameters are important, in turn changing the coordinate rankings, projection hierarchies, and potentially even the degree to which the resulting manifolds are hierarchical. These possibilities require their own attention in follow-up work.

### Limitations

It is important to discuss the limitations of the framework of modeling and analysis expounded upon in this paper. The primary hurdle to expanded usage of these methods is that the base calculations required for each step combined with the number of samples needed to visually resolve the embedded manifolds make it costly to increase the dimensionality of the parameter space or the number of network elements. The manifolds embedded here required a large number of sampled parameter points to resolve adequately, restricting the number of parameters considered. Similarly, calculating \(\sim N^{2}\) statistical parameters under the Gaussian process approximately would be computationally infeasible and high intractable. The first of these restrictions led to the choice of only two key parameters--the timescales--in the current work. The second restriction motivates the reduction of the model by population-averaging. The embedding of the inverse timescale sub-plane (see Sec. II.5) revealed that much of the manifold was comprised of points near the boundaries where the behavior became pathological (data not shown). This suggests that a principled or data-informed restriction of parameter space may lead to decrease in the necessary per-parameter sampling density and ease the restrictions presented here.

### Future directions

We conclude by discussing future directions for this work. As developed here, the methods discussed could be applied to models of particular neural circuits in the brain to understand their stationary behavior. For example, one could study the range of behavior of a cortical column when its internal timescales are subject to change. Further, one could study the conditioned range of behaviors in a network in response to a well-defined distribution of inputs as the statistics of the input distribution change. This latter example is meant to demonstrate that the general framework--marginalize, approximate, population-average, and then embed--can apply to modulated parameters other than those presented here.

Perhaps more interesting are the possible extensions of the methods themselves. Of primary interest is the extension of the iskL embedding methods to non-stationary systems. A first-pass way to do this would be to discretize time, apply the embedding procedure at each time-step, and trace points through the embedding space. While conceptually straightforward, this approach would involve significantly more computational investment and the interpretation of the results would be more complicated than in the system discussed in this paper. One could instead try to extend the iSKL embedding framework to apply directly to the path integral representations used to derive the Gaussian process approximations. This would require a proof that the desired properties of the iSKL embedding still hold in this functional context, a potentially harder barrier to clear. Together, these highlight the care with which these conceptual extensions of the current method must be carried out.

## Appendix A Population averaging of the Gaussian process approximated network

Here, we derive the reduced model from Eqn. 6 by first making a Gaussian process approximation on the full-network spiking model and then averaging the resulting dynamics by population. The stochastically-spiking full network, modeled using a nonlinear Hawkes process, is reproduced here:

\[\frac{dV_{i}}{dt}=-\tau_{m}^{-1}(V_{i}-\varepsilon_{I})+I_{i}+\tau_{s}^{-1} \left(\mu_{\text{ext}}-J_{\text{self}}\dot{n}_{i}(t)+\sum_{j}w_{ij}\dot{n}_{j }(t)\right)\]

\[\dot{n}_{i}(t)dt\sim\text{Poiss}[\phi(V_{i}(t))dt]\]

. Recall that the lowercase subscripts (\(i\), \(j\), etc.) denote individual neurons within the network. \(V_{i}\) is the membrane potential of neuron \(i\), \(\varepsilon_{i}\) is the leak reversal potential, \(w_{ij}\) is the strength of a synaptic connection from neuron \(j\) to neuron \(i\), and \(-J_{\text{self}}\) is an inhibitory self-coupling. The two currents \(\mu_{\text{ext}}\) and \(I_{i}\) represent an average current from an external network and an experimentally injected current, respectively. The process \(\dot{n}_{i}(t)\) is the spike train of neuron \(i\), and \(\phi(\cdot)dt\) is the instantaneous firing rate nonlinearity, here given by \(\phi(x)=\frac{1}{2}(x+\sqrt{x^{2}+1/2})\). Finally, \(\tau_{m}\) and
\(\tau_{s}\) are modulated membrane and synaptic timescales, respectively. The mean-field equations for the steady state membrane potentials can be obtained directly from these equations by using the fact that the approximation neglects fluctuations, and hence \(\langle n_{i}(t)\rangle=\langle\phi(V_{i}(t))\rangle\approx\phi(\langle V_{i}(t)\rangle)\), yielding

\[V_{i}^{\rm mf}=\varepsilon_{I}+\tau_{m}I_{i}+\frac{\tau_{m}}{\tau_{s}}\left( \mu_{\rm ext}-J_{\rm self}\phi(V_{i}^{\rm mf})+\sum_{j}w_{ij}\phi(V_{j}^{\rm mf })\right). \tag{100}\]

To obtain the dynamics of fluctuations around the mean-field predictions, and to set up for future calculations that go even beyond the Gaussian approximation, it is useful to introduce a path integral representation of this stochastic process, using techniques from statistical physics [55]. In discrete time, we can write the joint probability for the membrane potential \(\mathbf{V}(t)\) and the spike trains \(\dot{\mathbf{n}}(t)\) as follows:

\[P[\mathbf{V}(t),\dot{\mathbf{n}}(t)]=\prod_{t,i}P[V_{i}(t)|\dot{\mathbf{n}}(t -dt)(t-dt)]P[\dot{n}_{i}(t-dt)|\mathbf{V}(t-dt)],\]

where the dynamics of the membrane potential are deterministic given a particular history of the spike trains,

\[P[V_{i}(t)|\dot{\mathbf{n}}(t-dt)]\propto\delta\left(\frac{dV_{i}}{dt}+\tau_{ m}^{-1}(V_{i}-\varepsilon_{I})-I_{i}-\tau_{s}^{-1}\left(\mu_{\rm ext}-J_{\rm self }\dot{n}_{i}(t)+\sum_{j}w_{ij}\dot{n}_{j}(t)\right)\right).\]

Here, the proportionality hides a Jacobian factor that arises from a change of variables from \(V_{I}(t)\) to \(\dot{V}_{I}(t)\); this factor is constant for an Ito time discretization, which we assume here.

Next, we take the spike train process to be conditionally Poisson given the current value of the membrane potentials

\[P[\dot{n}_{i}(t-dt)|V_{i}(t-dt)]=\frac{\phi(V_{i}(t-dt))^{\dot{n}_{i}(t-dt)dt} }{(\dot{n}_{i}(t-dt)dt)!}e^{-\phi(V_{i}(t-dt))dt},\]

giving an overall representation

\[P[\mathbf{V}(t),\dot{\mathbf{n}}(t)]=\prod_{t,i} \delta\left(\frac{dV_{i}}{dt}+\tau_{m}^{-1}(V_{i}-\varepsilon_{I})-I_{i}- \tau_{s}^{-1}\left(\mu_{\rm ext}-J_{\rm self}\dot{n}_{i}(t)+\sum_{j}w_{ij} \dot{n}_{j}(t)\right)\right)\] \[\times\left[\frac{\phi(V_{i}(t-dt))^{\dot{n}_{i}(t-dt)dt}}{( \dot{n}_{i}(t-dt)dt)!}e^{-\phi(V_{i}(t-dt))dt}\right].\]

In order to cast this in a path integral representation, the standard approach is to represent the probability distributions in terms of a Fourier space representation. For the \(\delta\)-distribution we have

\[\delta(x)=\int_{-i\infty}^{i\infty}\frac{d\tilde{x}}{2\pi}e^{-\tilde{x}x},\]

and for a Poisson distribution with rate \(\lambda\) we have

\[p(x)=\int_{-i\infty}^{i\infty}\frac{d\tilde{x}}{2\pi}e^{-\tilde{x}x+W(\tilde{x })}=\int_{-i\infty}^{i\infty}\frac{d\tilde{x}}{2\pi}e^{-\tilde{x}x+\lambda( \varepsilon^{\tilde{x}}-1)},\]

where \(W(\tilde{x})=\lambda(\exp(\tilde{x})-1)\) is the cumulant generating function for the Poisson process. We have adopted the standard notation from physics of writing the auxiliary variables this process introduces with tildes, and absorbing the factor of the imaginary unit \(i\) into the notation (giving imaginary units of integration). The path integral representation of the spiking process above is then given by

\[P[\mathbf{V}(t),\dot{\mathbf{n}}(t)]=\int\mathfrak{D}[\tilde{\mathbf{V}}, \tilde{\mathbf{n}}]e^{-S[\tilde{\mathbf{V}},\mathbf{V},\tilde{\mathbf{n}}, \dot{\mathbf{n}}]},\]

where \(S[\tilde{\mathbf{V}},\mathbf{V},\tilde{\mathbf{n}},\dot{\mathbf{n}}]\) is referred to as the "action" of the process. We take the continuous-time limit, converting the product over time into an integral over time in the exponent. For this particular model, the action is given by

\[S[\tilde{\mathbf{V}},\mathbf{V},\tilde{\mathbf{n}},\dot{\mathbf{ n}}]= \int dt\ \sum_{i=1}^{n}\Bigg{\{}\tilde{V}_{i}\Bigg{[}\dot{V}_{i}+\frac{V_{ i}-\varepsilon_{i}}{\tau_{m}}-I_{i}-\tau_{s}^{-1}\left(\mu_{\rm ext}-J_{\rm self }\dot{n}_{i}(t)+\sum_{j}w_{ij}\dot{n}_{j}(t)\right)\Bigg{]}\] \[+\tilde{n}_{i}(t)\dot{n}_{i}(t)-\left(e^{\dot{n}_{i}(t)}-1\right) \phi(V_{i})\Bigg{\}}.\]
For our purposes, it will be convenient to marginalize out the dynamics of the spiking process \(\dot{\mathbf{n}}(t)\) and its conjugate variable \(\tilde{\mathbf{n}}(t)\) to obtain a representation for the stochastic dynamics of the membrane potentials (along with their auxiliary variables \(\tilde{\mathbf{V}}(t)\)). The spike-marginalized action is

\[S[\tilde{\mathbf{V}},\mathbf{V}]=\int dt\ \sum_{i=1}^{n}\left\{\tilde{V}_{i} \left[\dot{V}_{i}+\frac{V_{i}-\varepsilon_{I}}{\tau_{m}}-I_{i}-\tau_{s}^{-1} \mu_{\text{ext}}\right]-\left(e^{\tau_{s}^{-1}\left(-J_{\text{self}}\tilde{V}_ {i}+\sum_{j}\tilde{V}_{j}w_{ji}\right)}-1\right)\phi(V_{i})\right\}.\]

The Gaussian process approximation is derived by expanding this action around the mean-field solution, retaining only terms up to quadratic order in \(\mathbf{V}(t)-\mathbf{V}^{\text{mf}}\) and \(\tilde{\mathbf{V}}(t)\). The mean-field solution is obtained by the saddle-points of the action with respect to \(\mathbf{V}(t)\) and \(\tilde{\mathbf{V}}(t)\), which reproduce Eqn. (A2) for \(\mathbf{V}^{\text{mf}}\) and yield \(\tilde{\mathbf{V}}^{\text{mf}}=\mathbf{0}\). We thus perform a functional Taylor series expansion of the action around \((\tilde{\mathbf{V}},\mathbf{V})=(\mathbf{0},\mathbf{V}^{\text{mf}})\), keeping only terms to the second order in \(\delta\mathbf{V}=\mathbf{V}-\mathbf{V}^{\text{mf}}\) and \(\tilde{\mathbf{V}}\). The result is

\[S[\tilde{V},V] =\frac{1}{2}\int dtdt^{\prime}\ \sum_{ij}\tilde{V}_{i}(t)\left[- \tau_{s}^{-2}\sum_{k}\left(-\delta_{ik}J_{\text{self}}+w_{ik}\right)\left(- \delta_{jk}J_{\text{self}}+w_{jk}\right)\phi(V_{k}^{\text{mf}})\right]\tilde{ V}_{j}(t^{\prime})\] \[\quad+\int dtdt^{\prime}\ \sum_{ij}\tilde{V}_{i}(t)\left[\delta_{ij} \delta(t-t^{\prime})\frac{d}{dt}+\delta_{ij}\left(\tau_{m}^{-1}+\tau_{s}^{-1} J_{\text{self}}\phi^{\prime}(V_{j}^{\text{mf}})\right)-\tau_{s}^{-1}w_{ij}\phi^{ \prime}(V_{j}^{\text{mf}})\right]\delta V_{j}(t^{\prime}).\]

The form of the truncated action is the same as the path integral representation of an Ornstein-Uhlenbeck process derived explicitly by Chow and Buice [55]. We may therefore match terms to identify the effective stochastic process described by this action:

\[\frac{d\delta V_{i}}{dt}=-\sum_{j=1}^{n}\left[\delta_{ij}\left(\tau_{m}^{-1}+ \tau_{s}^{-1}J_{\text{self}}\phi^{\prime}(V_{j}^{\text{mf}})\right)-\tau_{s}^ {-1}w_{ij}\phi^{\prime}(V_{j}^{\text{mf}})\right]\delta V_{j}+\xi_{i}(t)\quad \text{for }i=1,2,...,n,\]

where \(\xi_{i}(t)\) is a zero-mean Gaussian noise with covariance

\[\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\tau_{s}^{-2}\sum_{k}\left(-\delta _{ik}J_{\text{self}}+w_{ik}\right)\left(-\delta_{jk}J_{\text{self}}+w_{jk} \right)\phi(V_{k}^{\text{mf}})\delta(t-t^{\prime}).\]

Casting this as a proper Ito stochastic differential equation, we get

\[d\delta\mathbf{V}=-\mathbf{A}\delta\mathbf{V}dt+\mathbf{\Sigma}d\mathbf{W}_{t}\]

or equivalently

\[d\mathbf{V}=\mathbf{A}\left(\mathbf{V}^{\text{mf}}-\mathbf{V}\right)dt+ \mathbf{\Sigma}d\mathbf{W}_{t},\]

where

\[A_{ij} =\delta_{ij}\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{ \prime}(V_{j}^{\text{mf}})\right)-\tau_{s}^{-1}w_{ij}\phi^{\prime}(V_{j}^{ \text{mf}})\] \[\left(\mathbf{\Sigma}\mathbf{\Sigma}^{T}\right)_{ij} =\tau_{s}^{-2}\sum_{k}\left(-\delta_{ik}J_{\text{self}}+w_{ik} \right)\left(-\delta_{jk}J_{\text{self}}+w_{jk}\right)\phi(V_{k}^{\text{mf}}).\]

In deriving the reduced dynamics for the population averages, we begin with the Langevin dynamics derived for the full network. We consider the network to have weakly heterogeneous populations in which the connections \(w_{ij}=w_{IJ}x_{ij}\) are given by Bernoulli variables, i.e. \(w_{ij}=w_{IJ}x_{ij}\) where \(w_{IJ}\) is a constant depending on the pre- and post-synaptic population identities (\(J\) and \(I\), respectively). We take each connection variable \(x_{ij}\) to be independent:

\[x_{ij}\sim\text{Bernoulli}(p_{IJ}).\]

We formally define the average of variable \(A_{i}\) across population \(I\) as

\[\langle\langle A_{i}\rangle\rangle_{I}=A_{I}\equiv\frac{1}{N_{I}}\sum_{i\in I}A_ {i}(t).\]
At this point, we write the population-averaged connection weights as follows:

\[\langle\langle w_{ij}\rangle\rangle_{I}\approx p_{IJ}w_{IJ}.\]

We will derive the effective equations for \(V_{I=0}\equiv V_{i=0}\) (the test neuron) and the population averages

\[V_{I=1} \equiv\frac{1}{N_{1}}\sum_{i\in 1}V_{i},\] \[V_{I=2} \equiv\frac{1}{N_{2}}\sum_{i\in 2}V_{i}.\]

We make mean-field-like approximations on the population-average of terms like \(\langle\langle f(A_{i})\rangle\rangle_{I}\approx f(\langle\langle A_{i}\rangle \rangle_{i}=f(A_{I})\), and we additionally assume approximate independence between the distributions of the synaptic connections, the stationary mean-field potentials \(V_{i}^{\text{mf}}\), and the potentials \(V_{i}\). We thus have

\[\frac{d}{dt}\left(\frac{1}{N_{I}}\sum_{i\in I}V_{i}\right) =\left\langle\left\langle\sum_{j}\left[\delta_{ij}\left(\tau_{m}^ {-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{\prime}(V_{j}^{\text{mf}})\right)- \tau_{s}^{-1}w_{ij}\phi^{\prime}(V_{j}^{\text{mf}})\right]\left(V_{j}^{\text{ mf}}-V_{j}\right)+\xi_{i}(t)\right\rangle\right\rangle_{I}\] \[=\left\langle\left\langle\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{ \text{self}}\phi^{\prime}(V_{i}^{\text{mf}})\right)V_{i}^{\text{mf}}\right\rangle _{I}-\left\langle\left\langle\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}} \phi^{\prime}(V_{i}^{\text{mf}})\right)V_{i}\right\rangle\right\rangle_{I}\right\rangle _{I}\] \[\approx\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{ \prime}\left((V_{i}^{\text{mf}})_{I}\right)\right)\langle(V_{i}^{\text{mf}})_{ I}-\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{\prime}\left((V_{i}^{ \text{mf}})_{I}\right)\right)\left\langle(V_{i}\right)\rangle_{I}\] \[\approx\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{ \prime}\left(V_{I}^{\text{mf}}\right)\right)\left(V_{I}^{\text{mf}}-V_{I} \right)-\tau_{s}^{-1}\left\langle\left\langle\sum_{j}N_{J}\langle\langle w_{ ij}\phi^{\prime}(V_{j}^{\text{mf}})_{J}\right\rangle\left\rangle\langle(V_{j}^{ \text{mf}})_{J}\right\rangle\right\rangle_{I}\] \[\qquad\qquad\qquad-\tau_{s}^{-1}\left\langle\left\langle\sum_{j} N_{J}\langle\langle w_{ij}\rangle\rangle_{J}\phi^{\prime}\left((V_{j}^{\text{mf}})_{J} \right)\langle(V_{j})_{J}\right\rangle\right\rangle\right\rangle_{I}+\Xi_{I}(t)\] \[\approx\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\text{self}}\phi^{ \prime}\left(V_{I}^{\text{mf}}\right)\right)\left(V_{I}^{\text{mf}}-V_{I} \right)-\tau_{s}^{-1}\sum_{J}N_{J}p_{IJ}w_{IJ}\phi^{\prime}\left(V_{J}^{\text{ mf}}\right)\left(V_{J}^{\text{mf}}-V_{J}\right)+\Xi_{I}(t)\] \[\Rightarrow\frac{dV_{I=1}}{dt} \approx\sum_{J}\left[\delta_{IJ}(\tau_{m}^{-1}+\tau_{s}^{-1}J_{ \text{self}}\phi^{\prime}(V_{I}^{\text{mf}}))-\tau_{s}^{-1}w_{IJ}p_{IJ}N_{J} \phi^{\prime}(V_{J}^{\text{mf}})\right](V_{J}^{\text{mf}}-V_{J})+\Xi_{I}(t).\]

In the last line above, the population-averaged effective noise processes are defined by \(\Xi_{I}(t)=\frac{1}{N_{I}}\sum_{i\in I}\xi_{i}(t)\), and the sum over \(J\) is over an arbitrary definition of sub-populations. In our particular case, we have \(J\in\{0,1,2\}\) as defined in Sec. II.4 with \(N_{0}=1\).

Next, we calculate the covariance of the population-averaged noise processes \(\Xi_{I}(t)\). We make the mean-field-like approxima
tions as before:

\[\langle\Xi_{I},\Xi_{J}\rangle= \Bigg{\langle}\frac{1}{N_{I}}\sum_{i\in I}\xi_{i},\frac{1}{N_{J}} \sum_{j\in J}\xi_{j}\Bigg{\rangle}\] \[=\frac{1}{N_{I}N_{J}}\sum_{i\in I,j\in J}\Big{[}\langle\xi_{i}\xi _{j}\rangle-\langle\xi_{i}\rangle\langle\xi_{j}\rangle\Big{]}\] \[\approx\frac{\tau_{s}^{-2}}{N_{I}N_{J}}\sum_{i\in I,j\in J} \Bigg{[}\sum_{k}\left(-\delta_{ik}J_{\rm self}+w_{ik}\right)\left(-\delta_{jk} J_{\rm self}+w_{jk}\right)\phi(V_{k}^{\rm mf}\right]\biggr{]}\delta(t-t^{ \prime})\] \[=\frac{\tau_{s}^{-2}}{N_{I}N_{J}}\sum_{i\in I,j\in J}\Big{[} \delta_{ij}J_{\rm self}^{2}\phi(V_{i}^{\rm mf})-J_{\rm self}w_{ji}\phi(V_{i}^ {\rm mf})-w_{ij}J_{\rm self}\phi(V_{j}^{\rm mf})+\sum_{K}\sum_{k\in K}w_{ik}w_ {jk}\phi(V_{k}^{\rm mf})\Big{]}\delta(t-t^{\prime})\] \[=\frac{\tau_{s}^{-2}}{N_{I}N_{J}}\Big{[}\delta_{IJ}N_{I}\langle \langle J_{\rm self}^{2}\phi(V_{i}^{\rm mf})\rangle\rangle_{I}-N_{I}N_{J} \Big{\langle}\Big{\langle}J_{\rm self}w_{ji}\phi(V_{i}^{\rm mf})\Big{\rangle} \Big{\rangle}_{I,J}-N_{I}N_{J}\Big{\langle}\Big{\langle}w_{ij}J_{\rm self}\phi (V_{j}^{\rm mf})\Big{\rangle}\Big{\rangle}_{I,J}\] \[\qquad\qquad+N_{I}N_{J}\Biggl{\langle}\Biggl{\langle}\sum_{K} \sum_{k\in K}\left\langle\left\langle w_{ik}w_{jk}\phi(V_{k}^{\rm mf})\right \rangle\right\rangle_{K}\Biggr{\rangle}\Biggr{\rangle}_{I,J}\biggr{]}\delta(t-t ^{\prime})\] \[=\tau_{s}^{-2}\Big{[}\delta_{IJ}\frac{J_{\rm self}^{2}}{N_{I}} \phi(V_{I}^{\rm mf})-J_{\rm self}\langle\langle w_{ji}\rangle\rangle_{I,J} \phi(V_{I}^{\rm mf})\] \[\qquad\qquad-\langle\langle w_{IJ}\rangle\rangle_{I,J}J_{\rm self }\phi(V_{J}^{\rm mf})+\sum_{K}\langle\langle w_{jk}\rangle\rangle_{K,I,J}N_{K }\phi(V_{K}^{\rm mf})\Big{]}\delta(t-t^{\prime})\] \[=\tau_{s}^{-2}\Bigl{[}\delta_{IJ}\frac{J_{\rm self}^{2}}{N_{I}} \phi(V_{I}^{\rm mf})-J_{\rm self}p_{JI}w_{JI}\phi(V_{I}^{\rm mf})-p_{JI}w_{IJ} J_{\rm self}\phi(V_{j}^{\rm mf})+\sum_{K}p_{JK}w_{JK}N_{K}\phi(V_{K}^{\rm mf}) \Bigr{]}\delta(t-t^{\prime})\] \[=\tau_{s}^{-2}\sum_{K}\left(-\delta_{IK}\frac{J_{\rm self}}{N_{K}} +p_{IK}w_{IK}\right)\left(-\delta_{JK}\frac{J_{\rm self}}{N_{K}}+p_{JK}w_{JK} \right)N_{K}\phi(V_{K}^{*})\delta(t-t^{\prime}).\]

Note that in this derivation we are assuming an equivalence between the temporal mean-field membrane potential for each individual neuron \(V_{i}\) (used in the previous section) with the mean-field value of the population-averaged membrane potential \(V_{I}\). This amounts to saying the network is sufficiently large and thus the mean of the membrane potential \(V_{i}\) for \(i\in I\) tends toward the mean of \(V_{I}\). This yields stochastic differential equation of the form

\[d{\bf V}={\bf A}\left({\bf V}^{\rm mf}-{\bf V}\right)dt+{\bf\Sigma}d{\bf W}_{t},\]

where

\[A_{ij}=\delta_{ij}\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\rm self}\phi^{\prime}( V_{j}^{\rm mf})\right)-\tau_{s}^{-1}w_{ij}\phi^{\prime}(V_{j}^{\rm mf})\]

\[\left({\bf\Sigma}{\bf\Sigma}^{T}\right)_{ij}=\tau_{s}^{-2}\sum_{k}\left(- \delta_{ik}J_{\rm self}+w_{ik}\right)\left(-\delta_{jk}J_{\rm self}+w_{jk} \right)\phi(V_{k}^{\rm mf}).\]

## Appendix B Gaussian process approximation of a population-averaged network

In this appendix, we derive the reduced model from Eqn. 6 by first averaging the Hawkes process dynamics across sub-populations and then making a Gaussian approximation, reversing the order of operations in Appendix A. We begin with the base model:

\[\frac{dV_{i}}{dt}=-\tau_{m}^{-1}(V_{i}-\varepsilon_{I})+I_{i}+\tau_{s}^{-1} \left(\mu_{\rm ext}-J_{\rm self}\dot{n}_{i}(t)+\sum_{J}\sum_{j\in J}w_{ij} \dot{n}_{j}(t)\right)\]

\[\dot{n}_{i}(t)dt\sim{\rm Poiss}[\phi(V_{i}(t))dt].\]

The population-averaged membrane potential dynamics are given by

\[\frac{d}{dt}V_{I}= \frac{d}{dt}\langle\langle V_{i}\rangle\rangle_{I}\] \[=-\frac{\langle\langle V_{i}\rangle\rangle_{I}-\varepsilon_{I}}{ \tau_{m}}+\langle\langle I_{i}\rangle\rangle_{I}+\frac{\mu_{\rm ext}}{\tau_{s}} -\tau_{s}^{-1}J_{\rm self}\langle\langle\dot{n}_{i}(t)\rangle\rangle_{I}+\tau_{s }^{-1}\Biggl{\langle}\Biggl{\langle}\sum_{J}\sum_{j\in J}w_{ij}\dot{n}_{j}(t) \Biggr{\rangle}\Biggr{\rangle}\] \[=-\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\frac{\mu_{\rm ext}} {\tau_{s}}-\tau_{s}^{-1}J_{\rm self}\langle\langle\dot{n}_{i}(t)\rangle\rangle_{I}+ \tau_{s}^{-1}\sum_{J}\sum_{j\in J}\langle\langle w_{ij}\rangle\rangle\dot{n}_{j}(t).\]
As before, we take the connections \(w_{ij}\) to be scaled Bernoulli variables, i.e. \(w_{ij}=w_{IJ}x_{ij}\) where \(w_{IJ}\) is a constant depending on the pre- and post-synaptic population identities (\(J\) and \(I\), respectively) and \(x_{ij}\sim\text{Bernoulli}(p_{IJ})\). The population-averaged connections are again given by \(\langle\langle w_{ij}\rangle\rangle_{I}\approx p_{IJ}w_{IJ}\). We next re-cast the spiking processes into population spiking processes using the following definition

\[\dot{m}_{I}(t)\equiv\sum_{i\in I}\dot{n}_{i}(t)=N_{I}\langle\langle\dot{n}_{i} \rangle\rangle_{I}.\]

As each \(\dot{m}_{I}(t)\) is a sum of conditionally-Poisson processes, it is also a conditionally-Poisson process. Using the same mean-field-esque approximation as before, we may approximate the conditional rate of each \(\dot{m}_{I}(t)\) as follows:

\[\dot{m}_{I}=\sum_{i\in I}\dot{n}_{i}(t) \sim\text{Poiss}\left(\sum_{i\in I}\phi(V_{i}(t))dt\right)=\text{ Poiss}\left(N_{I}\langle\langle\phi(V_{i}(t))\rangle\rangle_{I}dt\right)\] \[\approx\text{Poiss}\left(N_{I}\phi(\langle\langle V_{i}(t)\rangle \rangle_{I}dt)\right.\] \[=\text{Poiss}\left(N_{I}\phi(V_{I}(t))dt\right).\]

With this, the population-averaged Hawkes process dynamics become

\[\frac{d}{dt}V_{I}= -\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\frac{\mu_{\text{ ext}}}{\tau_{s}}-\tau_{s}^{-1}J_{\text{self}}\frac{\dot{m}_{I}(t)}{N_{I}}+\tau_{s}^{- 1}\sum_{J}p_{IJ}w_{IJ}\dot{m}_{J}(t)\] \[=-\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\tau_{s}^{-1}\left( \mu_{\text{ext}}+\sum_{J}\left(-\delta_{IJ}\frac{J_{\text{self}}}{N_{I}}+p_{ IJ}w_{IJ}\right)\dot{m}_{J}(t)\right)\] \[\dot{m}_{I}(t)dt\sim\text{Poiss}[N_{I}\phi(V_{I}(t))dt].\]

After deriving the population-averaged dynamics for the nonlinear Hawkes process, we apply the Gaussian-process approximation scheme to the new dynamics. We begin by applying a mean-field-like approximation to the average of the population-spiking processes, namely \(\langle\dot{m}_{I}(t)\rangle\approx N_{I}\phi(V_{I})\). This is used to find the stationary mean-field solution for the population-averaged membrane potential dynamics, given by a set of transcendental equations

\[V_{I}^{\text{mf}}=\varepsilon_{I}+\tau_{m}I_{I}+\frac{\tau_{m}}{\tau_{s}} \left(\mu_{\text{ext}}+\sum_{J}\left(-\delta_{IJ}\frac{J_{\text{self}}}{N_{I} }+p_{IJ}w_{IJ}\right)N_{J}\phi(V_{J}^{\text{mf}})\right).\]

As in Appendix A, we represent the joint probability distribution \(P[\mathbf{V}(t),\dot{\mathbf{m}}(t)]\) as a path integral by discretizing time, making appropriate Fourier transforms, and taking a continuous-time limit. This yields the expression

\[P[\mathbf{V}(t),\dot{\mathbf{m}}(t)]=\int\mathfrak{D}[\tilde{\mathbf{V}}, \tilde{\mathbf{m}}]e^{-S[\tilde{\mathbf{V}},\mathbf{V},\dot{\mathbf{m}},\dot{ \mathbf{m}}]},\]

with

\[S[\tilde{\mathbf{V}},\mathbf{V},\tilde{\mathbf{m}},\dot{\mathbf{ m}}]= \int dt\;\sum_{I=0,1,2}\Bigg{\{}\tilde{V}_{I}\left[\dot{V}_{I}+\frac{V_{I}- \varepsilon_{i}}{\tau_{m}}-I_{I}-\tau_{s}^{-1}\left(\mu_{\text{ext}}+\sum_{J} \left(-\delta_{IJ}\frac{J_{\text{self}}}{N_{I}}+p_{IJ}w_{IJ}\right)\dot{m}_{J} (t)\right)\right]\] \[+\tilde{m}_{I}(t)\dot{m}_{I}(t)-\left(e^{\dot{m}_{I}(t)}-1\right) N_{I}\phi(V_{I})\Bigg{\}}.\]

We marginalize out the explicit spiking dynamics as before by finding the zeros of the derivatives of the action w.r.t. \(\dot{\mathbf{m}}(t)\) and its conjugate variables \(\tilde{\mathbf{m}}(t)\). This yields the following marginalized action:

\[S[\tilde{\mathbf{V}},\mathbf{V}]=\] \[\int dt\;\sum_{I=0,1,2}\left\{\tilde{V}_{I}\left[\dot{V}_{I}+ \frac{V_{I}-\varepsilon_{I}}{\tau_{m}}-I_{I}-\tau_{s}^{-1}\mu_{\text{ext}} \right]-\left(e^{\tau_{s}^{-1}\left(-\frac{J_{\text{self}}}{N_{I}}\tilde{V}_{I} +\sum_{J}\tilde{V}_{I}p_{JI}w_{JI}\right)}-1\right)N_{I}\phi(V_{I})\right\}.\]

We expand this action around the mean-field solution \((\tilde{\mathbf{V}},\mathbf{V})=(\mathbf{0},\mathbf{V}^{\text{mf}})\) to quadratic order. Evaluating individual terms and derivatives at the mean-field solution, we get

\[S[\mathbf{0},\mathbf{V}^{\text{mf}}]=0,\] \[S_{V_{I}}[\mathbf{0},\mathbf{V}^{\text{mf}}]=0,\] \[S_{\tilde{V}_{I}}[\mathbf{0},\mathbf{V}^{\text{mf}}]=\int dt\; \left[\dot{V}_{i}+\frac{V_{i}^{\text{mf}}-\varepsilon_{I}}{\tau_{m}}-\frac{V_ {i}^{\text{mf}}-\varepsilon_{I}}{\tau_{m}}\right]=\int dt\;\left[\dot{V}_{i}- \dot{V}_{i}^{\text{mf}}\right]=\int dt\;\delta V_{i},\] \[S_{\tilde{V}_{I}\tilde{V}_{J}}[\mathbf{0},\mathbf{V}^{\text{mf}}] =\int dt\;\Bigg{[}-\tau_{s}^{-2}\sum_{K}\left(-\delta_{IK}\frac{J_{\text{self} }}{N_{I}}+p_{IK}w_{IK}\right)\left(-\delta_{JK}\frac{J_{\text{self}}}{N_{J}}+ p_{JK}w_{JK}\right)N_{K}\phi(V_{K})\Bigg{]},\]
\[S_{\tilde{V}_{I}V_{J}}[\mathbf{0},\mathbf{V}^{\mathrm{mf}}]=\int dt\ \left[\delta_{IJ}\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\mathrm{self}}\phi^{ \prime}(V_{I}^{\mathrm{mf}})\right)-\tau_{s}^{-1}p_{IJ}w_{IJ}N_{J}\phi^{\prime} (V_{J}^{\mathrm{mf}})\right],\]

and

\[S_{V_{I}V_{J}}[\mathbf{0},\mathbf{V}^{\mathrm{mf}}]=0.\]

Again defining fluctuations in the membrane potential as \(\delta V_{I}:=V_{I}-V_{I}^{\mathrm{mf}}\), approximated action can be written as

\[S[\tilde{V},V] =\int dt\ \sum_{I}\left\{\delta V_{I}+\sum_{J}\left[\delta_{IJ} \left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\mathrm{self}}\phi^{\prime}(V_{I}^{ \mathrm{mf}})\right)-\tau_{s}^{-1}p_{IJ}w_{IJ}N_{J}\phi^{\prime}(V_{J}^{ \mathrm{mf}})\right]\delta V_{J}\right\}\tilde{V}_{I}(t)\] \[\quad+\frac{1}{2}\int dtdt^{\prime}\ \sum_{ij}\left[-\tau_{s}^{-2} \sum_{K}\left(-\delta_{IK}\frac{J_{\mathrm{self}}}{N_{I}}+p_{IK}w_{IK}\right) \left(-\delta_{JK}\frac{J_{\mathrm{self}}}{N_{J}}+p_{JK}w_{JK}\right)N_{K} \phi(V_{K})\right]\tilde{V}_{I}(t)\tilde{V}_{J}(t^{\prime}).\]

As before, we can identify the GPA dynamics of the population-averaged Hawkes process as corresponding to an Ornstein-Uhlenbeck process. We may therefore match terms to identify the effective stochastic process described by this action:

\[\frac{d\delta V_{I}}{dt}=-\sum_{J=0,1,2}\left[\delta_{IJ}\left(\tau_{m}^{-1}+ \tau_{s}^{-1}J_{\mathrm{self}}\phi^{\prime}(V_{I}^{\mathrm{mf}})\right)-\tau_ {s}^{-1}p_{IJ}w_{IJ}N_{J}\phi^{\prime}(V_{J}^{\mathrm{mf}})\right]\delta V_{J} +\xi_{I}(t)\quad\text{for }I=0,1,2\]

where \(\xi_{I}(t)\) is a zero-mean Gaussian noise with covariance

\[\langle\xi_{I}(t)\xi_{J}(t^{\prime})\rangle=\tau_{s}^{-2}\sum_{K}\left(-\delta _{IK}\frac{J_{\mathrm{self}}}{N_{I}}+p_{IK}w_{IK}\right)\left(-\delta_{JK} \frac{J_{\mathrm{self}}}{N_{J}}+p_{JK}w_{JK}\right)N_{K}\phi(V_{K})\delta(t-t^ {\prime}).\]

Casting this as a proper Ito stochastic differential equation, we get

\[d\mathbf{V}=\mathbf{A}\left(\mathbf{V}^{\mathrm{mf}}-\mathbf{V}\right)dt+ \mathbf{\Sigma}d\mathbf{W}_{t},\]

where

\[A_{ij}=\delta_{IJ}\left(\tau_{m}^{-1}+\tau_{s}^{-1}J_{\mathrm{self}}\phi^{ \prime}(V_{I}^{\mathrm{mf}})\right)-\tau_{s}^{-1}p_{IJ}w_{IJ}N_{J}\phi^{ \prime}(V_{J}^{\mathrm{mf}})\]

We note that this is consistent with the form derived in Appendix A.

## Appendix C Population averaging for the linear non-spiking model

We also construct a simpler model of networked, linear non-spiking (or "graded potential") neurons. We assume the neurons are injected with large numbers of synaptic input that sum together to be approximately Gaussian, with non-zero mean \(\mu_{\mathrm{ext}}\), creating a stochastic system with dynamics described by

\[\frac{dV_{i}}{dt}= -\tau_{m}^{-1}(V_{i}-\varepsilon_{I})+I_{i}+\tau_{s}^{-1}\left( \mu_{\mathrm{ext}}-J_{\mathrm{self}}\phi(V_{i})+\sum_{j}w_{ij}\phi(V_{j}) \right)+\xi_{i}(t).\]

We begin this derivation by assuming the connections \(w_{ij}=w_{IJ}x_{ij}\) are scaled Bernoulli variables as in Appendices A,B. Here, the transfer function \(\phi(\cdot)\) is a simple linear function (i.e. \(\phi(x)=x\)). The processes \(\xi_{i}(t)\) are zero-mean Gaussian noise synaptic input from neurons external to the network being examined, and thus they scale with \(\tau_{s}^{-1}\) (i.e. \(\xi_{i}(t)\sim\tau_{s}^{-1}\)). We define the covariance of the noise processes \(\xi_{i}(t)\) as follows:

\[\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\tau_{s}^{-2}\mu_{\mathrm{ext}} \delta(t-t^{\prime}).\]

Here, \(k_{J}\) is a constant potentially depending on the identity of the receiving population \(J\). We wish to derive a population-averaged model for the membrane potential dynamics for comparison to the Gaussian-process-approximated spiking models. Again, we define

\[\langle\langle A_{i}\rangle\rangle_{I}\equiv\frac{1}{N_{I}}\sum_{i\in I}A_{i}( t).\]

We thus derive the population-averaged dynamics for population \(I\):

\[\frac{d}{dt}\left(\frac{1}{N_{I}}\sum_{i\in I}V_{i}\right) =-\left\langle\left\langle\frac{V_{i}-\varepsilon_{I}}{\tau_{m}} \right\rangle\right\rangle_{I}+\langle\langle I_{i}\rangle\rangle_{I}+\tau_ {s}^{-1}\left(\mu_{\mathrm{ext}}-\langle\langle J_{\mathrm{self}}\phi(V_{i}) \rangle\rangle_{I}+\left\langle\left\langle\sum_{J}N_{J}\langle\langle w_{ij} \phi(V_{j})\rangle\rangle_{J}\right\rangle\right)_{I}\right)+\langle\langle \langle\xi_{i}(t)\rangle\rangle_{I}\] \[\Rightarrow\frac{dV_{I}}{dt}\approx -\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\tau_{s}^{-1}\left( \mu_{\mathrm{ext}}-J_{\mathrm{self}}\phi(V_{I})+\sum_{J}p_{IJ}w_{IJ}N_{J}\phi( V_{J})\right)+\Xi_{I}(t),\]
where we have defined \(\Xi_{I}(t)\equiv\frac{1}{N_{I}}\sum_{i\in I}\xi_{i}(t)\) for \(I=0,1,2\). The means and covariances of the population-averaged noise processes are as follows:

\[\langle\Xi_{I}(t)\rangle=\left\langle\frac{1}{N_{I}}\sum_{i\in I}\xi_{i}(t) \right\rangle=\frac{1}{N_{I}}\sum_{i\in I}\langle\xi_{i}(t)\rangle=0,\]

and

\[\langle\Xi_{I}(t),\Xi_{J}(t)\rangle= \Bigg{\langle}\frac{1}{N_{I}}\sum_{i\in I}\xi_{i}(t),\frac{1}{N_{ J}}\sum_{j\in J}\xi_{j}(t)\Bigg{\rangle}=\frac{1}{N_{I}N_{J}}\sum_{i\in I,j\in J} \langle\xi_{i}(t)\xi_{j}(t)\rangle-\langle\xi_{i}(t)\rangle\langle\xi_{j}(t)\rangle\] \[=\frac{1}{\tau_{s}^{2}N_{I}N_{J}}\sum_{i\in I,j\in J}\delta_{ij} \mu_{\text{ext}}\delta(t-t^{\prime})\] \[=\frac{\delta_{IJ}}{\tau_{s}^{2}N_{I}^{2}}\sum_{i\in I}\mu_{\text {ext}}\delta(t-t^{\prime})=\frac{\delta_{IJ}}{\tau_{s}^{2}N_{I}^{2}}N_{I}\mu_ {\text{ext}}\delta(t-t^{\prime})\] \[=\frac{\delta_{IJ}}{\tau_{s}^{2}N_{I}}\mu_{\text{ext}}\delta(t-t ^{\prime}).\]

We can then rewrite the population dynamics as

\[dV_{I}= \left(-\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\tau_{s}^{-1} \left(\mu_{\text{ext}}-J_{\text{self}}\phi(V_{I})+\sum_{J}p_{IJ}w_{IJ}N_{J} \phi(V_{J})\right)\right)dt+\Xi_{I}(t)dt\] \[\to d\mathbf{V}=\mathbf{A}\left(\mathbf{A}^{-1}\left(\tau_{s}^{-1} \mu_{\text{ext}}+I_{I}\right)-\mathbf{V}\right)dt+\mathbf{\Sigma}d\mathbf{W}_{t}\] \[= \mathbf{A}\left(\mu-\mathbf{V}\right)dt+\mathbf{\Sigma}d\mathbf{W }_{t},\]

where

\[\mathbf{A}_{IJ}= \delta_{IJ}\left(\tau_{m}^{-1}-\tau_{s}^{-1}J_{\text{self}} \right)+\tau_{s}^{-1}p_{IJ}w_{IJ}N_{J}\] \[=\delta_{IJ}\tau_{m}^{-1}+\tau_{s}^{-1}w_{IJ}^{*},\] \[w_{IJ}^{*}= -\delta_{IJ}J_{\text{self}}+p_{IJ}w_{IJ}N_{J},\] \[\left(\Sigma\Sigma^{T}\right)_{IJ}= \frac{\delta_{IJ}}{\tau_{s}^{2}N_{I}}\mu_{\text{ext}}.\]

## Appendix D Balance equations

To derive the balanced state conditions for the network, we begin with the population-averaged spiking network as derived in Appendix B:

\[\frac{d}{dt}V_{I}=-\frac{V_{I}-\varepsilon_{I}}{\tau_{m}}+I_{I}+\tau_{s}^{-1} \left(\mu_{\text{ext}}+\sum_{J=0,1,2}\left(-\delta_{IJ}\frac{J_{\text{self}} }{N_{I}}+p_{IJ}w_{IJ}\right)\dot{m}_{J}(t)\right)\]

\[\dot{m}_{I}(t)dt\sim\text{Poiss}[N_{I}\phi(V_{I}(t))dt],\]

where \(p_{IJ}w_{IJ}\) came from the population-averaged synaptic connection \(\langle\langle w_{ij}\rangle\rangle_{J}\) and the effective spike count processes are \(\dot{m}_{J}(t)=\sum_{j\in J}\dot{n}_{j}(t)\). The total external input to "neuron" \(I\) is \(I_{I}+\tau_{s}^{-1}\left(\mu_{\text{ext}}-\frac{J_{\text{self}}}{N_{I}}\dot{m} _{I}(t)+\sum_{J}p_{IJ}w_{IJ}\dot{m}_{J}(t)\right)\). We want to estimate the mean and variance of this input, taken over the stochastic process. The mean is straightforward, yielding

\[\tau_{s}^{-1}\kappa_{I}\equiv I_{I}+\tau_{s}^{-1}\left(\mu_{\text{ext}}-J_{ \text{self}}\phi(V_{I})+\sum_{J}p_{IJ}w_{IJ}N_{J}\phi(V_{J})\right).\]

Note that the correction term \(-J_{\text{self}}\dot{m}_{I}(t)\) is always going to be smaller than the \(\sum_{J}p_{IJ}w_{IJ}\dot{m}_{J}\) term, so for the purposes of the balanced condition calculation we will neglect it. For the current work, we take the injected currents \(I_{I}\) to be constants.

Calculating the covariance of the total input at times \(t\) and \(t^{\prime}\) yields

\[\sum_{JK}p_{IJ}w_{IJ}p_{IK}w_{IK}\Big{[}\langle\dot{m}_{J}(t)\dot{m}_{K}(t^{ \prime})\rangle-\langle\dot{m}_{J}(t)\rangle\langle\dot{m}_{K}(t^{\prime}) \rangle\Big{]}.\]

We make a Poisson approximation to replace the covariance of the \(\dot{m}\)'s with \(\langle\dot{m}_{J}(t)\rangle\delta_{JK}\delta(t-t^{\prime})\). Hence, the covariance becomes

\[\sum_{J}(p_{IJ}w_{IJ})^{2}N_{J}\phi(V_{J})\delta(t-t^{\prime}).\]
We want the variance of the synaptic input to be \(\mathcal{O}(N^{0})\), which means that to leading order we want

\[\sum_{J}(p_{IJ}w_{IJ})^{2}N_{J}\phi(V_{J})\approx(p_{I1}w_{I1})^{2}N_{1}\phi(V_{ 1})+(p_{I2}w_{I2})^{2}N_{2}\phi(V_{2})\sim\mathcal{O}(N^{0}).\]

We neglect the contribution from the test neuron because it is sub-leading here, i.e. \(N_{0}=1\ll N_{1},\,N_{2}\). In order for this expression to be order 1, we see that we need \(w_{IJ}\) to scale like \(1/\sqrt{N}\) as implemented in Eqns. 12 & 13.

We return to the mean input to neuron \(I\), which we will write as

\[\tau_{s}^{-1}\kappa_{I}\approx\sqrt{N}\left(\frac{I_{I}+\tau_{s}^{-1}\mu_{ \text{ext}}}{\sqrt{N}}+\tau_{s}^{-1}\left\{p_{I1}w_{I1}\frac{N_{1}}{\sqrt{N}} \phi(V_{1})+p_{I2}w_{I2}\frac{N_{2}}{\sqrt{N}}\phi(V_{2})\right\}\right).\]

For a balanced network \(\kappa_{I}\) should be \(\mathcal{O}(1)\) for all \(I\), which means that the terms in brackets must vanish faster than \(1/\sqrt{N}\). We assume that \(I_{I},\mu_{\text{ext}}\propto\sqrt{N}\), and because \(w_{IJ}\sim 1/\sqrt{N}\) and \(N_{I}\propto N\) (for \(I\neq 0\)), the terms in brackets are \(\mathcal{O}(1)\).

As \(N\to\infty\), the terms in brackets must vanish in order for \(\kappa_{I}\) to be finite. This yields a linear system of equations that uniquely determines the means \(\mu_{I}\equiv\phi(V_{I})\), and allows us to place constraints on the parameters:

\[-\begin{bmatrix}I_{1}+\tau_{s}^{-1}\mu_{\text{ext}}\\ I_{2}+\tau_{s}^{-1}\mu_{\text{ext}}\end{bmatrix}=\frac{1}{\tau_{s}}\begin{bmatrix} p_{11}w_{11}N_{1}&p_{12}w_{12}N_{2}\\ p_{21}w_{21}N_{1}&p_{22}w_{22}N_{2}\end{bmatrix}\begin{bmatrix}\phi(V_{1})\\ \phi(V_{2})\end{bmatrix}.\]

Solving this system of equations for the spike rates \(\phi(V_{I}^{\text{mf}})\), we get

\[\phi(V_{1}) =\frac{\tau_{s}}{N_{1}}\frac{p_{12}w_{12}\left(I_{2}+\tau_{s}^{- 1}\mu_{\text{ext}}\right)-p_{22}w_{22}\left(I_{1}+\tau_{s}^{-1}\mu_{\text{ext }}\right)}{p_{11}p_{22}w_{11}w_{22}-p_{12}p_{12}w_{21}w_{12}},\] \[\phi(V_{2}) =\frac{\tau_{s}}{N_{2}}\frac{p_{21}w_{21}\left(I_{1}+\tau_{s}^{- 1}\mu_{\text{ext}}\right)-p_{11}w_{11}\left(I_{2}+\tau_{s}^{-1}\mu_{\text{ext }}\right)}{p_{11}p_{22}w_{11}w_{22}-p_{12}p_{12}w_{21}w_{12}}.\]

In the case of our particular models, we can further reduce this expression by noting that \(I_{I}=0\) for \(I=1,2\) and \(p_{IJ}=p\,\,\forall\,\,I,J\):

\[\phi(V_{1}) =\frac{1}{pN_{1}}\frac{w_{12}-w_{22}}{w_{11}w_{22}-w_{21}w_{12}} \mu_{\text{ext}},\] \[\phi(V_{2}) =\frac{1}{pN_{2}}\frac{w_{21}-w_{11}}{w_{11}w_{22}-w_{21}w_{12}} \mu_{\text{ext}}.\]

We highlight here that \(\phi(V_{I})>0\) by its definition as a firing rate. Additionally, \(\mu_{\text{ext}}\) is assumed to by synaptic input projected into the local network and is thus positive (i.e. excitatory) here. Taken together, these two points mean the synaptic parameters must satisfy one of the two following sets of inequalities to be in a balanced regime:

\[\begin{cases}w_{11}w_{22}>w_{12}w_{21}\\ w_{12}>w_{22}\\ w_{21}>w_{11}\end{cases} \tag{12}\]

or

\[\begin{cases}w_{11}w_{22}<w_{12}w_{21}\\ w_{12}<w_{22}\\ w_{21}<w_{11}\end{cases}. \tag{13}\]

With this, we have derived the appropriate scaling for the various parameters in the model and found constraints for the synaptic strengths in order satisfy the necessary properties of a balanced network.
