Firing rate distributions in a feedforward network of neural oscillators with intrinsic and network heterogeneity

1.   Introduction

The underlying mechanisms of how cortical neural networks process stimuli is an important and complex topic. In a population of neurons, understanding how sensory signals are efficiently encoded and transmitted to higher areas of the brain is especially challenging given that neurons are known to be stochastic and have heterogeneous attributes. For this reason, and due to the nature of quantifying information content and efficient coding, many advances stem from a combination of experiments and computational modeling ^{[1,2,3,4,5,6]}. Firing rate heterogeneity is the most commonly studied heterogeneous attribute; it has been shown to have consequences on neural coding in the olfactory bulb ^[6,59], on the visual system (i.e., orientation tuning curves ^[7,53]), and on a variety of other systems including the auditory ^[18] and electrosensory systems ^[4,29]. In general, the population firing rate and the level of its heterogeneity directly affect important information-theoretic measures of coding such as the Fisher information and mutual information ^[35]. Firing rate heterogeneity is a fundamental entity for systems that code signals based on rate, or the total number of spikes. Although we do not focus on the potential impact of higher order spiking statistics (i.e., covariance ^[2]) or spike timing in this paper, firing rate heterogeneity has deep implications on the way cells encode sensory signals.

Several studies have shown that increased firing rate heterogeneity is a signature of more efficient coding (e.g., Fisher Information ^[6,53], discrimination error of distinct trial-averaged firing rates ^[4]). However, the relative level of heterogeneity is important because other studies have shown that too much or too little heterogeneity can be detrimental. This has been shown, for example, in spiking model networks ^[8], in a generalized linear model (GLM) fit to olfactory bulb neurons ^[59], and in orientation tuning curves ^[9]. Thus, the actual level of heterogeneity, not just whether neurons are heterogeneous, is important to quantify. Analyzing firing rate heterogeneity in model networks has the potential benefit of enabling predictions of how neural attributes relate to efficient coding.

An understudied aspect of firing rate heterogeneity is the genesis of this network statistic from underlying neural attributes. Many theoretical studies specify the firing rate heterogeneity a priori, but interactions of cellular and network attributes are nonlinear and can lead to unexpected spiking dynamics ^{[3,10,11,12,13]}. However, others have focused on how intrinsic attributes affect neural coding ^[1,14]. Importantly, experimental measurements have shown that both intrinsic and network heterogeneity are material. At the cellular level, many intrinsic factors influence the firing rate of a cell such as ion channel composition or cell morphology. Our modeling centers on neural oscillators, a class of cells whose intrinsic properties lead to repetitive action potentials. We primarily consider intrinsic heterogeneity via a neuron's phase-response curve ( PRC) ^[15], an experimentally measurable entity that quantifies how inputs advance/delay time to next spike in oscillating neurons. Experiments have shown that PRCs are heterogeneous in olfactory bulb cells ^[16] as well as in mice visual cortex ^[55]. The researchers ^[55] demonstrate that the PRC type can modulate with acetylcholine (carbachol) rather than altering stimulus input. At the network level, the same presynaptic input strength can result in heterogeneous postsynaptic responses. Effective synaptic input strength is determined by many physiological parameters: presynaptic firing rate, postsynaptic potential size ^[23,38] for each presynaptic spike, or in the total number of inputs ^[38,48,49] (in addition to these cortical examples, experiments in the cerebellum show postsynaptic targets are heterogeneous ^[25,65]). We do not distinguish here between these different aspects of network heterogeneity; rather we use a simple parameter for input strength to model network heterogeneity (see ^[7,12], who modeled network heterogeneity in a similar way).

In this study, we focus on the nonlinear interactions of intrinsic and network heterogeneity because such theoretical studies are relatively rare (i.e, ^{[6,8,16,59,66]} considered intrinsic heterogeneity alone, and ^[7,33,49] considered network heterogeneity alone). Some exceptions are Marder and colleagues ^[11,12], who have considered many heterogeneous attributes in detail, including intrinsic and network, with a focus on the detailed structure of voltage traces with relatively small numbers of neurons. Both intrinsic and synaptic diversity in two coupled oscillators have been considered in dynamic clamp experiments combined with theoretical modeling in ^[32,54] (see ^[30] for a more thorough review). At least in the context of altering the spiking statistics of large networks of neurons with noise, the only other studies that we are aware of that consider the effects of both intrinsic and network heterogeneity are: ^[3,10].

