Model-agnostic neural mean field with a data-driven transfer function

Abstract As one of the most complex systems known to science, modeling brain behavior and function is both fascinating and extremely difficult. Empirical data is increasingly available from ex vivo human brain organoids and surgical samples, as well as in vivo animal models, so the problem of modeling the behavior of large-scale neuronal systems is more relevant than ever. The statistical physics concept of a mean-field model offers a tractable way to bridge the gap between single-neuron and population-level descriptions of neuronal activity, by modeling the behavior of a single representative neuron and extending this to the population. However, existing neural mean-field methods typically either take the limit of small interaction sizes, or are applicable only to the specific neuron models for which they were derived. This paper derives a mean-field model by fitting a transfer function called Refractory SoftPlus, which is simple yet applicable to a broad variety of neuron types. The transfer function is fitted numerically to simulated spike time data, and is entirely agnostic to the underlying neuronal dynamics. The resulting mean-field model predicts the response of a network of randomly connected neurons to a time-varying external stimulus with a high degree of accuracy. Furthermore, it enables an accurate approximate bifurcation analysis as a function of the level of recurrent input. This model does not assume large presynaptic rates or small postsynaptic potential size, allowing mean-field models to be developed even for populations with large interaction terms.


Introduction
The brain is one of the most complex systems known to science, which makes the problem of computationally modeling its behavior and function both fascinating and extremely difficult.The computational substrate of the brain consists of billions of neurons, each with tens of thousands of connections spanning neurons from the local region to distant areas [1].Computational neuroscientists are thus faced with a problem of scale: while the biophysical behavior of the individual neuron is understood from first principles [2,3], and large numbers of neurons have been modeled phenomenologically as a dynamical system [4,5], unifying different scales of description remains a fundamental difficulty and the subject of ongoing research [6,7].Today, as neuromorphic systems are emerging for computation and robotic control [8][9][10][11][12][13] and human brain organoids and surgical samples are becoming key systems for studying disease states and neuronal computation [14][15][16][17][18][19], the question of modeling the behavior of large-scale neuronal systems is becoming all the more crucial.
One of the most popular approaches to unifying model scales is to make the mean-field assumption, namely that all the neurons in the population under consideration, as well as their inputs, are independent and statistically identical [20,21].This assumption is applicable to a variety of systems, ranging from randomly connected sparse networks [20] to modular networks of biological interest [22].When Figure 1.Schematic illustration of a single neuron in a typical network, where the total number of incident presynaptic inputs is large but firing rates and postsynaptic potential sizes are not.The firing rate of the neuron is given by a transfer function F of the total input rate R = Nr.
connections are not strongly selective, neurons within a population are essentially equivalent and can be represented by one population firing rate.Even networks with all-to-all connectivity can be treated, provided the neurons are stochastic enough to remain approximately independent [23].Under the mean-field assumption, the members of a neuronal population act as independent samples of a shared statistical distribution; this reduces many problems to the study of the firing rate of a single representative neuron receiving N inputs at an average rate r as depicted in figure 1.Similar approaches have been used extensively in statistical physics, for example in the study of coupled oscillators [24,25].
Analytical mean-field models can also be derived for some neuronal systems, such as networks of rate neurons [26] or of stochastic spiking neurons characterized by a renewal process [27].For dynamical spiking neurons, on the other hand, no such analytical approach is known.Thus current mean-field models examine the limit of an infinite number of presynaptic neurons with infinitesimal postsynaptic potentials, and apply the central limit theorem to replace these inputs with a single white Gaussian noise process.This so-called diffusion approximation transforms the single-neuron dynamics into a drift-diffusion problem, for which many techniques have been developed [28,29].Perturbation analysis can model complex interactions, but is only valid for small effects [30,31].For more general study of effects such as spike frequency adaptation (SFA) or synaptic dynamics, the Fokker-Planck equation can be solved to get the time evolution of the probability distribution of the neuron's membrane voltage under the drift-diffusion process [23,32].However, each solution is particular to a single neuronal model, and calculating solutions is computationally expensive even when mathematical analysis has revealed an explicit model [33].
One popular basis for constructing mean-field models is a neuron's 'transfer function' , which represents its firing rate as a function of the amount of input noise, analogous to the firing rate curve used to study the deterministic dynamics of an individual neuron [20,21].Analytical transfer functions are known for a variety of different neuron types and noise models [34][35][36].However, we are not aware of any mean-field model making use of a more general noise assumption than the simple diffusion approximation.
Mean-field models based on the diffusion approximation and the approximate transfer function of [20] have been applied to a variety of different model neurons [37], have been studied in terms of bifurcation theory [38,39], and have been compared to biological data [22].However, the sigmoidal shape of this transfer function does not match analytical results [36,40].It also contains an arbitrary timescale T, which must be large enough for a quasistatic assumption to hold, but smaller than the minimum interspike interval (ISI) in the system [41].Furthermore, it is explicitly assumed through moment closure that the membrane potential itself is normally distributed as well [42], which is inaccurate for realistic noise sources [43].
In this paper, we propose a novel, explicitly phenomenological approach based on numerically fitting a simple parametrized transfer function that can match both the shape of analytical and observed transfer functions and simulation results.It is also inexpensive to evaluate, as it has only four parameters, which are fitted with a good level of accuracy based on only a few simulations of single uncoupled neurons.Being based on a transfer function, this model is still limited by the assumption of a stationary firing rate distribution, but it makes no assumption about the form or dynamics of the membrane potential distribution.In particular, we do not require the diffusion approximation and can describe populations of neurons with finite numbers of connections and finitely large postsynaptic potentials.We numerically solve a consistency condition to find stable firing rates of a neuronal population, achieving good accuracy in a wide range of simulation examples.We demonstrate the effectiveness of our model by predicting a saddle-node bifurcation to bistability as a function of the connectivity parameter.Furthermore, we implement a popular first-order dynamical model and demonstrate its predictive power alongside its limitations.

