Emergent stochastic oscillations and signal detection in tree networks of excitable elements

TL;DR

This paper investigates how stochastic oscillations and signal detection emerge in tree networks of excitable elements, revealing how network topology influences spike train statistics and sensory signal processing.

Contribution

It demonstrates that in strongly-coupled tree networks, spike train statistics can be predicted from a single excitable element with rescaled parameters, and shows how topology tuning affects signal detection.

Findings

Spike train statistics can be predicted from an isolated element with rescaled parameters.

Network topology influences firing rate and variability.

Optimal input discrimination can be achieved by tuning network topology.

Abstract

We study the stochastic dynamics of strongly-coupled excitable elements on a tree network. The peripheral nodes receive independent random inputs which may induce large spiking events propagating through the branches of the tree and leading to global coherent oscillations in the network. This scenario may be relevant to action potential generation in certain sensory neurons, which possess myelinated distal dendritic tree-like arbors with excitable nodes of Ranvier at peripheral and branching nodes and exhibit noisy periodic sequences of action potentials. We focus on the spiking statistics of the central node, which fires in response to a noisy input at peripheral nodes. We show that, in the strong coupling regime, relevant to myelinated dendritic trees, the spike train statistics can be predicted from an isolated excitable element with rescaled parameters according to the network topology. Furthermore, we show that by varying the network topology the spike train statistics of the central node can be tuned to have a certain firing rate and variability, or to allow for an optimal discrimination of inputs applied at the peripheral nodes.

Equations

$$N=\frac{d^{G+1}-1}{d-1}\text{and}{N}_p=d^G,$$

$$C{\dot{V}}_k=-I_{\text{ion}}[V_k,{\mathbf{u}}_k]+\kappa \sum _{j=0}^{N-1}{A}_{k,j}(V_j-V_k)+I_{{\mathrm{ext}}_k},$$

$${\dot{\mathbf{u}}}_k=\alpha (V_k)(1-{\mathbf{u}}_k)-\beta (V_k)(1-{\mathbf{u}}_k),$$

$$\kappa =\frac{a}{4lL\rho },$$

$$I_{{\mathrm{ext}}_k}={\delta }_{k,p}[I+\sqrt{2D}{\xi }_p(t)],$$

$$r=\frac{1}{⟨\Delta t_j⟩},C_V=r\sqrt{⟨(\Delta t_j-⟨\Delta t_j⟩{)}^2⟩},$$

$$d^{\prime }=2\frac{|{\mu }_T(I)-{\mu }_T(I+\Delta I)|}{{\sigma }_T(I)+{\sigma }_T(I+\Delta I)},$$

$$J_{\text{LB}}(I)=\frac{1}{{\sigma }_T^2(I)}{\left(\frac{d{\mu }_T}{dI}\right)}^2,$$

$$d^{\prime }\approx \Delta I\sqrt{J_{\text{LB}}(I)}.$$

$$I_{\text{th}}/I_{\text{b}}=R^{-1}(G,d),$$

$$R(G,d)=\frac{N_p}{N}=\frac{d^G(d-1)}{d^{G+1}-1},$$

$$C{\dot{V}}_0\approx f[V_0,{\text{u}}_0]+I_{\text{eff}}+\sqrt{2D_{\text{eff}}}{\xi }_{\text{eff}}(t),$$

$${\dot{\mathbf{u}}}_0\approx g[V_0,{\text{u}}_0].$$

$$C{\dot{V}}_k=f[V_k,{\text{u}}_k]+\kappa \sum _{j=0}^{N-1}{A}_{k,j}(V_j-V_k)+I_{{\mathrm{ext}}_k},k=0,1,...,N-1.$$

$$⟨V{⟩}_g:=\frac{1}{d^g}\sum _{g\text{th shell}}{V}_k.$$

$$C⟨\dot{V}{⟩}_g$$