In ^[10], we analyzed how the relationship between intrinsic (spike threshold) and network heterogeneity, along with the operating regime, alter firing rate distributions in a generic recurrent network model of leaky integrate-and-fire (LIF) neurons. In ^[3], we made stimulus-dependent predictions about the relationship between threshold and synaptic input, as well as effective network connectivity, with modeling and in vivo recordings of the hindbrain in weakly electric fish. However, these works do not include realistic intrinsic dynamics when modeling the relationship between intrinsic and network heterogeneity. Simple spiking models such as LIF do not account for ion channel dynamics that can have significant and counter-intuitive effects on firing rate statistics (i.e., blocking a potassium current does not always increase firing rates ^[41,43]). In addition, higher dimensional models that incorporate these features can better characterize the intrinsic dynamics of neurons, as opposed to LIF. LIF and many spiking models only integrate network input without properly accounting for how the current state affects responses.

Here we consider a feedforward network of heterogeneous Morris-Lecar cells ^[45], where "feedforward" indicates that target neurons do not affect the presynaptic population. We find the relationship between intrinsic (PRC) and network heterogeneity can dramatically alter firing rate heterogeneity or distribution. We use a phase reduction method to calculate asymptotic firing rates assuming weak noise and coupling. Our first order analytic formulas qualitatively and succinctly capture the observed firing rate statistics of the high-dimensional model. We also demonstrate the value of the second order analytic formulas in capturing how firing rates change as the noise level varies. Our work reveals which aspects of the PRC control firing rate dynamics in the weak coupling and noise regime.

2.   Materials and method

2.1.   General feedforward oscillator model

Consider a population of $N$ distinct neural oscillators receiving independent noise, coupled to a presynaptic population providing feedforward input. Let $X_j\in\mathbb{R}^n$ denote the state variables of the $j^{\text{th}}\in\{1, 2, \dots, N\}$ oscillator. The uncoupled and noise-free oscillatory system would be characterized as $\frac{dX_j}{dt} = F_j(X_j)$ with a period $T_j$. With the addition of noise and coupling, the equation for $X_j$ is of the form:

where ${\bf \varsigma}\ll1$ is the power of the noise, $q_j w_{j, k}G\ll 1$ is the coupling, and ${\vec \xi}_j(t)$ are independent white noise processes with zero mean and variance $\langle \xi(t)\xi(t') \rangle = \delta(t-t')$. Combining noise and coupling perturbations, we define the system as $\frac{dX_j}{dt} = :F_j(X_j)+\epsilon_j$ and assume the models are It$\hat{\hbox{o}}$ stochastic differential equations (SDEs). We also assume that the unperturbed system $\frac{dX_j}{dt} = F_j(X_j)$ has an asymptotically stable limit cycle, $X_0(t) = X_0(t+T_j)$.

We use a phase reduction method applied to model above (Eq. 2.1) to obtain a network of phase oscillator models ^[17,26]. The phase reduction employed here is standard and has been described previously by many authors (see Chapter 8 of ^[28] and Chapter 1 of ^[15]). Note that there are other phase reduction methods; for example, see ^[56] for addressing how the time-scale of the noise and return time to limit cycle could yield a different phase reduction; see ^[52,57] for phase descriptions with noise and without a stable limit cycle; with moderate to strong perturbations, see ^[21,46,47]. For our purposes, we will focus on the case where the oscillators return to the limit cycle very quickly with weak perturbations.

Near the limit cycle, there is a function $\phi: \mathbb{R}^n \rightarrow S^1$ mapping a neighborhood of the limit cycle to the one-dimensional phase along the limit cycle, $\Theta_j \in [0, T_j)$, where $T_j$ is the period of the uncoupled oscillator. Defining $\Theta_j = \phi(X_j)$, the variable $\Theta_j$ satisfies

This analytically exact equation has to be further approximated because $X_j$ is not captured with the phase variables. However, since $\epsilon_j$ is small, we can approximate $X_j$ by $X_0(\Theta_j)$ to get a simpler equation for $\Theta_j$:

where $Z(\Theta) : = \nabla_X \phi(X_0(\Theta))$. The function $Z(\Theta)$ is called the adjoint and satisfies the linear equation

A normalization condition uniquely determines the solution to this equation (see ^[15]). We assume that the noise and coupling only affect the voltage component (i.e., voltage is the first component so $\epsilon_j$ has 0 except in the first), which holds for a wide class of neuron models. Thus, for the vector product (i.e., last term of RHS of Eq. 2.3), the only relevant quantity is the first component of $Z(\Theta)$, which we call $\Delta_j$, the (infinitesimal) Phase Resetting Curve or PRC of the neuron. The PRC is a periodic function that has negligible value at the end points because in neurons, perturbations have negligible effect on the dynamics at the moment of a spike, and is also an experimentally measurable entity. The PRCs here are are calculated with the program XPPAUT ^[27].

The result of applying the phase reduction and scaling Eq. 2.3 by $T_j$ is the following set of stochastic differential equations:

where again the noise is independent, frequency $\omega_j = \frac{1}{T_j}$ because the phase variables $\Theta_j$ have been scaled to take on values in [0, 1), and