Method
There exists a wide variety of dynamical spiking neural models.All are excitable, in the sense that increasing membrane voltage past a threshold causes a rapid upstroke followed by a return to equilibrium, called a spike [44].However, their underlying dynamics vary significantly, and this variation can have important computational implications [45].This paper includes results for three popular neuron models: the leaky integrate-and-fire (LIF), Izhikevich (Izh), and Hodgkin-Huxley (HH) models.

Leaky integrate-and-fire model
The dynamics of the LIF neuron are given in equation (1).The membrane potential of this neuron undergoes a drift back towards its resting state, while being driven by a time-dependent input ξ(t).It is also common to introduce a refractory period T ref , a duration for which the state of the neuron does not change after every spike. ( Synaptic input to the neuron is represented by a random process ξ(t).It is common in models to assume that each presynaptic neuron fires as a Poisson point process [46].However, it is well known that individual neurons frequently deviate from Poisson statistics [47][48][49][50].Fortunately, the superposition of N independent renewal processes with total rate R tends in the large-N limit to a Poisson process, without requiring that the input renewal processes be Poisson themselves, or even identically distributed [51].Therefore, our model makes the somewhat more general assumption that each presynaptic neuron fires as an independent renewal process.
Each presynaptic firing event is modeled as a Dirac delta which, when integrated, instantaneously increases the membrane voltage by its fixed postsynaptic potential (PSP) value q.A typical order of magnitude for these parameters in the neocortex is 10 4 neurons [1] firing at a rate of several hertz [20], with PSP magnitudes of about one millivolt [52], as in figure 1.

Izhikevich and Hodgkin-Huxley models
Besides the LIF neuron, we also consider the Izhikevich and Hodgkin-Huxley models.The Izhikevich model is a quadratic integrate-and-fire neuron augmented with an adaptation variable u [53].The result is a model not substantially more computationally expensive than LIF, but with much more realistic spike shapes.On the other hand, the HH model is based on biophysically modeling the dynamics of the membrane potential [2].This results in a more realistic but much more complex model.See appendix D for dynamical equations as well as a more thorough comparison of the dynamical features of all three models.
The Izhikevich model uses Dirac delta synapses.However, delta synapses are not supported by the HH neuron available in the simulator.Instead, alpha postsynaptic currents are used, where the weight w of a synapse gives the peak value of the postsynaptic current α(t) = wte 1−t/τsyn /τ syn , which is integrated by a membrane capacitance C m .To keep inputs comparable between neuron models, we compute w as a function of an equivalent PSP size q by fixing its integral so that it does not depend on the synaptic time constant and the zero-τ syn limit corresponds to a delta PSP of size q:

Excitatory-inhibitory balance
Although it is convenient to assume that the inputs to the neuron are as homogeneous as possible, models generally assume a distinction between excitatory and inhibitory inputs, which differ at least in the sign of their PSPs.In our simulations, a fraction η of the input connections are from excitatory neurons with PSP amplitude q e > 0, and the remainder of the input is inhibitory, with a negative PSP with absolute value q i .This applies both to the background input, where the input to each neuron is drawn from two independent Poisson sources with rates ηR and (1 − η)R corresponding to its own excitatory and inhibitory background, as well as to simulations with recurrent connectivity, where each neuron is either excitatory or inhibitory and produces the corresponding PSP in all of its postsynaptic neurons.
An exact transfer function is still unknown for even the simple LIF neuron under Poisson inputs [54].However, it has long been known that in the diffusion limit of infinitely large total input rate R and small PSPs, the membrane voltage of the LIF neuron follows an Ornstein-Uhlenbeck process with time constant equal to the membrane time constant τ [34].For a given excitatory fraction η, the drift rate µ and the diffusion constant D of this Ornstein-Uhlenbeck process would be given by: ( We are interested in modeling neurons in the asynchronous irregular state observed in biological systems, which computational experiments have shown occurs when spiking is driven by fluctuations rather than drift [39,47,55,56].The condition where mean synaptic input µ = 0, is called loose excitatory-inhibitory (EI) balance and is widely observed in experiments and commonly assumed in models [41,57].Setting µ = 0 fixes the inhibitory PSP size q i as a function of q e : In the diffusion limit, any inputs that yield the same value of D will behave identically, regardless of the value of η.We do not expect this to hold exactly because we are using finite values of R and q, but it remains a convenient parametrization.Therefore, we use equation ( 4) and (3) to define an effective diffusion coefficient for the balanced input regime as follows: Throughout this paper, we use the value η = 0.8, as is commonly assumed in computational models [53,55].This results in excitatory and inhibitory PSP sizes that differ by a factor of 4. For example, for q = 1 mV, we have q e = 0.5 mV and q i = 2 mV.

Refractory SoftPlus
The statistics of interspike intervals (ISIs) in the diffusion limit of equation ( 3) are found by solving a first passage time problem [34], but these results are mathematically complex, and do not apply to other neuron models, or even outside the diffusion limit (see appendix A).
For this reason, it is appealing to derive a mean-field model from an approximate transfer function that can be fitted numerically to multiple neuron models.In this section, we propose one called 'Refractory SoftPlus' after a function occasionally used in machine learning [58].SoftPlus is a smooth rectified linear function of one variable whose sharpness is controlled by a shape parameter β.It generalizes the ramp function x → max(x, 0), to which it converges in the infinite-β limit.
To define the Refractory SoftPlus transfer function, we augment SoftPlus with scale and shift parameters α and σ 0 .The input to SoftPlus is taken to be √ D = q √ R because the firing rate of LIF neuron in the diffusion limit scales as the square root of the diffusion coefficient [34].If the value of q is not available, it can be taken equal to unity, which will change the fitted parameter values, but not the results.For a transfer function which handles variable q, see appendix B.
We additionally model an absolute refractory period of duration T ref by inverting the output of SoftPlus to produce an ISI, adding the refractory period, then inverting again.This yields the following transfer function as a function of R with four shape parameters: For any given neuron model, and for fixed values of q and η, we compute the firing rate of the neuron over 100 simulated seconds for a range of 500 different values of R. For maximum comparability between different values of q, we use evenly spaced R calculated so that D ranges from 0 to 100 mV 2 s −1 .This cutoff point is arbitrary and our results do not depend on it; it is chosen only because it tends to drive our simulated neurons to fire in the tens of hertz, so that most of the shape of their transfer functions can be observed.The four parameters α, β, σ 0 , and T ref are then optimized to fit this simulated transfer function via nonlinear least squares.To make performance more comparable across models and conditions, we report fitting error normalized by the maximum firing rate observed in the simulation.For further details on numerical methods, see appendix C. Figure 2 shows an example of the performance of this transfer function in comparison to transfer functions from the literature.The simulation data is from the LIF neuron of equation ( 1), with V thresh − V rest = 15 mV, τ = 10 ms, and T ref = 2 ms, subject to Poisson inputs with η = 0.8 and q = 1 mV (giving q e = 0.5 mV and q i = 2 mV).It is compared to the analytical solution for the diffusion limit [34] as well as ReLU and sigmoidal transfer functions.These are simply univariate functions of √ D augmented with scale and shift parameters (see appendix C); for the ReLU, this is the ramp function defined above [58], whereas for the sigmoidal transfer function it is the hyperbolic tangent [21].In each case except the analytical solution, the transfer function is fitted to the left half of the data, i.e.the region where R<50 kHz.The full range of values is then used to show not only the performance of the original fit but also its ability to extrapolate.
The curve for the diffusion limit (black dashed) represents an analytical solution in the limit of infinitesimal PSP size (with constant D) and simulation time step.However, in this simulation, where both quantities are finite, the analytical solution loses substantial accuracy, resulting in a noticeable bias in the reported 'error' in the lower panel.The fitted ReLU (green) performs better, but suffers from an overestimate in the extrapolation regime, as well as a small downward bump caused by its sharp elbow.On the other hand, Refractory SoftPlus (red) is capable of capturing the functional form of the simulated transfer function, with virtually no error in the interpolation regime, and lower error in extrapolation as well.

Mean-field modeling
Next, we apply the Refractory SoftPlus fit of a neuronal transfer function to derive a mean-field model for the equilibrium firing rate of a recurrently connected population of M neurons, each of which fires as a renewal process at the rate F(R) when subject to a Poisson input with total rate R. We assume that each neuron is subject to background synaptic input at a total rate R bg , in addition to receiving inputs from a random set of N other neurons within the population.We exclude multiple and self connections from the present analysis.Note that R bg can be quite large compared to physiological firing rates (even in the tens of kHz) because it represents the combined activity of a large number of input neurons.The neurons within the population are grouped into an excitatory and an inhibitory subpopulation such that the excitatory fraction η is the same for the population as for the background input, and the PSP weights q e and q i are set as described in section 2.3 as well.
Under these conditions, the recurrent input and the background input will behave similarly, such that the neuron can be described by its transfer function using a total input rate r eff = R bg + Nr.This enables writing a consistency condition which can be solved to find the equilibrium firing rate of the population: Figure 3. Solutions to the mean-field consistency condition of equation ( 8) for LIF neurons with q = 5 mV and background noise rate R bg = 0.1 kHz.Solid curves correspond to different amounts N of recurrent connectivity, and give the postsynaptic firing rate as a function of the presynaptic firing rate; the mean-field system is at equilibrium when both are identical (black dashed line).
The stability of equilibria is indicated by the fill of the marker (stable = filled, unstable = unfilled, half-stable = half-filled).The model predicts a saddle-node bifurcation, where sufficient levels of recurrent connectivity result in a bistable system.
The solutions to this equation are visualized for several concrete values of N in figure 3.In this example, we numerically fitted the parameters of equation ( 7), as described in section 2.4, for LIF neurons with the large PSP size q = 5 mV and sparse Poisson background input R bg = 0.1 kHz.In each case, stability is defined by interpreting equation ( 8) as a map.If the absolute value of the derivative of the map is greater than 1 at the fixed point, iterated application will diverge and the fixed point is unstable [59].Stable fixed points have derivatives with magnitude less than 1, and half-stable fixed points are those where the absolute value of the derivative is exactly 1.
As shown in the figure, the number and location of solutions to the consistency condition depend on the value of the connectivity parameter N. When there is very little recurrent connectivity, the only stable firing rate is the spontaneous rate caused by the background input, close to r = 0 (blue).As N increases, the consistency curve gradually slopes upwards to contact the r = F(r eff ) line.This causes a saddle-node bifurcation at N = 51, which creates a half-stable fixed point at r = 33.9Hz (orange).As N continues to increase, the bifurcation continues, with the half-stable fixed point becoming one stable and one unstable fixed point (green).

Dynamics model
So far we have described only the equilibrium firing rate of a population of neurons, but it is typically of greater interest to investigate the response of a population to a time-varying external input.We do this using a simplified version of the Master equation approach [41].As in the prior work, we assume that the mean firing rate evolves according to linear dynamics with a fixed characteristic timescale τ [60].This yields the following simple dynamical model: The dynamics of equation ( 9) have stationary points at the zeros of its right hand side, which are exactly the fixed points of equation (8).In prior work on the master equation formalism, the timescale τ was considered somewhat arbitrary, but it was tightly coupled with the transfer function, and needed to be chosen carefully to yield accurate results [41,60].
However, in our model, τ appears only as the timescale of convergence, so model performance is insensitive to its value.As a result, τ can be set approximately based on a simple experiment such as a step response (as we do in this paper), or, when studying slow phenomena, simply specified to be faster than the largest frequency of interest.The only limitation is that τ should be large enough to enable the quasistatic assumption.The fitted transfer function describes the stationary firing rate of the neuron, which means that higher-order dynamical features such as SFA can lead to errors in the transient response of the model if the dynamics of interest are not fast compared to τ .

Transfer function fit performance
In order to assess the practicality of a data-driven mean-field model based on a Refractory SoftPlus fit, we first considered the amount of simulation data required to achieve an accurate fit.The computation time required to simulate a neuron's transfer function scales as O(ST), i.e. bilinearly with the number S of different input rates considered and the total simulated time T over which firing rates are calculated.Increasing either of these quantities increases the quality of the fit by increasing the precision or number of data points used.However, the discrepancy between Refractory SoftPlus and the true transfer function of the neuron also imposes a limit on the achievable accuracy.Therefore, we sought to minimize the computational costs of achieving the best possible quality of fit to the simulated transfer function.
Figure 4 demonstrates convergence of the fit for the LIF (blue), Izhikevich (orange), and the Hodgkin-Huxley (green) neurons.For each model, we simulated the transfer function once for large numbers S of input rates and once for large T, then considered subsets of the data.Fits were performed to subsets of the full simulation output, and the error of the fitted function was assessed on the final simulated transfer function.
The left side of figure 4 demonstrates convergence with the number S of input rates sampled.The PSP amplitude was fixed at q = 1 mV, and input rates ranged from 0 to R = 100 kHz.Fits were then calculated for evenly spaced subsets of this range.The fit error is already only a few percent at S = 5, and decreases to a final plateau value by S = 100.This can be attributed to the small number of parameters included in the fitted transfer function, which allows it to match observations without a large number of different data points.
The right side of figure 4 shows the convergence of the fit error as the duration of simulation considered for the calculation of firing rate for each neuron increases from T = 1 s to T = 1 × 10 4 s.A single simulation was performed of S = 100 different rates for a duration T = 1 × 10 4 s.Shorter simulation times were handled by considering only the beginning of the simulation data.This ensures that any startup transients and nonstationary effects, e.g.due to SFA, are included to the same extent that they would be in disjoint simulations.Again, fit error begins at only a few percent, and inclusion of more data in the fit soon reaches a point of diminishing returns.
Interestingly, for smaller T, the error of the fitted curve was substantially lower than the error of the short-time firing rate estimates to which the curve was fitted.The advantage continues for T up to 2.7 × 10 3 s for LIF neurons, 1.1 × 10 2 s for Izhikevich, and 7.4 × 10 2 s for Hodgkin-Huxley, with error below 0.5% in each case.(After this point, the error in the short-time firing rate estimates continue downwards as T −1/2 as suggested by appendix E.) This demonstrates that the fitted Refractory SoftPlus transfer function is able to reject a significant amount of noise in the empirical transfer function by providing averaging between adjacent input rates.

Generalization to other neuron models
In the majority of this paper, we use the simulator's default parameters, but varying these does not affect our conclusions.Indeed, the results can be generalized not only to other parameter values, but also to other  neuron models.To demonstrate this, we assessed the quality of the fit as the parameters of each neuron model were varied randomly according to the distributions given in table 1.In each case, dynamics of a single postsynaptic neuron were numerically simulated for 100 seconds subject to balanced Poisson input with η = 0.8 and q = 1.0 mV at 500 different total input rates R from 0 to 100 kHz.Performance is evaluated on the LIF neuron discussed above, together with the efficient quadratic integrate-and-fire neuron of Izhikevich [53] and the biophysical Hodgkin-Huxley neuron [2].The results of this comparison are shown in figure 5.
Crucially, Refractory SoftPlus is a good fit for the transfer function of all these models, despite their different dynamics and qualitative behavior.
For the LIF and Hodgkin-Huxley neurons, the distributions were derived from the simulator defaults by giving each time constant an exponential distribution with its default value as the mean, while the membrane capacitance was normally distributed about its default.For the Izhikevich neuron, the distribution of each parameter was made uniform over the range of values considered in the original description of the model [53].These choices of distribution are arbitrary, but they cover a broad parameter range for each model.For each of the models shown in figure 5, the fitting error under default parameters is indicated with a colored star in the histogram.This point is not an outlier, indicating that our results generalize well across neuron models for a broad range of parameters.
Note, however, that these results should be expected to depend on the dynamical regimes considered.In particular, bursting and strong adaptation are incompatible with any transfer function approach.The Figure 6.Agreement between numerical solutions to the consistency condition (lines) and observed equilibrium firing rates (markers) for M = 10 4 LIF neurons under three different combinations of connectivity parameters.For the monostable system, the model (black dashed line) and simulation (blue triangles) agree.For the bistable system, the discrepancy between the model (black solid line) and simulation (orange circles) is likely due to temporal correlations between the input neurons, as it can be eliminated by using annealed-average connectivity (green) The unstable fixed points are indicated with the red line.
qualitative dynamics of the LIF neuron are unaffected by its parameters, so the performance of the model is quite consistent in this case.On the other hand, the Izhikevich and HH neurons experience multiple dynamical regimes throughout their parameter space.In their default parametrization, all neurons considered are in a dynamical regime with weak or no adaptation, and do not experience bursts.This is addressed in more detail in appendix D.
Also note that for simplicity, the notation with which equation ( 2) describes the parametrization of synaptic weights to match results between PSP-based and PSC-based neurons ignores the fact that excitatory and inhibitory PSCs may have different parameters.In our simulations, not only were the PSP amplitudes q e and q i different because η = 1  2 , but also the synaptic time constants in the Hodgkin-Huxley neuron differ for excitatory and inhibitory synapses.However, despite these asymmetries and significant deviations from the input regime in which the approximations were derived, the flexibility of our assumptions has allowed fits to the simulated transfer function of the Hodgkin-Huxley neuron to be quite good as well.

Bifurcation analysis
Next we compared the bifurcation and equilibrium locations between the mean-field model and simulation.Activity fixed points were calculated for three different networks of M = 10 5 LIF neurons, as shown in figure 6. Simulations were performed at a range of values of N, each corresponding to a different network where each neuron receives connections from N presynaptic neurons sampled without replacement from the remainder of the population.In this and all other simulations where recurrent connectivity was included, recurrent synapses included an axonal delay chosen uniformly at random between 1 and 10 ms to represent the possibility of neurons spread out in space [1,61], but results were similar for small fixed delays as well.The effect of delays on population dynamics is discussed in appendix F.
Each displayed fixed point is the average firing rate across 10 simulations, each 2.5 s in duration.The lower stable population rate is computed by simply simulating the network from a zero initial condition subject to its background Poisson input.To find the higher stable population rate in the bistable network, we transiently raised the firing rate of the population using a 'warm-up' input intended to shift it into the basin of attraction of the upper fixed point.We did this by simulating the population without counting spikes for a 'warm-up time' of 1 s while the population was subject to additional excitatory Poisson input.To minimize the effects of this drive on the firing rate, we linearly ramped it down from R bg to 0 over that time.This was done for all values of N, not only those expected to be bistable.
The first condition pictured is a high-noise case which is monostable for the full range of N. Here, q = 3 mV, and R bg = 10 kHz.In this situation, the single firing rate fixed point increases slightly from the baseline spontaneous activity level at N = 0, but no interesting dynamics are observed.The model agrees well with simulation for the full simulated range, which is expected because recurrent input contributes under half of the total (ranging from 17% at N = 30 to 42% at N = 90).This means that there is enough background input to ensure that the assumptions of the model are satisfied even if the recurrent input deviates from Poisson statistics.
The second (orange) uses q = 5 mV and R bg = 0.1 kHz, the combination of parameters which produced a saddle-node bifurcation in figure 3.In this case, the accuracy of the mean-field model suffers somewhat.The network transitions from monostable to bistable between N = 50 and N = 55 as predicted, but the equilibrium firing rate in the active state is typically noticeably underestimated.This is likely caused by deviations in the input statistics from those used in calibrating the mean-field model.Across the upper branch of this condition, the fraction of input to the neuron which comes from the rest of the network rather than the Poisson background ranges from 97% to 99%.At the same time, firing rates increase, meaning that the absolute refractory period will make up a larger part of the interspike interval, and deviations from Poisson behavior will be more significant.
Indeed, the overestimation disappears entirely for annealed average connectivity.In this condition, rather than generating a fixed network where each neuron has N fixed presynaptic neurons, all-to-all connectivity is combined with probabilistic synapses, whose transmission probability N/M is chosen to match the input rate due to N presynaptic neurons on average [62].This breaks up temporal correlations and noticeably improves model predictions of equilibrium firing rates.
An additional caveat must be made specifically for the bistable networks with N just above the bifurcation.At N = 55, the unstable fixed point, which marks the boundary between the basins of attraction of the two equilibria, is at a relatively high firing rate.As a result, population activity fluctuations are often sufficient to move the network into the basin of attraction of the low-frequency equilibrium.Therefore, for N = 55 in particular, we rejected the simulations where this occurred, so the upper fixed point is based on only 1 data point for the second condition and 5 for the third.We did not observe the reverse basin hopping effect, where the network spontaneously jumps to the higher fixed point, because fluctuations near zero firing rate are much smaller.

Finite size effects
Although our mean-field model is based on finite connectivity, it does not explicitly depend on the total number of neurons M in the population.However, a sufficiently large M is implicitly assumed by requiring the neurons to be independent of each other.Two deterministic neurons will only fire independently if their inputs are independent.This can happen in two ways: (1) when M is very large, the inputs to each neuron will be disjoint, and therefore conditionally independent given the population firing rate; and (2) when R bg Nr, the input to each neuron is predominantly from the background, which is independent by construction.In this section, we investigate the effect of varying M on the dynamical and steady-state behavior of the model.
Figure 7 explores a moderate case with R bg and Nr on the same order of magnitude, by setting q = 3 mV, R bg = 10 kHz, and N = 75.For each of a range of values of M, ten different networks were simulated, each with the same overall parameters but different randomly generated connectivity.The binned firing rate was calculated as the number of spikes in each 1 ms bin of the 2 s recording, divided by M as well as the bin size.The result is a time-varying approximation of the firing rate in hertz, which follows a random walk with regression to the mean predicted by the model, as pictured.The standard deviation of the binned firing rate decreases proportionally to M −1/2 , as might be expected from the central limit theorem.This is illustrated using the function M → kM −1/2 with k fitted to the data by least squares (black dashed).The variation in the mean between simulations also decreases, but does not follow such a clear functional form.
Particularly surprising in light of past results showing substantial finite-size effects in population bursts [63] is the accuracy of the fixed point even when M = 100.For small M, input is highly correlated between neurons because the 75 synapses are chosen from only 100 neurons.This highly correlated input makes up about 37% of the total input to the neuron, but the location of the firing rate fixed point remains quite close to the model predictions.Evidently, the simple presence of background input on the same order of magnitude as the recurrent input adequately rescues the independence assumption even when the recurrent input itself is not close to independent.
Past mean-field models have assumed a specific form for the variance in the population rate, based on the binomial distribution [41,64].However, this assumption only applies to the fluctuating regime, where population activity is primarily input-driven; the results of these works consistently used less than 10% recurrent input.On the other hand, in the biologically common 'reverberating' regime, where population activity is driven substantially by recurrent inputs, correlated fluctuations tend to appear, which cannot be modeled in this way [62,65].Indeed, although at N = 0 the binomial prediction is recovered (not pictured), in this parameter regime, we observe an order of magnitude greater variance than would be predicted by such models.
Despite the insensitivity of the equilibrium location to M, an important finite size effect appears in the bifurcation analysis of figure 6.The simulations with N = 55 are close to a bifurcation, which means that the unstable fixed point marking the watershed between the lower and upper fixed points is 16 Hz.As a result, at M = 10 5 , the fluctuations in population rate are large enough that most simulations 'fall off ' the fixed point within a few seconds of simulation.In other words, stable equilibria predicted by equation ( 9) may be only metastable for finite M. Interestingly, it was recently proved that all-to-all coupled probabilistic integrate-and-fire neurons driven by Poisson noise have a unique stationary firing rate distribution [66].Thus, the appearance of metastability in our results suggests that finite-size effects are exacerbated by either hard spiking thresholds or sparse connectivity.
To investigate further, we define practical stability as the probability that a network will stay near the upper fixed point for the entire duration of the 2.5 s simulation used to calculate the firing rate in the previous section.This finite-size effect was recently investigated in detail in networks of quadratic integrate-and-fire neurons [67].Table 2 presents the practical stability of this particular case across a range of values of M.There is a clear trend towards higher practical stability as population rate fluctuations are reduced, with better practical stability for annealed-average than for static random connectivity.50 simulations were performed for each case in the table.

Dynamical response
Finally, we investigated the ability of the simple dynamics described in section 2.6 to model the sinusoid-following behavior of a recurrently connected network.In each case, a network of M = 10 3 neurons with q = 3 mV, N = 50, and R bg = 10 kHz was simulated subject to independent per-neuron Poisson background input as above, but in this case it is not stationary.Instead, the rate of the background input varies sinusoidally with time: We keep the frequency f bg = 1 Hz fixed and consider two different oscillation amplitudes A bg = 5 kHz and 10 kHz.The left part of figure 8 considers these two amplitudes (columns) across three different neuron models (rows).In the first case, the population firing rate also oscillates roughly sinusoidally, whereas when the background rate drops to zero, the effects appear much less linear.
The right side of figure 8 shows the RMS error of the left experimental condition as a function of the frequency with which the firing rate of the background input varies over time.The mean-field prediction is compared to the average of 50 independent simulations of the same experimental condition in order to avoid confounding by the substantial variation between independent simulations.The result is an error which remains relatively flat with frequency below 1 Hz, but which becomes substantially larger as frequency increases.This is to be expected, as the mean-field model is based on stationarity assumptions which cannot be expected to hold as the rate of population fluctuations approaches the slowest time constants of the individual neuronal dynamics.
It has been observed that the response of the population firing rate to fluctuations in the input is much faster than the response of any individual neuron to those changes [68].Indeed, we find that although the autonomous population rate fluctuations in figure 7 have a timescale in the tens of milliseconds, the population rate responds to a step input virtually instantaneously (not pictured).Indeed, past work suggests that the perturbation response of the population rate is limited by synaptic, not membrane, time constants [68].Therefore, we set τ = 1 ms in these simulations.Appendix F shows that changing τ has minimal effect for slowly fluctuating inputs.
Interestingly, the firing rate of the Izhikevich neurons appears to lead the prediction of the mean-field model, as if a derivative term is missing from the true dynamics.The derivative effect disappears under a set of parameters where SFA is absent (not pictured).This suggests that the discrepancy is caused by a failure of the model's stationarity assumption on short timescales.Since the recovery variable u is directly driven by firing, there can be a significant difference between the steady-state and transient response of the neuron.Indeed, a recent detailed mean-field model of all-to-all coupled Izhikevich neurons observed a similar effect which disappeared when SFA terms were replaced with their population mean [69].This limitation is not absolute, however.First, the small-deviation example is less problematic, because the population continues firing throughout the simulation in this case, so the adaptation level varies less.Also, the Hodgkin-Huxley neuron has larger-timescale effects due to dynamical variables which persist across multiple firings, and yet is modeled quite accurately in both examples.

Conclusion
We have developed an approximate mean-field model capable of finding the activity fixed points of a randomly connected network of neurons under the influence of synaptic inputs from other neurons in the network together with background Poisson input.This method is applicable to a wide variety of different neuronal models, and parameter values for the transfer function used in the mean-field approximation can be derived from small-scale and short-duration simulations.The resulting model predicts equilibrium and dynamical population firing rates in a randomly connected network as well as the location of a bifurcation from monostability to bistability.
The method of fitting a transfer function to simulations of the neuron under Poissonian input allows relaxing two key assumptions of past mean-field models.The diffusion limit is avoided by explicitly using the rate of finite-magnitude presynaptic events as input to the transfer function.Furthermore, since the transfer function is fitted numerically using only the neuron's spike times, our model does not utilize any information about the neuron's subthreshold dynamics.This is not entirely hypothesis-free, as some subthreshold dynamics may result in transfer functions that cannot be approximated well by Refractory SoftPlus, but our approach can just as easily incorporate a different transfer function in such cases.
Mean-field models often assume a particular distribution of the postsynaptic membrane potential, generally that it is normally distributed [20,41,70].However, realistic membrane potential distributions can be far from normal [43], making this so-called 'moment closure' assumption hard to justify.Another approach instead assumes a Lorentzian distribution, but this is incompatible with master-equation-type models due to the pathological behavior of the moments of the Cauchy distribution [26,69].Our approach avoids moment closure by instead assuming that the input to the neuron is sufficiently similar to the input used in calculating the original transfer function.
We observed an improvement in model accuracy in the most difficult parameter regime when the recurrent input within the population was made more Poisson-like by using annealed-average connectivity instead of a static random network.However, annealed-average connectivity is a strong assumption not typical in modeling; if accuracy in this regime is essential, it would be more useful to fit the transfer function under inputs more similar to those which will later be observed.Although the present work assumes Poisson inputs with loose excitatory-inhibitory balance, this assumption is not fundamental, and simulated transfer functions can in principle be calculated for other input models.This possibility is a key advantage of our approach.
The main limitations of our model are in the simplicity of the dynamics.Past mean-field models of firing rate variance have approached only the fluctuating regime [41,64], and generalization to the reverberating dynamics considered here is challenging.For this reason, we used only a first-order master equation, which does not quantify fluctuations from the mean, and so requires somewhat larger populations for full accuracy in some parameter regimes [60].Furthermore, the use of a transfer function requires assuming a stationary distribution of neuronal firing rate [60], but by making a quasistatic assumption, such models can also describe the dynamics of the firing rate [41].Furthermore, mean-field models can be augmented to represent other population dynamics such as adaptation variables or population heterogeneity, either by replacing dynamical terms with explicit population averages [71,72], or by expanding a model's existing usage of distributional models of single-neuron quantities [73].Incorporating these effects into the dynamics is a key direction of future work.
Like the results of previous authors, this approximation can be viewed as a bridge between computational neuroscience and machine learning perspectives on neural networks.Given that a similar functional form appears to be conserved across multiple spiking neural models, it is conceivable that it serves a function in biological neural networks as well as in artificial neural networks.In any case, theoretical developments may benefit from the existence of such a simple functional form which is robust across such a wide variety of neuronal dynamics.
had no role study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Figure 2 .
Figure 2. Fit and extrapolation performance of four potential transfer functions for the LIF neuron.Sigmoid, ReLU, and Refractory SoftPlus transfer functions are all fitted to half of the data (left side).Refractory SoftPlus performs well even in extrapolation (right side).Top: the four transfer functions plotted on top of simulated data (blue dots).The sigmoid appears adequate on the training data, but has high error in the extrapolation regime due to saturating.Bottom: the three fitted transfer functions and the simulated data plotted as 'error' from the diffusion limit.

Figure 4 .
Figure 4. Left: convergence of relative RMS error in the SoftPlus approximation with the number S of input rates sampled at a fixed simulation duration T = 100 s.The simulated transfer function for the full simulation is compared to Refractory SoftPlus functions fitted to subsets of the data.Right: convergence for S = 100 different rates sampled with increasing sample duration T.

Figure 5 .
Figure 5. Performance of Refractory SoftPlus fits across neuron models and parameter values.Top: comparison of fits to empirical transfer functions of M = 100 LIF (left), Izhikevich (middle), and Hodgkin-Huxley (right) neurons with default parameters, subject to Poisson inputs for T = 100 s with q = 1 mV.Bottom: distribution of the normalized residual of the fit over 1000 randomly parametrized instances of each neuron.Error for default parameters is indicated with a star.

Figure 7 .Table 2 .
Figure 7.The equilibrium firing rate does not depend on M, but time-correlated fluctuations scale as M −1/2 .Top: a single example of time-correlated fluctuations in population firing rate at M = 10 3 .Left: mean firing rate as a function of M, compared to the model prediction (gray).Averaged over 10 simulations, with one standard deviation shaded.Right: fluctuations in firing rate as a function of M, exhibiting clear scaling as large populations cancel finite size effects.A fitted inverse square root law is pictured for comparison.

Figure 8 .
Figure 8. mean-field model (black) is capable of predicting the dynamical response of recurrently connected networks to a time-varying stimulus.The leaky integrate-and-fire (blue) and Hodgkin-Huxley (green) networks are modeled quite accurately, whereas the Izhikevich network (orange) changes its firing rate in response to stimulus slightly ahead of the model prediction.

Table 1 .
Distributions of random parameter values used to generate figure 5.The notation Exp(µ) represents an exponential distribution by its mean, N (µ, σ) represents a normal distribution by its mean and standard deviation, and U [a,b] represents a uniform distribution over an interval.