$$C\Delta {\dot{V}}_g=\left(⟨f(V,\text{u}){⟩}_g-⟨f(V,\text{u}){⟩}_{g+1}\right)+\Delta I_g+\Delta {\xi }_g(t)+\{\begin{array}{ll}\kappa \left(d\Delta V_1-(d+1)\Delta V_0\right),& g=0,\\ & \\ \kappa \left(\Delta V_{g-1}-(d+1)\Delta V_g+d\Delta V_{g+1}\right),& 0<g<G-1,\\ & \\ \kappa \left(\Delta V_{G-2}-(d+1)\Delta V_{G-1}\right),& g=G-1,\end{array}$$

$$\Delta f_g^0=0,\kappa \gg |\Delta f_g^1|,g=0,1,...,G-1,$$

$$C\frac{d}{dt}\mathbf{\Delta }\text{V}\approx \kappa \mathbf{B}\mathbf{\Delta }\text{V}+\mathbf{\Delta }\text{I}+\mathbf{\Delta }\mathbf{\xi }(t).$$

$$\mathbf{\Delta }\text{V}=(\Delta V_0,\Delta V_1,...,\Delta V_{G-1}{)}^T,\mathbf{\Delta }\text{I}=(\Delta I,\Delta I_1,...,\Delta I_{G-1}{)}^T,\mathbf{\Delta }\mathbf{\xi }(t)=(\Delta {\xi }_0(t),\Delta {\xi }_1(t),...,\Delta {\xi }_{G-1}(t){)}^T,$$

$$\mathbf{B}=\left(\begin{array}{ccccc}-\left(d+1\right)& d& 0& ...& 0\\ 1& -\left(d+1\right)& d& ...& ...\\ 0& 1& -\left(d+1\right)& ...& 0\\ ...& ...& ...& ...& d\\ 0& ..& 0& 1& -\left(d+1\right)\end{array}\right).$$

$$\mathbf{\Delta }\text{V}\approx -\frac{1}{\kappa }{\mathbf{B}}^{-1}\left(\mathbf{\Delta }\text{I}+\mathbf{\Delta }\mathbf{\xi }(t)\right),$$

$$C⟨\dot{V}{⟩}_0=C{\dot{V}}_0=f[V_0,{\text{u}}_0]+d{\mathbf{B}}^{-1}{\left(\mathbf{\Delta }\text{I}+\mathbf{\Delta }\mathbf{\xi }(t)\right)}_1.$$

$$I_{\text{eff}}=d{\left({\mathbf{B}}^{-1}\mathbf{\Delta }\text{I}\right)}_1,D_{\text{eff}}=d{\left({\mathbf{B}}^{-1}\mathbf{\Delta }\mathbf{\xi }(t)\right)}_1.$$

$$I_{\text{eff}}=\frac{N_p}{N}I=R(G,d)I,$$

$$D_{\text{eff}}=\frac{N_p}{N^2}D=\frac{R(G,d)}{N}D,$$

$$\dot{m}$$

$${\dot{h}}_k$$

$${\alpha }_m(V)$$

$${\beta }_m(V)$$

$${\alpha }_h(V)$$

$${\beta }_h(V)$$

$$I_X$$

Full text

Introduction

Coupled noisy excitable systems serve as relevant models for a wide range of natural phenomena, including pattern formation in chemical reactions [1, 2] and in social networks [3, 4, 5, 6], dynamics of gene regulatory networks [7] and of single and networked neurons [8, 9, 10]. Networks of noisy excitable elements exhibit a rich variety of spatio-temporal dynamics, depending on the strength and topology of coupling and the noise intensity [11, 12, 13]. For example, the coherence of emergent network oscillations can be controlled by modifying the noise intensity, the coupling strength, or by changing the network size or topology [14, 15, 16, 17, 18, 19]. The dynamic range and sensitivity of complex networks of excitable elements to external stimuli can by optimized for critical topologies [20, 21, 22].

In the present paper, we focus on the dynamics of regular tree networks of strongly coupled excitable elements which receive random and independent excitations to their peripheral nodes, as sketched in Fig. 1. Our study is motivated by the morphology of certain peripheral sensory neurons, which possess branched myelinated dendritic terminals at their receptive fields, with multiple nodes of Ranvier. Their extended terminal branching resembles the dendrite structure of neurons in the central nervous system (CNS) [23, 24]. Myelinated segments form a tree-like structure with nodes of Ranvier at each branching point. Myelination terminates at peripheral nodes of Ranvier, called heminodes, which receive sensory signals. Thus, such sensory neurons may possess multiple spike initiation zones at heminodes which encode a local sensory signal into a stream of action potentials (APs) which are merged into a single output spike train transmitted to the CNS [25]. Examples for such neurons are the afferent innervation of muscle spindles [26, 27, 28, 29], pain receptors [30], cutaneous mechanoreceptors [31, 32], and lung receptors [33]. Interestingly, sensory neurons with myelinated dendrites may exhibit spontaneous activity characterized by coherent periodic spiking, despite that their peripheral heminodes presumably receive uncorrelated noisy excitations [28]. Figure 1 then can be viewed as a model for a branched myelinated dendritic terminal, where peripheral nodes receive uncorrelated stochastic inputs and are linked by myelinated segments. Due to the high density of Na+ ion channels at the nodes of Ranvier, APs may be excited independently at different peripheral nodes. High electrical conductivity of myelinated segments, which link the individual nodes, result in a strong coupling between the nodes. Therefore, their stochastic dynamics synchronizes. This may result in noisy periodic spiking of the primary branching (central) node, as we have shown for star networks of excitable elements [34].

Here we use a biophysical model for nodes of Ranvier connected by myelinated links on regular trees and show numerically and analytically that the collective response of the network can be deduced from the stochastic dynamics of a single effective node with parameters scaled according to the network size and topology. Thus, our study allows for the prediction of the stochastic network dynamics from the tree topology. We then discuss how the tree topology affects the firing statistics of the central node and the discriminability of input signals.

Model and Methods

Discrete Cable Model

In the present paper, we study the stochastic dynamics of excitable elements linked on a regular tree (see Fig. 1). Branching starts at the primary (central) node (number 0 in Fig. 1) and continues through several generations. Only the peripheral nodes receive external inputs. Referring to a model of branched myelinated dendrites, these peripheral nodes are called heminodes and receive inputs from thin unmyelinated processes (neurites). APs are initiated at the heminodes and then propagate on the tree towards the primary branching node and eventually to the CNS.

Here we consider regular trees whose topology is characterized by two parameters: the branching, $d$, and the number of generations, $G$. Given these two parameters, the total number of nodes, $N$, and the number of peripheral nodes, $N_{p}$, are given by

[TABLE]

respectively. The dynamics of the membrane potential is approximated by a discrete cable model [35] in which nodes of Ranvier are connected by passive resistive links according to the network topology. All active nodes and passive links are assumed to be identical, except that peripheral nodes receive external inputs. The membrane potential $V_{k}(t)$ of the $k$th node obeys the dynamics

[TABLE]

where the index $k=0,1,2,...,N-1$ marks the respective node. In particular, $k=0$ refers to the central node. In Eq. (2) the term $I_{ion}[V_{k},{\bf u}_{k}]$ stands for nodal ionic currents and ${\bf u}_{k}(t)$ is a vector whose components are the gating variables of the nodal ion channels and $C=2\leavevmode\nobreak\ \mu$F/cm2 is the nodal capacitance per area. In the following we use two particular models for the nodes of Ranvier: a Hodgkin-Huxley-type (HH) model with Na+ and leak currents [34] and the Frankenhaeuser-Huxley (FH) model which includes additional K+ and persistent Na+ currents. The HH nodal model includes two gating variables, $m$ and $h$, for Na+ channels, i.e. $I_{\text{ion}}[V,\mathbf{u}]=I_{\text{ion}}[V,m,h]$. The FH model includes two additional gating variables, $n$ for K+, and $p$ for persistent Na+ channels: $I_{\text{ion}}[V,\mathbf{u}]=I_{\text{ion}}[V,m,h,n,p]$. The detailed equations and parameters of the nodal models are provided in the Supplementary Material.

The coupling term in Eq. (2), $\kappa\sum_{j=0}^{N-1}A_{k,j}(V_{j}-V_{k})$, contains the adjacency matrix $\mathbf{A}$ of the tree graph and the coupling strength, $\kappa$, in units of Siemens per area. Its value can be calculated from the sizes of the node and myelinated links, and the axoplasmic resistivity:

[TABLE]

where $a$ is the diameter of the node (and of links), $l$ is the nodal length, $L$ is the length of connecting links and $\rho$ is the axoplasmic resistivity. For example, for $\rho=100$ $\Omega$cm, the nodal diameter and length $a=10$ $\mu$m, $l=1$ $\mu$m, and the length of myelinated segment $L=200$ $\mu$m, the coupling strength is $\kappa=1250$ mS/cm2. This provides a biophysically-plausible range for $\kappa$, which we use as a control parameter in the following.

The external current $I_{ext}$ is applied only to the peripheral nodes and consists of a constant part $I$ and noisy part, i.e.

[TABLE]

where $p$ denotes indicies of peripheral nodes; $\delta_{k,p}$ is the Kronecker delta; $D$ scales the intensity of the Gaussian white noise $\xi_{p}(t)$, which is uncorrelated for different peripheral nodes, $\langle\xi_{i}(t)\xi_{j}(t+\tau)\rangle=\delta_{i,j}\,\delta(\tau)$. Thus, peripheral nodes receive random uncorrelated inputs.

Equations (2) were integrated numerically using explicit Euler–Maruyama methods with timestep of 0.1$\mu$s.

Variability of Generated Sequences of Action Potentials

Our primary interest is the statistics of a spike train generated by the central node. A spike is identified as a full-size AP with a magnitude of at least 60 mV. We extracted a sequence of spike times, $t_{j}$, at the central node from 60 – 120 s long simulation runs. The corresponding sequence of interspike intervals (ISIs) $\Delta t_{j}=t_{j+1}-t_{j}$, is characterized by the mean firing rate, $r$ and the coefficient of variation, $C_{V}$ as,

[TABLE]

where the average is taken over all ISIs in the spike train of the central node.

Signal Detection

To characterize the signal detection capacity of a tree network, we considered a small constant stimulus, $\Delta I$, applied to the peripheral nodes in addition to the stimulus $I$, and calculated a normalized distance between resulting spike count distributions of the central node with and without this addition. Such a measure of distance is given by the discriminability, $d^{\prime}$, defined as [36],

[TABLE]

where $\mu_{T}$ and $\sigma_{T}$ are the mean and standard deviation of the spike count in a time interval $T$, respectively. The discriminability quantifies how well the network responses to two different stimuli, $I$ and $I+\Delta I$, can be distinguished by observing corresponding spike count statistics at the central node.

The discriminability is related to the Fisher information, which provides the theoretical limit of how accurately a stimulus $I$ can be estimated by observing a spike train [37]. For the spike count statistics, a lower bound of the Fisher information can be written as [36],

[TABLE]

and is related to the discriminability, $d^{\prime}$ by [36],

[TABLE]

Larger values of the Fisher information refer to more accurate estimation of the stimulus from the spike train and so to better discrimination between two stimuli $I$ and $I+\Delta I$.

The discriminability Eq. (6) was calculated by collecting spike counts of the central node for $5000$ independent time intervals of lengths $T=200$ ms, and calculating the mean and standard deviation for two values of the stimulus, $I$ and $I+\Delta I$, applied to the peripheral nodes of a tree network [36]. We also calculated the lower bound of the Fisher information Eq. (7) for the single uncoupled node as a function of the input current (stimulus) $I$ and the noise intensity, $D$, using a similar numerical procedure.

Results

Emergence of Periodic Firing in Deterministic Tree Networks

At first, we consider the case of a deterministic input, $D=0$. In the absence of the external input, $I_{\rm ext}=0$, an isolated node is in the excitable regime. A sufficiently high constant current, $I_{\text{AH}}$ results in a subcritical Andronov-Hopf bifurcation of the equilibrium state rendering an isolated node to fire a periodic sequence of APs. The corresponding limit cycle disappears in a saddle-node bifurcation for a lower external current, $I_{\text{SN}}$. For the HH nodal model the saddle-node bifurcation occurs at $I_{\text{SN}}\approx 28.15\mu$A/cm2 and the subcritical Andronov-Hopf bifurcation at $I_{\text{AH}}\approx 29.06\mu$A/cm2, so in a narrow range $I_{\text{SN}}<I<I_{\text{AH}}$ an isolated node is bistable, possessing a stable equilibrium and a stable limit cycle. When the nodes are coupled on a tree network and external currents are applied to the peripheral nodes, the dynamics of the network may become quite complex. For example, in case of weak coupling, peripheral nodes fire APs, which fail to propagate to the central node, so that nodes in the inner generations of the network exhibit small-amplitude spikes. For a stronger coupling, nodes in the inner generations may fire APs, but with skipping relative to APs in the periphery, demonstrating various $m:n$ synchronization patterns. However, for strong coupling and sufficiently high external currents the network shows fully synchronized periodic firing.

A comprehensive analysis of the deterministic dynamics is beyond the scope of this study. Instead, since our primary interest is in the emergence of periodic sequences of full-size APs at the central node, we address the following question: Given the tree topology, $G$ and $d$, and the coupling strength, $\kappa$, what is a threshold value $I^{\text{th}}$ of a constant current applied to peripherals, $I$, which makes the central node to generate repetitive firing of full-size APs? To this end, we perform simulations of tree networks with given $\kappa$, $G$ and $d$. Initially membrane potentials of the individual nodes are randomly distributed around the stable equilibrium of an isolated node for $I=0$. Then we apply a current $I>0$ and determined the minimal value, $I^{\text{th}}$ of $I$ at which the central node generated APs repetitively at steady state. Results are shown in Fig. 2.

At the central node periodic firing of APs occurs for values of $I$ and $\kappa$ above the corresponding curves in Fig. 2. Below these curves, the network is excitable in the sense that no repetitive firing of APs is observed at the central node. In the following, we refer to these two regimes as oscillatory (repetitive firing of full-size APs by the central node) and excitable (no repetitive firing of APs by the central node). The threshold value of the external current, $I^{\text{th}}$, increases for weak and moderate values of the coupling strength. Consequently the network needs stronger external input to the peripheral nodes to sustain periodic firing of the central node.

Figure  2 shows two distinct coupling regimes. For weak coupling, $\kappa<2$ mS/cm2, the threshold current $I^{\text{th}}$ is independent of the network size, i.e. the number of generations, $G$, and branching, $d$. In contrast, for strong coupling, $\kappa>60$ mS/cm2, the threshold current saturates, and its value increases with increasing number of generations. This is illustrated further in Fig. 3 showing the threshold current vs the number of generations for strong coupling. Note that the strong coupling regime spans the range of realistic coupling strengths for models of branched myelinated dendrites. As can be seen in Fig. 3, the threshold current follows a characteristic dependence saturating for trees with a large number of generations, $G$, and decreases with the increase of branching, $d$. Apparently, this dependence follows the scaling relation:

[TABLE]

where $I_{\text{b}}$ is a bifurcation value of the constant current in the isolated single node and the scaling factor $R(G,d)$ is the ratio of the number of peripheral nodes to the total number of nodes,

[TABLE]

whith $N_{p}=d^{G}$ and $N=\sum_{k=0}^{G}d^{k}$, for regular trees. This scaling relation, Eq. (9), holds for strongly-coupled trees and is derived below.

Deterministic trees with the FH nodal model show similar dynamics with the same scaling as in Fig. 3.

Stochastic Dynamics

The addition of uncorrelated noise to the peripheral nodes allows for the generation of APs in the excitable regime. Fig. 4 shows an example of the stochastic dynamics for a tree with $G=5$ generations and $d=3$ branching. In the excitable regime ($I=20\leavevmode\nobreak\ \mu$A/cm2) noise of sufficient intensity induces APs in peripheral nodes. For weak coupling ($\kappa=0.3$ mS/cm2) noise-induced APs in adjacent generations are not synchronized (superimposed spikes for peripheral nodes fill densely corresponding generation panels) and do not propagate beyond the 2-nd generation, which shows only sparse APs. Increasing the coupling strength leads to progressive synchronization of nodes in adjacent generations and finally results in the generation of APs in the central node. For strong coupling the whole network fires almost in synchrony. We note, however, that even for strong coupling outer generations show some spike jitter. We also note that strong coupling leads to slower and more random firing of APs.

As observed for star networks[34], the dynamics of the central node in a tree network depends non-monotonously on the coupling strength. As shown in Fig. 5, there exist optimal, rather small values of the coupling strength for excitable and oscillatory trees at which fastest (maximum firing rate) and most coherent (minimal coefficient of variation, $C_{V}$) firing is observed, respectively. For extremely weak coupling APs, which are fired by different peripheral nodes, are not synchronized and fail to propagate to the central node (Fig. 4, upper left panel). Increasing the coupling strength leads to stronger interaction between the branch nodes and results in synchronous firing of all nodes.

However, the size of a tree, i.e. the number of generations, is critical for the firing statistics of the central node. Furthermore, excitable and oscillatory trees demonstrate qualitatively different behaviour in the biologically-relevant strong coupling regime. In excitable trees, firing of APs becomes slower and more irregular if the coupling is strengthened and trees with more generations are considered. For large $G$ and strong coupling firing stops [Fig. 5(a1)] since excitatory inputs to peripheral nodes are too weak to sustain firing of APs. In contrast, in oscillatory trees, the firing rate saturates for strong coupling [Fig. 5(a2)] and firing becomes more regular if strongly-coupled trees with more generations are considered [Fig. 5(b2)].

Scaling of Effective Current and Noise intensity

In the strong coupling regime the dynamics of the central node of a tree network can be described by the dynamics of a single isolated node with membrane potential $V_{0}(t)$ and with effective input current $I_{\text{eff}}$ and an effective Gaussian noise with intensity $D_{\text{eff}}$, i.e. the influence of the coupling term on the dynamics can be approximated by a constant current and a white Gaussian noise. Then the dynamics of the membrane potential of the central node in eq. (2) can be approximated by

[TABLE]

In the following, we derive those effective parameters for regular trees of diffusively-coupled nodes.

In the network model, the dynamics of the membrane potential of $k$-th node is given by

[TABLE]

In order to derive approximations for the scaling of the effective current $I_{{\rm ext}}$ and the noise intensity $D_{\text{eff}}$, we extend the approach of Kouvaris et al. [38], who considered the propagation of excitable waves in a tree network of identical Fitz-Hugh Nagumo nodes in the absence of noisy inputs. Following their approach, we consider the dynamics of the average membrane potential $\langle V\rangle_{g}$ (termed density by Kouvaris et al.) in each shell in a tree. Here and in the following $\langle\rangle_{g}$ denotes averaging over all nodes of the $g$th shell,

[TABLE]

The dynamics of those densities can be obtained by averaging the respective equations for the dynamics of the membrane potentials, Eq. (12), over all nodes in one shell. Since the total number of connections between nodes in shell $g$ and $g-1$ is $d^{g}$, we obtain

[TABLE]

Note that, since peripheral nodes are subject to independent white Gaussian noises, the corresponding equation for the averaged membrane potentials of the peripheral generation contains white Gaussian noise $\xi_{G}(t)$ with reduced intensity $\frac{D}{N_{p}}$.

Since the coupling terms depend only on the difference between densities of the membrane potentials in adjacent generations $\Delta V_{g}:=\langle V\rangle_{g}-\langle V\rangle_{g+1}$, we consider the dynamics of those differences next. Subtracting equations for $\langle\dot{V}\rangle_{g}$ yields,

[TABLE]

where $\Delta I_{g}=-\delta_{g,G-1}\,I$ and $\Delta\xi_{g}(t)=-\delta_{g,G-1}\,\sqrt{2\frac{D}{N_{p}}}\,\xi_{G}(t)$ is Gaussian white noise.

Next, we consider the case of strong coupling. In that case, $\Delta V_{g}$ becomes small, and the membrane potentials of individual nodes approach the average potentials of the corresponding shell. Thus, we can approximate $\langle f(V,\textbf{u})\rangle_{g}-\langle f(V,\textbf{u})\rangle_{g+1}$ by a Taylor expansion around $\Delta V_{g}=0$, i.e. $\langle f(V,\textbf{u})\rangle_{g}-\langle f(V,\textbf{u})\rangle_{g+1}\approx\Delta f_{g}^{0}+\Delta f_{g}^{1}\Delta V_{g}+h.o.$. It then follows for strong coupling, i.e.

[TABLE]

that the dynamics of the averaged potential is dominated by the coupling term and $\Delta V_{g}$ can be approximated by a multidimensional Ornstein$-$Uhlenbeck process,

[TABLE]

Here we introduced the $G$-dimensional vectors,

[TABLE]

and the $G\times G$ tridiagonal Toeplitz matrix,

[TABLE]

In the strong coupling limit (16), deviations of $\boldsymbol{\Delta}\textbf{V}$ from its mean value decay extremely fast and we can use an adiabatic elimination [39] to approximate $\boldsymbol{\Delta}\textbf{V}$ by its mean value plus a white Gaussian noise. Both, the mean voltage difference and the intensity of the Gaussian white noise in the strong coupling limit can be obtained by setting the left-hand side of Eq. (17) to zero. This yields

[TABLE]

where $\mathbf{B}^{-1}$ is the inverse of the matrix $\mathbf{B}$. In order to obtain an approximation for the dynamics of the central node, we can use Eq. (19) to replace $V_{0}-\langle V\rangle_{1}$ by $\Delta V_{0}$ in Eq. (Scaling of Effective Current and Noise intensity) for the central node, $g=0$. This yields

[TABLE]

Here and in the following the index ”$1$” denotes the first component of a $G$-dimensional vector. Next, the effective parameters $I_{\text{eff}}$ and $D_{\text{eff}}$ can be obtained by comparing Eqs. (20) and (11). This yields the effective input current and the intensity of the effective white Gaussian noise,

[TABLE]

For the special case, considered in this study, that only peripheral nodes are subject to noisy inputs, i.e. $\boldsymbol{\Delta}\textbf{I}=(0,0,...,-I)^{T}$, and $\boldsymbol{\Delta}\boldsymbol{\xi}(t)=(0,0,...,-\sqrt{2D/N_{p}}\,\xi_{G}(t))^{T}$, the calculation of the effective parameters $I_{\text{eff}}$ and $D_{\text{eff}}$ requires only a single component, $(1,G)$, of the inverse matrix, $\mathbf{B}^{-1}$. Since $\mathbf{B}$ is a tridiagonal Toeplitz matrix, we can apply the results of Ref. [4] to calculate this component (see Supplemental Material for details on calculations) and find for the effective current,

[TABLE]

and for the effective noise intensity,

[TABLE]

where the scaling factor $R(G,d)$ is given by Eq. (10).

Investigating the scaling of the effective parameters in more detail, we first note that our theory yields the scaling relation, Eq. (9), observed for the deterministic threshold current in Fig. 3. In fact, the same scaling relation applies to the bifurcation values of $I$ in the deterministic model, e.g. the subcritical Andronov-Hopf bifurcation of the equilibrium or the saddle-node bifurcation of the limit cycles. Second, in Fig. 6, we demonstrate the validity of the theoretical scaling predictions by comparing results for the mean firing rate and the CV from direct simulation of tree networks with those from a single node (11) with input current and noise intensity scaled according to Eqs. (22) and (23), respectively. As illustrated in Fig. 7, we find an excellent correspondence of both results. This indicates that in the strong coupling limit the response of the network can be predicted from the stochastic dynamics of the effective central node. The statistics of interspike intervals for a single isolated node versus input current parameters, i.e. constant component, $I$, and noise intensity, $D$, can be easily computed numerically yielding two-dimensional maps, such as shown in Fig. 7. Then for a given size (number of generations, $G$) and branching, $d$, of a tree, the scaled parameters, Eqs. (22) and (23), set an operation point for the tree on the parametric map of a single element. Thus, predictions of the firing statistics of the central node of a tree of strongly-coupled excitable elements can be deduced by superimposing parametric dependencies of $I_{\text{eff}}$ and $D_{\text{eff}}$ on the parameters of the network. Figure 7 demonstrates this for trees of strongly-coupled HH nodes in the oscillating regime. In trees with more generations $G$ the operation point is shifted towards smaller currents and lower noise intensities, resulting, for oscillatory trees, in slower and more coherent firing of the tree’s central node. Finally, in the strong coupling limit the scaling relations (22, 23) are independent of the particular choice of the nodal model, e.g. they are expected to work for either Hodgkin-Huxley or Frankenhaeuser-Huxley nodal models.

Signal Detection

The signal detection efficiency of a neuron can be quantified using the discriminability and the Fisher information [41, 36, 42, 43, 44]. In case of our model of coupled excitable elements on a tree, we use these measures to characterize how the tree topology affects its ability to distinguish between two stimuli, $I$ and $I+\Delta I$, applied to the peripheral nodes.

The preceding section showed that in the strong coupling limit, the stochastic dynamics of the network could be predicted from the dynamics of a single node with appropriately scaled parameters of the input current. Thus, we first analyze the lower bound of the Fisher information of a single node. Equation (7) indicates that the Fisher information is determined by two factors: the term $\left(\frac{d\mu_{T}}{dI}\right)^{2}$, which is related to the slope of the so-called $f-I$ curve (mean firing rate vs input current curve) and determines the sensitivity of a neuron to small variations of the input current. The sensitivity is largest in the vicinity of the bifurcation point, where the limit cycle is born, and where the slope of the $f-I$ curve is the steepest. In this region, the Fisher information is high. However, the second factor in Eq. (7), the variance of the spike count, may degrade the Fisher information. In the excitable regime, when the input current is below its bifurcation value and APs are induced by noise, the phenomenon of stochastic resonance is observed [45], i.e. due to the competition of two factors, the sensitivity and the spike count variance, the Fisher information possesses a maximum at an optimal noise intensity [36].

Figure 8(a1) shows the lower bound of the Fisher information, $J_{\text{LB}}$, for a single HH node as a function of input current and noise intensity. The Fisher information is maximal for an input current $I$, which brings the system close to the transition to periodic spiking, i.e. 28 – 29 $\mu$A/cm2. In Fig. 8(a1), a vertical section across the map corresponds to the dependence of the Fisher information on noise intensity. As can be seen, such a dependence is non-monotonous in the excitable regimes, e.g. for $I=15$ or $I=20$ $\mu$A/cm2, indicating the phenomenon of stochastic resonance [45], reported before for the original Hodgkin-Huxley neuron model in Ref.[36]. Indeed, stochastic resonance is a generic phenomenon in excitable systems [45, 12] and so the Frankenhaeuser-Huxley (FH) nodal model demonstrates qualitatively similar parameter dependence, shown in Fig. 8(a2). In the absence of noise, a stable equilibrium of the single FH node passes through a subcritical Andronov-Hopf bifurcation at $I_{\text{AH}}\approx 60.19$ $\mu$A/cm2. Consequently, the Fisher information in Fig. 8(a2) is maximal around this value, similar to the HH node. In the excitable regime, e.g. $I=40$  $\mu$A/cm2, the Fisher information vs. noise intensity passes through a maximum, demonstrating stochastic resonance, again, qualitatively similar to the HH node.

The scaling relations for the input current, Eq. (22), and noise intensity, Eq. (23), enable us to predict the signal detection capability of a tree network in the regime of strong coupling. Given the branching, $d$ and the input current to the peripheral nodes, an increase in the number of generations (i.e. tree size) results in a decrease of the effective input current and noise. Then, depending on the particular values of $I$ and $D$, signal detection by the tree may show distinct dependencies on the tree size, $G$. This is illustrated in Fig. 8(a1,a2) by superimposing the scaling of the input current and noise intensity on the Fisher information map of the single node. In particular, our theory predicts that in the excitable regime, i.e. when the network does not produce sustained periodic firing in the absence of stochastic inputs, the scaling of $I$ and $D$ may bring an effective operating point of the network across the local maximum of the Fisher information. As can be seen, for instance, for the input current $I=30$ or $40$ $\mu$A/cm2 for HH nodal model, an increase of the number of generations to $G=$2–4 brings the effective operating point to regions of higher values of the Fisher information; further growth of the tree size eventually suppresses AP firing and thus small signals cannot be detected. In contrast, in the oscillatory regime (e.g. for $I=60$ $\mu$A/cm2 for HH nodal model), the increase of the network size moves the operation point always to regions of higher values of the Fisher information and so the discriminability increases monotonously with the tree size, $G$. Interestingly, one could predict the input current to the network which for a tree with large enough generations would result in an effective operating point close to bifurcation value of the single node. For example, for the HH nodes, such a value of the external current is $I=55$ $\mu$A/cm2 and for the FH node, $I=115$ $\mu$A/cm2. For such currents increasing the tree size should result in a higher degree of signal discrimination.

To test these predictions we computed the discriminability (6) for trees with different numbers of generations and for the single node with the scaled values of constant input current and noise intensity according to Eqs. (22, 23). Figure 8(b1,b2) shows excellent correspondence between the respective discriminabilities.

Conclusion

We have studied the emergence of noisy periodic spiking in regular tree networks of coupled excitable elements. Using biophysical models of excitable nodes, we showed that noisy periodic network spiking can be generated, although the periphery of the tree is excited by random and independent inputs (Fig. 4). The firing rate and coherence of spiking can be maximized by varying the coupling strength and is altered by changing the network topology (Figs. 5,6).

We put special emphasis on the strong coupling regime, which refers to the case of excitable nodes of Ranvier linked by myelinated (dendrite or axon) fibers of a neuron.

It is intuitively clear that in the strong coupling limit, the collective dynamics of the network could be described by a single effective excitable system. We have derived the corresponding scaling relations for random inputs Eqs. (22, 23) which allows for reliable predictions of the collective network response based on the stochastic dynamics of a single isolated node with scaled input parameters. Stochastic excitable systems demonstrate non-trivial behaviour versus the noise intensity. Examples include the phenomena of coherence resonance [46], whereby the variability of spiking events (e.g. coefficient of variation) is minimal for non-zero noise intensity, and stochastic resonance, characterized by non-monotonous dependence of a response to an external signal on the noise intensity [45]. Similar phenomena have been observed in networks of excitable elements. In particular, the phenomena of system size stochastic [47] and coherence resonance [14], which are also observed in strongly-coupled star networks of excitable elements [34]. As we have shown in the present paper, the phenomenon of system size stochastic resonance also occurs in strongly coupled tree networks, i.e. the number of generations in a tree network of excitable elements can be tuned in order to optimize the network ability to discriminate between different input signals. In particular, our analytical approach allows for the prediction of optimal tree sizes and branching ratio.

The analytical approach developed here can be extended to random trees [48] in which the branching ratio varies among different generations, yielding similar scaling relations in the strong coupling limit. While we considered networks of identical nodes, our approach can be readily extended to the inhomogeneous case, as long as the condition for strong coupling Eq. (16) is satisfied.

Our results suggest a mechanism for the emergence of noisy periodic firing and information coding by peripheral sensory neurons which possess branched tree-like myelinated dendrites [28]. Such neurons may possess multiple spike initiation zones at peripheral nodes (heminodes) and nodes of Ranvier at branching points. Examples of the muscle spindles [27] and cutaneous mechanoreceptors [31] indicate that myelinated dendritic trees extend to up to 7 generations. Myelin provides low-resistance links between nodes and fast saltatory conduction of APs, which corresponds to strong coupling between the nodes of Ranvier. For example, the average diameter of a cat muscle spindles afferents ranges from $3$ to $13$ $\mu$m, while links between nodes are relatively short, $50$ – $200$ $\mu$m [27]. An estimate of the coupling strength from Eq. (3) yields values well within the range of the strong coupling regime used in our study. The collective noisy periodic firing then may occur due to the synchronized noise-induced generation of APs by stimulating the peripheral heminodes, as described by our model. Given the biophysical properties of the nodes of Ranvier and the sensory inputs, the variability of interspike intervals and the stimulus discrimination capability of a neuron are determined by the ratio of the number of signal-receiving peripheral heminodes to the total number of nodes in the network.

Acknowledgements

We thank D.F. Russell, E. Schoell, T. Isele for fruitful discussions. AN acknowledges support by the Lobachevsky University of Nizhny Novgorod through the Russian Science Foundation grant 14-41-000440. LSG thanks Ohio University for hospitality and support.

Author contributions. AN formulated the problem. JK performed analytical calculation. AKN and AN performed numerical simulations. JK, AKN, LSG and AN wrote and reviewed the manuscript.

Competing financial interests. The authors declare no competing financial interests.

Supplementary material.

Emergent stochastic oscillations and signal detection in regular tree networks of strongly coupled excitable elements

Justus A. Kromer, Ali Khaledi-Nasab, Lutz Schimansky-Geier, Alexander B. Neiman

Models for Nodes of Ranvier

A Hodgkin-Huxley type model (HH) for a node of Ranvier contains only sodium and leak ionic currents. Thus, in eq.(2) of the main paper, the ionic current becomes $I_{\text{ion}}=I_{\text{Na}}+I_{\text{L}}$. For the sodium current we used the Hodgkin-Huxley (HH) type kinetics [1, 2], $I_{\text{Na}}=g_{\text{Na}}m^{3}h(V-V_{\text{Na}})$, where $g_{\text{Na}}=1100$ mS/cm2 is the maximal value of the sodium conductance and $V_{\text{Na}}=50$ mV is the Na reversal potential. The gating activation and inactivation variables obey the dynamics

[TABLE]

with the following rate functions:

[TABLE]

The leak current is $I_{\text{L}}=g_{\text{L}}(V_{k}-V_{\text{L}})$ with $g_{\text{L}}=20$ mS/cm2 and $V_{\text{L}}=-80$ mV.

The Frankenhaeuser-Huxley (FH) model [3] uses four ionic currents: sodium, potassium, persistent sodium and leak: $I_{\text{ion}}=I_{\text{Na}}+I_{\text{K}}+I_{\text{p}}+I_{\text{L}}$. The currents are given by

[TABLE]

where $E=-70+V$; $P_{X}$ are the permeabilities of sodium (Na), potassium (K) and persistent (p) ionic currents:

[TABLE]

where $\bar{P}_{X}$ are maximal values of permeabilities of corresponding channels. In Eq. (32) $[X]_{0}$ and $[X]_{i}$ are corresponding extracellular and intracellular ionic concentrations, of Na and K ions, respectively; for the persistent current, $X=p$, we have $[X]\equiv[\text{Na}]$. The gating variables follow,

[TABLE]

and the rate functions are:

[TABLE]

where $s$ is a scaling factor. The parameters for the model were taken from the original Frankenhaeuser-Huxley paper [3] with three modifications to reduce the frequency of periodic spiking when a sufficiently-high constant current is injected: (i) the rate equations for the gating variables included a scale factor $s=0.3$; (ii) the maximal permeability of potassium channels is reduced to $\bar{P}_{\text{K}}=3.6\times 10^{-4}$ cm/sec; and (iii) the leak conductance was reduced to $g_{\text{L}}=6$ mS/cm2. Other parameters are the same as in the original FH model: $[Na]_{0}=114.5$ mM, $[Na]_{i}=13.74$ mM, $[K]_{0}=2.5$ mM, $[K]_{i}=120$ mM; $\bar{P}_{\text{Na}}=8\times 10^{-3}$ cm/sec; $\bar{P}_{\text{p}}=54\times 10^{-5}$ cm/sec; $V_{\text{L}}=0.026$ mV. Constants $R$ and $F$ are the universal gas constant and the Faraday constant, and $T=293.15$ K.

Inverse of Tridiagonal Toeplitz Matrix

In order to calculate the effective current $I_{\text{mod}}$ eq. (20) and noise intensity $D_{\text{mod}}$ eq. (21), we need to evaluate the inverse matrix $\mathbf{B}^{-1}$ of $\mathbf{B}$. To this end, we apply results from Ref. [4] in which the components of the inverse of a tridiagonal Toeplitz matrix are given in terms of two sequences $\{v_{i}\}$ and $\{u_{i}\}$ with $i=1,2,...,G$. If we apply their results to the matrix $\mathbf{B}$, eq. (17), we obtain for its $(i,j)$th component

[TABLE]

The sequences $\{u_{i}\}$ and $\{v_{i}\}$ are only defined up to a multiplicative constant and it is convenient to set $u_{1}=1$. Then the sequences can be calculated from

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

Since we are interested in the dynamics of the first order node ($g=0$) when only peripherals are subject to noisy currents, we only need to evaluate a single component, $(1,G)$, of $\mathbf{B}^{-1}$. For this component we find from eq. (44), $\mathbf{B}^{-1}_{1,G}=v_{G}$. Applying the definition of $v_{k}$, eq. (Inverse of Tridiagonal Toeplitz Matrix), multiple times, we obtain

[TABLE]

Next, we evaluate the product for which one can show by using mathematical induction that

[TABLE]

Finally, multiplication by $d_{G}$ and applying eq. (Inverse of Tridiagonal Toeplitz Matrix) yields

[TABLE]

This yields for the $(1,G)$ component of $\mathbf{B}^{-1}$

[TABLE]

Using this in eqs. (19) yields the effective current and noise intensity, respectively.