The first factor $\frac{1}{T_j}$ is from the scaling. $\overline{Y_k}(t)$ is the average value of presynaptic input. ${\bf \Delta}_j(t)$ is the time-scaled version of the PRC ${\bf \Delta}_j(t): = T_j\Delta_j\left(\frac{t}{T_j}\right)$. Here the noise level ${\bf \varsigma}$ is simply scaled by the period: $\sigma = \frac{{\bf \varsigma}}{T_j }.$ The term $\frac{\sigma^2}{2} \Delta_j(\Theta_j)\Delta'_j(\Theta_j)$ is a result of Itô's Lemma because the original model (Eq. 2.1) is an Itô stochastic differential equation (SDE). We model intrinsic heterogeneity with different $\Delta_j$ (although there is minuscule variation in $\omega_j$ in the models we consider), and network heterogeneity is captured by $q_j$ that scales the presynaptic input.

2.2.   Morris-Lecar model

The specific multi-dimensional model we focus on is the Morris-Lecar model ^[45]. Although this model has only 2 state variables (dynamics of voltage trace and potassium gating), it is an ideal model for our purposes because our central focus is on intrinsic heterogeneity manifested with different PRCs. Indeed, a wide variety of PRCs result with varying 4 parameters (Figure 3C, D). The model dynamics consist of multiple time-scales, three intrinsic ionic currents, and spike generation. At its core, the dynamics of the voltage trace are modeled by the difference between an applied background current and the intrinsic ionic currents, so $\frac{dV_j}{dt} \propto I_{app} - I_{ion}$. Meanwhile, potassium dynamics are modeled by the difference between the steady-state fraction of open potassium channel at the current voltage level and the actual current fraction of open channels, so $\frac{dW_J}{dt} \propto W_{inf}(V_j) - W_j$. The model is:

$\xi_j(t)$ is an independent white noise process with strength ${\bf \varsigma}$ and $C_m$ is the membrane capacitance. The parameters $g_L$, $g_K$, and $g_{Ca}$ are the maximal conductances for the leak current, potassium channels, and calcium channels, respectively; $E_L$, $E_K$, and $E_{Ca}$ are the corresponding reversal potentials. The functions in the model are:

where $m_\infty$ and $W$ model the steady-state fraction of open calcium and potassium channels, respectively, and $\tau_W$ represents the time scale of potassium; $V_1$ through $V_4$ are the parameters that affect the intrinsic ion channel dynamics. All of the intrinsic parameters are fixed except $V_3$, $V_4$, $\phi$ and the background current $I_{app}$ (see Table 1). These four parameters are varied to model intrinsic heterogeneity, resulting in different PRCs (see Figure 3C). In some circumstances, $\phi$ is calculated via the temperature ^[28], but here we assume it can vary within the population.

Finally, the sum ${q_j}\sum\limits_{k = 1}^M {{w_{j, k}}} {s_k}(t)({V_j} - {E_{sy{n_k}}})$ represents the total pre-synaptic input assumed to be feedforward from an unmodeled population. The parameter $q_j$ models network heterogeneity; the values are chosen independently from a uniform distribution from positive values: $q_j \sim \mathcal{U}(q_{min}, q_{max})$, with either $[q_{min}, q_{max}] = [0.005, 0.015]$ or $[q_{min}, q_{max}] = [0.001, 0.003]$. The same presynaptic input strength can result in heterogeneous postsynaptic responses (see ^[7,12] who modeled network heterogeneity in a similar way). The coupling matrix $w_{j, k}$ is an $N \times M$ matrix with $N = 730$ neurons (population of interest) and $M = 1000$ synapses, $k\in\{1, \dots, M\}$, consisting entirely of 0's and 1's. The coupling matrix is an Erdös-Rényi graph, with the probability of connection of 0.3 for all $w_{j, k}$ (independently chosen). The random number of presynaptic inputs to the $N = 730$ neurons is a minor source of variability but not one of our significant sources of heterogeneity. The synapse variable $s_k(t)$ is modeled by the following ODE system:

The times $t^*$ when $a$ instantaneously jumps $a(t^*)\rightarrow a(t^*)+A$ are random, governed by a homogeneous Poisson process with rate $\lambda$.

The synaptic kernel associated with the synapse model is:

where $H(x)$ is the Heaviside step function (1 if $x > 0$, and 0 if $x < 0$). Here we assume $\tau_d > \tau_r$, commonly observed in cortical synapses.

2.3.   Phase reduction of the feedforwad Morris-Lecar network

The result of applying the general phase reduction method in section 2.1 to the feedforward Morris-Lecar network with conductance-based synaptic input is described here.

The noise level ${\bf \varsigma}$ in the Morris-Lecar model is simply scaled by the period so that:

Recall that the time-scaled version of the PRC, ${\bf \Delta}_j(t)$, is defined as the following: ${\bf \Delta}_j(t): = T_j\Delta_j\left(\frac{t}{T_j}\right)$. As described previously for a general model (Eq. 2.5), we scale time by $T_j$ and use the dimensionless version of the PRC, $\Delta_j$. The feedforward coupling term (Eq. 2.6) is approximated by:

where $V^{LC}_j(t)$ is the unperturbed voltage in one cycle. The PRC is taken out of the integral as an ad-hoc approximation, although partially justified because the feedforward input is very noisy (synapses are driven by a Poisson process) and is effectively pulse-coupled inputs when $E_{syn_k}$ is much larger/smaller than the possible values of $V_j(t)$ (see ^[31,40] who use phase oscillator models of a similar form). Finally, the synaptic input $s_k(t)$ is approximated by applying the synaptic filter (Eq. 2.11) at all discretized points in the period and averaging because the synapses are driven by a Poisson process. Thus,

where $\langle \cdot \rangle_{t'}$ denotes the average over all $t'\in[0, T_j)$.

2.4.   Asymptotic approximation to the firing rate distribution

Here we describe the phase reduction theory to capture the firing statistics of the full Morris-Lecar feedforward network. The phase reduction (Eqs. 2.13–2.15) makes the subsequent asymptotic calculations feasible. Stochastic systems are often characterized by a probability density equation which is described by a Fokker-Planck ^[51] equation. Let

The corresponding Fokker-Planck equation of the j$^{\text{th}}$ neuron is:

with periodic boundary conditions: $\varrho_j(\theta = 0, t) = \varrho_j(\theta = 1, t)$ and normalization $\int_0^1 \varrho_j(\theta, t)\, d\theta = 1$. The firing rate of actions potentials is:

The probability flux $J_j(\theta, t)$ can re-written by differentiation to get:

We are only interested in the steady-state solution ($\frac{\partial \varrho_j}{\partial t} = 0$):

Since $J_j(\theta)$ is constant, we pick a convenient $\theta$, namely $\theta = 1$, which gives the firing rate of the $j^{\text{th}}$ neuron. This value of $\theta$ has the added benefit that the PRC vanishes there at the moment of firing: $\Delta(1) = 0$. Thus, the firing rate is simply:

We exploit the weak coupling and noise assumption by applying a standard asymptotic approximation to the P.D.F. $\varrho_j$ to obtain an approximation for the firing rate of the $j^{\text{th}}$ neuron, and consequently the firing rate heterogeneity (i.e., standard deviation across all $r_j$):

The P.D.F. has to integrate to 1 (i.e., $\int_0^1 \varrho_j(\theta)\, d\theta = 1$); we assume the following normalization conditions for the asymptotic approximation of $\varrho$:

so that any truncation from 0$^{\text{th}}$ order and onward results in the correct normalization. Substituting Eq. 2.21 into Eq. 2.19 while assuming $\sigma^2$ and $q_j\bar{P}_j$ are order $\epsilon$ gives:

This results in a hierarchy of asymptotic equations (after dividing by $\epsilon^l$). We impose the condition that each order of the approximation for $J_j$ has 0 derivative (with respect to $\theta$), equivalent to assuming the steady-state Fokker-Planck Eq. 2.19 is satisfied in all orders. The $0^{\text{th}}$ order equation is simply:

With the normalization condition: $\int_{0}^{1} \varrho_0(\theta)\ d\theta = 1$ and the fact that $J_j\equiv K_0$ is constant (i.e., $J_j$ has 0 derivative) results in $\varrho_0\equiv 1$, giving the following 0$^{\text{th}}$ order approximation to $r_j$:

First Order Formula

The $1^{\text{st}}$ order formula is derived via the $\mathcal{O}(\epsilon)$ equation in Eq. 2.23. Note that:

Since $\varrho_j'(\theta)$ is multiplied by an order $\epsilon$ term ($\frac{\sigma^2}{2}$), the resulting $1^{\text{st}}$ order approximation for $J_j$ does not contain derivatives of $\varrho_l$:

Setting $J_j$ to a constant $K_1$ via $0 = \frac{\partial}{\partial\theta}\{J_j(\theta)\}$ and substituting for $\varrho_0$ results in a simple equation for $\varrho_1$:

The constant $K_1$ is determined by integrating both sides from $\theta = 0$ to 1 and using the normalization condition: $\int_{0}^{1} \varrho_1(\theta)\ d\theta = 0$.

Thus we have:

where: $\bar{\Delta}: = \int_{0}^{1} \Delta (\theta) d\theta$. The resulting 1$^{\text{st}}$ order approximation is:

Second Order Formula

We derive a 2$^{\text{nd}}$ order approximation to the P.D.F. and firing rate analogously. Again we use Eq. 2.23 but only consider $\mathcal{O}(\epsilon^2)$ this time.

From the equation for $\varrho_1$ (Eq. 2.28), we have:

Substituting $\varrho_1(\theta)$ and $\varrho_1'(\theta)$ into the RHS of Eq. 2.30:

We skip some of the straight forward yet tedious calculations. Assuming all order of $\varrho_j$ satisfy the equation $0 = \frac{\partial}{\partial\theta}\{J_j(\theta)\}$ as before, the equation for $\varrho_2$ is:

To distinguish the contribution of noise, network inputs, etc., the terms are grouped by $(\frac{\sigma^2}{2})^2$, $(q_j\bar{P}_j)^2$, and the interaction between these two: $\frac{\sigma^2}{2} q_j \bar{P}_j$:

We solve for $\varrho_2(\theta)$ and $K_2$ by integrating both sides in $\theta$ from 0 to 1, using the normalization condition (Eq.2.22) and recalling that $\Delta_j(0) = \Delta_j(1) = 0$. Then, we have:

where $\overline{\Delta^2}: = \int_{0}^{1} \Delta^2(\theta) d\theta$. This equation can be further simplified with integration by parts:

Eq. 2.33 is equal to:

The final result for $\varrho_2$ using Eqs 2.32 and 2.34 is:

Evaluating at $\theta = 1$ greatly simplifies $\varrho_2$, the only value that matters for the firing rate.

The $2^{\text{nd}}$ order approximation to the firing rate then is:

Finally, we consider the standard deviation across the population predicted by our analytic theory: $\sigma(\vec{r}) = \frac{1}{N}\sum_{j = 1}^{N} (r_j-\mu{\vec{r}})^2$. We only use the $1^{\text{st}}$ order approximation (Eq. 2.29) in the firing rate heterogeneity approximation – in the Results section the $2^{\text{nd}}$ order formula is not useful for the large feedforward network.

(all variables in the standard deviation $\sigma$ are vectors).

3.   Results

We start with an illustrative example that highlights the importance of properly accounting for intrinsic dynamics. Consider the Morris-Lecar (ML) model with two different parameter sets (Figure 1). The voltage trajectories are similar and the phase planes have minor discrepancies that are predominant in the $W$ variable, so one would naively think that the population response would be similar. However, we will show that the population response can be different with two simple homogeneous networks in Figure 2. Although experimentalists interested in systems neuroscience and coding tend to focus primarily on network coupling structure and synaptic dynamics, we show that cellular attributes alone can shape response statistics and subsequently provide a mathematical explanation of these observations.

Consider two homogeneous populations of ML cells, each having different intrinsic dynamics (Figure 1 and Eqs. 2.7–2.8) which we call SNIC (saddle-node bifurcation on invariance circle) and Hopf. We use these names (SNIC and Hopf) because the stable limit cycles emerge from these respective types of bifurcations. Each homogeneous population, of size $N = 730$, receives feedforward input that is heterogeneous, modeled by a factor $q_j$ (i.e., network heterogeneity) that scales input from a feedforward population of $M = 1000$, ${q_j}\sum\limits_{k = 1}^{1000} {{w_{j, k}}} {s_k}(t)({V_j} - {E_{sy{n_k}}})$; cells in the target population are uncoupled. The coupling between the presynaptic and postsynaptic populations is randomly and independently chosen (Erdös-Rényi graph), with the probability of connection of 0.3 for all $w_{j, k}$; all feedforward networks in this paper have this structure (all except Figure 6). These two network models have vastly different population firing rate statistics with statistically identical input (Figure 2). Not only are the mean firing rates different, but the standard deviation of the population firing rate vary more with network type (SNIC or Hopf) than the actual type of input (excitatory versus inhibitory input, Figure 2A, B). The SNIC network has consistently larger firing rate heterogeneity than the Hopf network. The mean firing rate also changes more significantly between excitatory and inhibitory inputs for SNIC cells than it does for Hopf cells. This surprising result can be explained with our theoretical calculations (see below). Note that the bifurcation diagrams for each network (Figure 2C, D) are very different, providing a hint to a possible explanation of the observations; one network type has cells that undergo a SNIC bifurcation and the other following a Hopf bifurcation with increased current injection.

Biophysical neuron models often have many dimensions, and in a coupled network, analysis is often intractable because of the dimensions and nonlinear processing of presynaptic inputs. The assumption that both coupling and noise are weak, however, allows us to employ the phase reduction while still incorporating biophysically-realistic dynamics. The phase-resetting curve ( PRC) is a key entity in the phase reduction. We use the common convention for the PRC that an excitatory input at phase $\theta$ will delay the time until the next spike if $\Delta(\theta) < 0$ but will advance the time until the next spike if $\Delta(\theta) > 0$ (see Figure 3A). These advances or delays with perturbations naturally impact firing rate. We derive approximations to the individual firing rates: $1^{\text{st}}$ order Eq. 2.29. Specifically, the approximation to firing rate heterogeneity (standard deviation across the population) is:

This is derived from Eq. 2.38, noting that frequencies ($\omega_j$), integrated PRCs ($\bar{\Delta}: = \int_0^1\Delta_j(\theta)\, d\theta$), and average input strength times noise ($\sigma(\bar{P}\, q)$) are all the same within a population. Thus, independent of input type, the consistently lower firing rate heterogeneity in the Hopf population occurs because the integral of the PRC $\bar{\Delta}$ is smaller compared to the SNIC network; recall that all other components of the prior equation are fixed. The $\bar{\Delta}$ are shown in Figure 3D – the smallest value (dark blue) is the Hopf and the largest value (dark red) is the SNIC. Note that how limit cycles emerge has deep implications on the shape of the PRC, i.e., Hopf generally has large negative region versus SNIC has mostly positive PRC ^[26]. SNIC is often called Type I because the firing rate vs. input current curve is continuous, and Hopf is often called Type II because the firing rate vs. input current curve has a discontinuous jump.

Next we consider an ML network with both intrinsic heterogeneity via the PRC and network heterogeneity in the same manner as before with a parameter $q_j$. Experimentalists often focus on either measuring intrinsic cellular attributes (uncoupled) or synaptic dynamics and connectivity structure because of feasibility in recordings and precise control. However, we show here that not only do both cellular and circuit attributes matter, but they have to be taken together and not separately because they interact nonlinearly to affect firing rate statistics. To this end, we developed 68 parameter sets to model intrinsic heterogeneity by systematically varying four parameters: effective input current and three others related to recovery variable dynamics (see Table 1). The effective input current is varied to ensure all cells have an approximate period of $85\, \text{ms}$ when varying the other three parameters (see all limit cycles in Figure 3B and PRCs in Figure 3C). The PRC of a prototypical SNIC or Type I oscillator ^[26] would result from parameters similar to the following: $V_3 = 12\, \text{mV}$, $V_4 = 17\, \text{mV}$, and $\phi = 0.06667\, \text{ms}^{-1}$. Meanwhile, the PRC of a prototypical Hopf or Type II oscillator ^[26] would result from parameters on the opposite extremes: $V_3 = 1\, \text{mV}$, $V_4 = 31\, \text{mV}$, and $\phi = 0.04\, \text{ms}^{-1}$. We produced intermediate PRCs by uniformly selecting the three parameter values from the ranges bracketed by these prototypical values and then adjusting current to establish periodic uniformity. Figure 3 also shows the integrals of the 68 PRCs $\bar{\Delta}$ (Figure 3D), which as we have already seen, is important for firing rate heterogeneity. We remark that there are notable differences in the both the PRCs and its integral $\bar{\Delta}$.

The network in Figure 4 is from duplicating the unique set of 68 parameter sets to obtain a larger network of $N = 730$ with a wide range of intrinsic dynamics; the network is feedforward with unmodeled 'pre-synaptic' neurons that only provide excitatory presynaptic input. This network has both intrinsic heterogeneity via the PRCs (Figure 3C) and network heterogeneity via the $q_j$ parameter that scales the network input.

We perform Monte Carlo simulations to accurately capture the firing rate heterogeneity. With both sources of heterogeneity chosen independently, we see significant heterogeneity (measured by the standard deviation across the entire population (Figure 4A)) even though the intrinsic periods are roughly equal $\omega_j\approx\omega$ (Figure 3C). Recall the firing rate heterogeneity approximation Eq. 2.38 ($\sigma(\omega+q\, P\, \bar{\Delta})$) that has the term:

Our calculations predict that the relative level of firing rate heterogeneity can be modulated by changing the relationship, or correlation, between components of the intrinsic $\bar{\Delta}$ and network heterogeneity $q$ (neglecting $P$, average synaptic input). Specifically, if $q_j$ is negatively correlated with $\bar{\Delta}_j$, the firing rate heterogeneity will be relatively smaller since the product of large and small numbers will not deviate much cross the population $j = 1, \dots, N$. Contrast that when the correlation is positive: the product of large with large and small with small numbers willl result in larger heterogeneity across the population. Note that given a vector of $\bar{\Delta}_j$ values, we can generate $q_j$ so that the vector of $q$ and $\bar{\Delta}$ have any specified (Pearson's) correlation coefficient, see ^[3,10]. Indeed, this prediction is precisely what is shown in the Monte Carlo simulations (Figure 4B on a cell-by-cell basis and Figure 4D); these trends are captured by both the $1^{\text{st}}$ order (Eq. 2.29) and $2^{\text{nd}}$ order (Eq. 2.37) approximations.

Our $1^{\text{st}}$ theory is valuable, producing significantly less error than the $0^{\text{th}}$ order theory that just incorporates intrinsic frequency (Figure 4C). We can see from Figure 4A–B that the error of cell-to-cell firing rate approximations increases with actual firing rate, which is not surprising since the deviations from the uncoupled system ${\bf \varepsilon}$ are relatively large. However, the trends are still captured qualitatively with the phase reduction theory, and the cell-to-cell error appears to be largely independent of correlation (Figure 4B). From Figure 4C, we see that the firing rate heterogeneity approximation is better with $\rho(\bar{\Delta}, q) < 0$, and the error increases with $\rho$. The $2^{\text{nd}}$ order approximation (Eq. 2.37) does not improve the predictive power; in fact it is slightly less accurate yet more complicated than the $1^{\text{st}}$ order approximation.

We next applied our theory with the same heterogeneous ML neurons but with presynaptic input from both excitatory and inhibitory conductance-based inputs. The Poisson process input rates are scaled to three levels: $\lambda_E = 0.2$ and $\lambda_I = 0.1$, $\lambda_E = 0.4$ and $\lambda_I = 0.2$, and $\lambda_E = 0.8$ and $\lambda_I = 0.4$. In Figure 5A and B, we show the cell-to-cell comparisons for the $1^{\text{st}}$ order approximation (Eq. 2.29) with Monte Carlo simulations, including just the two extreme correlation values, $\rho = -0.7$ and $\rho = 0.7$. We see from Figure 5C–D that the asymptotic approximation becomes very inaccurate as firing rates increase, but it still qualitatively captures the mean firing rate trends (Figure 5C) and the firing rate heterogeneity trends (Figure 5D). In Figure 5C–D, with a fixed correlation, our approximation captures how the firing rate changes with presynaptic input rate (i.e., varying 1X, 2X, 3X with a fixed color); also, the approximation captures how the mean and standard deviation both increase with correlation (fixing synaptic input rate, gray lines connecting stars and squares).

Thus far, the $2^{\text{nd}}$ order approximation to the firing rate has not been any more informative than the $1^{\text{st}}$ order approximation. This result was initially quite puzzling because one would expect that the $2^{\text{nd}}$ order approximation to be more accurate than $1^{\text{st}}$ order, at least for small, perhaps minuscule parameters. However, even with very small parameters, the $2^{\text{nd}}$ order approximation does no better in any of the networks we considered. Upon further investigation, we find that both terms in the $2^{\text{nd}}$ order approximation (see Eq. 2.37 or Eq. 2.40) are positive for all of the networks thus far, so no matter how small the parameters are, this observation will hold since the 1$^{\text{st}}$ order approximation overestimates all firing rates. It would seem that the $2^{\text{nd}}$ order approximation yields no information, however note that it can likely account for differences in how the noise parameter $\sigma$ or other aspects of the PRC ($\Delta^2$, $\Delta'$, $\bar{\Delta^2}: = \int \Delta^2(\theta)\, d\theta$) affect firing rate. Thus, we demonstrate the value of this calculation on two single phase oscillator SDEs where we have more precise control over the components (i.e., $\omega_j, q_j\bar{P}_j, \sigma$), rather than having them endowed from the full coupled network. Consider the simpler equation (based on Eq. 2.13):

where we set the intrinsic frequency to 1 and treat the network heterogeneity (formerly $q_j\bar{P}_j$) as a specified parameter $Q$. The $2^{\text{nd}}$ order approximation is then (see Eq. 2.37)

We select two cells where the integrals of the PRCs ($\bar{\Delta}$) are similar but the dynamics of the PRCs are drastically different (Figure 6A). These two PRCs are obtained from the full ML network but scaled to have (arbitrary) units of phase. Consequently, while our $1^{\text{st}}$ order theory (Eq. 2.29) would predict similar firing rates, our $2^{\text{nd}}$ order theory (Eq. 2.40) might capture the discrepancies between the two cells via the additional terms. In the bar charts of Figure 6B, note that the $Q$ and $Q^2$ terms in Eq. 2.40 do not vary much, while there is a relatively large discrepancy in the $\sigma^4$ term. Thus, our theory indicates that the variation in the noise/coupling term $\sigma$, rather than the network strength $Q$, drives the difference in firing rate of the cells. As $\sigma$ increases, our $2^{\text{nd}}$ order calculation predicts that Cell II will then fire at a faster rate than Cell I (Figure 6D). Our Monte Carlo simulation corroborates the theory (Figure 6C), with Cell I's firing rate only mildly increasing with the increase in noise and coupling while Cell II's firing rate increases at a far more significant rate. The theory qualitatively matches the results for both cells up to approximately $\sigma = 25$, at which point our assumption of weak noise and coupling has been violated because $\sigma \neq \mathcal{O}{(\epsilon)}$. Therefore, the $2^{\text{nd}}$ order firing rate approximation is valuable in capturing firing rate dynamics despite quantitative inaccuracies in our models.

4.   Discussion

Analyzing the modulation of firing rate heterogeneity in model neural networks could potentially link how neural attributes (cellular and circuit) are related to efficient coding. As we have seen in Figure 2, seemingly similar populations of neurons receiving the same synaptic input can have significantly different firing rate means and standard deviations. This statistical discrepancy demonstrates that we cannot ignore intrinsic heterogeneity. We then consider a feedforward network with combined intrinsic heterogeneity and network heterogeneity in order to evaluate their nonlinear interaction. Our asymptotic calculations, developed from phase reduction analysis, succinctly capture the impact of how sources of heterogeneity lead to firing rate heterogeneity in several feedforward networks: with excitatory input only and with different levels of balanced excitatory and inhibitory input.

Using our theory, we can effectively modulate the level of firing rate heterogeneity by altering the correlation structure between sources of intrinsic and network heterogeneity in a biophysical (Morris-Lecar) model. The elegant $1^{\text{st}}$ order theory demonstrates the impact of phase-resetting curves on firing rate statistics and its nonlinear interaction with network input. Although the calculations are less accurate if the weak noise and coupling assumptions are violated, they still capture the overarching trends of intrinsically oscillating neurons in these feedforward networks. The $2^{\text{nd}}$ order theory is not any more informative than $1^{\text{st}}$ order in our large feedforward networks. However, the $2^{\text{nd}}$ order theory can account for the impact of noise on firing rate dynamics in single cell models (Figure 6), giving an overall more complete asymptotic theory. Unlike many common spiking models (e.g., leaky integrate-and-fire), we have incorporated biophysically realistic intrinsic dynamics with significant heterogeneity.

The spike threshold of the cell (or equivalently, the intrinsic frequency of the phase oscillator) is a crucial attribute that controls output firing rates, since low threshold (high intrinsic frequency) will directly cause high firing rates and vice-versa. Thus, variable spike threshold or frequency is a common way to induce intrinsic heterogeneity ^[3,8,10,44]. Threshold heterogeneity has been shown experimentally in cortical cells ^[22] and has crucial effects in the electrosensory system ^[44] and others ^[50]. We did not focus on spike threshold as an intrinsic heterogeneous attribute because the interaction between network heterogeneity and heterogeneity in these variables was unsurprising and had been analyzed before by ^[10]. A low effective network strength correlated with high spike threshold would of course lead to lower firing rate heterogeneity, while the opposite would lead to higher firing rate heterogeneity. These results are fairly intuitive and consistent with our results, so we focused on the relationship between effective network strength and intrinsic phase-dependent response to inputs because the results are not obvious. Thus, accounting for cellular response over a cycle, beyond just the spike frequency, should be a consideration to experimentalists to account for different neural responses.

Although the oscillating regime is not always applicable for neural dynamics, there are many times when neural networks can at least be well-modeled by oscillators (also see ^[64]). Neural oscillator models have already been successfully used to model the olfactory system ^[16,66] as well as breathing ^[24] and locomotion ^[58]. There are other mathematical frameworks that account for oscillatory dynamics ^[21,46,57] generated from excitable cells ^[39,48] that we have not used here because of our phase oscillator assumption, but may be effective in analyzing (intrinsic) firing rate heterogeneity.

Neural networks were the underlying motivation for this work, but the general framework of using phase reduction methods to study output heterogeneity from nonlinear interaction of heterogeneous components might generalize to other fields. Indeed, there is a long history of using phase reduction methods in biology and other fields. Winfree pioneered the approach for biological systems of synchronized oscillators, with a particular focus on circadian rhythms ^[60,61]. Kuramoto's model of coupled oscillators is usable for such phenomenon as the Belousov-Zhabotinsky chemical reaction ^[17,36,37]. The call and response of two tree frogs, an anti-phase synchronized system, can also be modeled via phase oscillators ^[19]. In engineering, the leg movements of passively walking robots ^[42], electric circuits powering flashing LEDs with periodicity ^[20], and injection locking of neighboring electric circuits and lasers ^[63] all have limit-cycle oscillations.

The exact relationship of various heterogeneous components, such as the correlation of $q$ and functions of $\Delta$, is not precisely known, so we systematically varied these entities in this modeling study. However, ^[34] labeled and classified 1600 neurons and their connectivity profiles in the mouse visual cortex in vitro. Some other experimental studies suggest that heterogeneous attributes are not random but have structure ^[11,12,62]. Thus, technological advances may lead to opportunities to experimentally measure these relationships. Our theoretical study explores the consequences on spiking statistics with idealized relationships to provide a better sense of how heterogeneity and relationships between forms of heterogeneity shape neural responses. This theoretical work complements studies of how entities that are difficult to measure experimentally (i.e. intrinsic and network diversity) affect firing rate activity, which is easier to measure. Such theories are necessary to ultimately understand how the relationship of neural attributes affect neural coding.

Acknowledgments

Cheng Ly is supported by a grant from the Simons Foundation # 355173.

Conflict of interest

All authors declare no conflict of interest in any capacity.