Variability of collective dynamics in random tree networks of strongly-coupled stochastic excitable elements

We assume that all nodes and passive links are identical, except for the leaf nodes, representing heminodes, which receive external inputs. Given a tree with $\mathcal{N}$ nodes, the dynamics of the nodes’ membrane potentials are approximated by a discrete cable model Ermentrout and Terman (2010) in which the membrane potential of the $n$th node is governed by

Here the index $n=1,2,...,\mathcal{N}$ marks the respective node. In particular, $n=1$ refers to the root node. In Eq. (1) $C$ is the nodal membrane capacitance and the term $I_{\text{ion}}$ represents the nodal ionic currents. In the following, we use a Hodgkin-Huxley type (HH) model for the ionic currents of nodes of Ranvier, which is a simplified version of the model used in McIntyre, Richardson, and Grill (2002). Ionic currents are represented by Na⁺ and leak currents, $I_{\text{ion}}=I_{\text{Na}}+I_{\text{L}}$, Kromer, Schimansky-Geier, and Neiman (2016); Kromer et al. (2017). The Na⁺ current is $I_{\text{Na}}=g_{\text{Na}}m^{3}h(V-V_{\text{Na}})$, where $g_{\text{Na}}=1100$ mS/cm² 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

with the following rate functions:

The leak current is given by $I_{\text{L}}=g_{\text{L}}(V-V_{\text{L}})$ with $g_{\text{L}}=20$ mS/cm², $V_{\text{L}}=-80$ mV, and the nodal capacitance is set to $C=2$ $\mu$F/cm².

In Eq. (1) the coupling between nodes is described by $\kappa\sum_{j=1}^{N}A_{n,j}(V_{j}-V_{n})$, where $\mathbf{A}$ is the adjacency matrix of the undirected rooted tree graph. It is a $\mathcal{N}\times\mathcal{N}$ symmetric matrix with elements $A_{i,j}=1$ for connected nodes $i$ and $j$, and $A_{i,j}=0$ for unconnected nodes, see Appendix A for more details. In the following the coupling strength, $\kappa$, is used as a control parameter. However, the physiologically-relevant values of $\kappa$ can be estimated from the sizes of a node, the myelinated links, and the axoplasmic resistivity Ermentrout and Terman (2010), giving the range of $\approx$ 125 – 1500 mS/cm² Kromer, Schimansky-Geier, and Neiman (2016); Kromer et al. (2017). The external currents $\mathcal{J}_{n}$ are applied to the leaves only and consist of a constant and a noisy part, i.e

where $l$ denotes indices of leaf nodes; $\delta_{n,l}$ is the Kronecker delta. The zero-mean Gaussian white noise $\xi_{l}(t)$ with intensity $S$ is uncorrelated for different leaves, i.e. $\langle\xi_{i}(t)\rangle=0$, $\langle\xi_{i}(t)\xi_{j}(t+\tau)\rangle=\delta_{i,j}\,\delta(\tau)$. Thus, leaf nodes receive random uncorrelated inputs. Other sources of noise, e.g. due to fluctuations of nodal ion channel conductances are neglected. In contrast to regular tree networks where inputs are administered only to nodes in the last generation of a tree Kromer et al. (2017), leaves can occur at any generation in random tree networks, see Fig. 1.

Eqs.(1–6) were integrated numerically using the explicit Euler - Maruyama method with time step of $0.1$ $\mu$s for 60 – 600  seconds long simulation runs.

We first discuss repetitive action potential generation at the root node for deterministic input currents, $S=0$ in Eq. (6). Qualitatively different dynamical operation modes of the root node are separated by a threshold value, $J_{\text{th}}$, of the constant current, $J$, applied to leaf nodes. While APs evoked at the leaf nodes do not fire up the root node for low currents, values of $J$ exceeding $J_{\text{th}}$ result into sustained periodic firing of the root node. This is reminiscent of the dynamic behavior of a single isolated node. The latter is at resting state in the absence of external input, $\mathcal{J}=0$. A sufficiently high constant current results in an Andronov-Hopf bifurcation of the equilibrium state rendering an isolated node to fire a periodic sequence of action potentials (APs). For a single node this Andronov-Hopf bifurcation occurs at $J_{\text{AH}}\approx 29.06~\mu$A/cm² Kromer et al. (2017).

As in Ref. Kromer et al. (2017), we numerically calculated the threshold current $J_{\text{th}}$, which is the minimum constant current applied to the leaf nodes, which for a given coupling strength, results in the repetitive sustained generation of full-size APs (a voltage spike of at least 60 mV magnitude) at the root node. As for regular trees Kromer et al. (2017), the threshold current depends on the coupling strength, $J_{\text{th}}(\kappa)$, as exemplified in Fig. 2b.

As the coupling strength increases, more input current is required to sustain firing of the root node. In consequence, the threshold current increases with $\kappa$ and, finally, saturates at the limiting value $J_{\infty}:=\lim_{\kappa\to\infty}J_{\text{th}}(\kappa)$ for strong coupling. The strong coupling regime spans the range of physiologically realistic values, $\kappa>100$ mS/cm², for branched myelinated terminals of sensory neurons Kromer, Schimansky-Geier, and Neiman (2016); Kromer et al. (2017).

For a given value of $\kappa$ the root node shows a sustained sequence of APs, if values of the input current are above $J_{\text{th}}(\kappa)$, shown in Fig. 2b, and no APs if the value of $\mathcal{J}_{n}$ is below that curve. These two regimes are referred to as oscillatory and excitable, respectively, in the following. Note that for weak coupling the dynamics of the tree can be quite complex, e.g. not every AP generated at leaf nodes may propagate all the way to the root node. Such regimes will be studied elsewhere.

For strong enough coupling and sufficient input current, all nodes in the tree are synchronized and fire periodic sequences of APs. As in the case of regular trees, the limiting value of the threshold current for strong coupling, $J_{\infty}$, matches the threshold value of the current of an effective single node with parameters re-scaled by the ratio of the number of inputs (leaf nodes) and the total number of nodes as, $J_{\infty}=(\mathcal{N}/\mathcal{H})J_{\text{AH}}$, see dashed line in Fig. 2b. This value is derived in Sec. III.2 and Appendix C.

In the presence of noise, the AP generation becomes stochastic. We are particularly interested in statistical properties of spike trains generated at the root node. We extracted a sequence of spike times of the root node, $\{t_{i}\}$, from $60-600$ s long simulation runs. The sequence of interspike intervals (ISIs), $\tau_{i}=t_{i+1}-t_{i}$, is characterized by the firing rate, $r$, and by the coefficient of variation (CV), $\mathcal{C}_{\tau}$,

where $\overline{\tau}$ and $\sigma_{\tau}$ are the mean and the SD of the sequence of ISIs, respectively. The bar stands for the averaging over all ISIs in the spike train of the root node. The CV quantifies the ISI variability.

Figure 3 exemplifies the stochastic firing dynamics for the tree shown in Fig. 2a. The constant current was chosen such that the tree is in the excitable regime for $\kappa>20$ mS/cm². As in regular trees Kromer et al. (2017), the firing rate, $r(\kappa)$, depends non-monotonously on the coupling strength. For weak coupling, spikes generated at the leaves often fail to propagate to the root, leading to its sparse firing, Fig. 3a. This results in low values of the firing rate and large values of the CV. Furthermore, the root and leaf nodes fire asynchronously. Increasing the coupling strength leads to stronger interaction between nodes. In consequence, synchronous coherent firing with large firing rates and small CVs emerges, see middle panel of Fig. 3a. As the coupling increases, more input current is required to sustain the network firing, see Fig. 2b. Consequently, for constant input current the firing rate decrease again. Competition of these two tendencies results in a maximum of the root node firing rate as a function of the coupling strength. For strong coupling, nodal firing is perfectly synchronized, see rightmost traces in Fig. 3a, and the firing rate and the CV saturate at their limiting values, indicated by dashed lines in Fig. 3b,c. In Sec. III.2 we show that in the case of strong coupling, a tree can be well represented by its root node which receives effective input with rescaled constant current and noise intensity.

An ensemble of random trees can be constructed as a set of realizations of a stochastic branching process Drmota (2009). In this paradigm, the tree shown in Fig. 2a is just one possible realization of a random tree network. Each network realization causes certain statistical properties of its root node’s firing characteristics, such as those depicted in Fig. 3.

In each realization, branching starts at the root and continues up to a prescribed maximal number of generations, $G$, or until all branches end in leaf nodes. Each node in a certain generation $g$ is located at the same path distance, $g$, from the root. The latter defines the generation $g=0$. Furthermore, each node in the $g$-th generation, $g<G$, is a parent to a random number of offspring in the $(g+1)$-st generation. However, each child node has only one parent.

The $k$-th realization of a tree network consists of several generations each containing a random number of nodes, $\mathcal{D}_{g,k}$, which we model by using the Galton-Watson process Harris (2002),

Here $g$ indicates the generation and $k=1,2,\ldots$ indicates a particular realization of a tree network. The sum runs over all nodes in generation $g$. The number of offspring, $d_{g,i,k}$, of a certain parent node $i$ in generation $g$ is an independent random variable generated from the probability mass function (PMF), $p_{g}(d)$. In the following sections, we will consider several examples of branching PMFs.

We consider basic properties of a single tree network realization, first. The total number of nodes $\mathcal{N}_{k}$ of a tree network realization $k$ is a random variable and obtained by summing $\mathcal{D}_{g,k}$ over all generations,

where $\mathcal{G}_{k}\leq G$ is the height of the tree and is given by the actual number of generations in the current tree network realization. Thus, there might be no nodes in the outer generations for particular realizations of the Galton-Watson process, Eq. (8). The number of leaf nodes in a particular generation $g$ is also a random variable given by

By construction, branching terminates at the maximum generation $G$ and all nodes in that generation are leaf nodes. In general, however, leaf nodes can be found in any but the $0$-th generation. Thus, a node $i$ in generation $0<g\leq G$ becomes a leaf node with probability $p_{g}(0)$, i.e. if $d_{g,i,k}=0$. The total number of leaf nodes, $\mathcal{H}_{k}$, is a random variable, obtained by summing $h_{g,k}$ over all generations, i.e.

In the special case of no extinction, i.e. if $p_{g}(0)=0$ for $0\leq g<G$, all leaf nodes are in the peripheral generation, $g=\mathcal{G}_{k}=G$ and $\mathcal{H}_{k}=\mathcal{D}_{G,k}$.

Next, given the branching PMF, $p_{g}(d)$, and the maximum allowed number of generations, $G$, some statistical properties of random trees formed by the Galton-Watson process, such as the mean number of total nodes and leaves, their variances, etc. can be determined by the standard method of probability generating functions Harris (2002); Drmota (2009); Newman (2010). In the present paper, we are mainly interested in the influences of structural variability, as arising from different tree network realizations for the same branching PMF, on the firing statistics of the root node. We characterize that statistics for a single tree network realization $k$ by calculating the firing rate, $r_{k}(J,S,\kappa)$, and the CV, $\mathcal{C}_{k}(J,S,\kappa)$ according to Eqs. (7). Note that these quantities depend not only on the network realization with the particular coupling structure, but also on input current, i.e. $J$ and $S$, and the coupling strength $\kappa$.

The variability of a certain quantity $Q_{k}$, of the $k$-th tree realization, can be assessed by calculating its ensemble average mean and standard deviation (SD):

This yields the ensemble averaged firing rate $\langle r(J,S,\kappa)\rangle$ and its normalized standard deviation:

The latter provides a measure of variability for the root nodes’ firing rates resulting from tree network realizations from the same branching PMF with respective input currents and coupling strengths.

By defining sets of identical or isomorphic trees in an ensemble, averaging can be performed using the probability distribution of sets of identical trees. Denoting the set of parameters, which uniquely defines trees by $\{T_{k}\}$, the $k$-th realization of the quantity $Q$ is denoted as $Q_{k}:=Q(\{T_{k}\})$. Its average can be formally written as, $\langle Q\rangle=\sum_{\{T\}}Q(\{T\})\,P(\{T\})$, where $P(\{T\})$ is the PMF of non-identical trees and the summation run over all possible values of the parameters set, $\{T\}$. Various levels of coarse-graining methods can be used to simplify the ensemble averaging.

In the present paper we restrict to identical nodes and interconnections except for the external input, which is only applied to leaf nodes. As a first way of coarse-graining, we consider trees with the same number of nodes and leaf nodes in each generation as identical. We get $(2G-1)$-tuple $\{T_{k}\}=(\mathcal{D}_{1,k},...,\mathcal{D}_{G,k},h_{1,k},...,h_{G-1,k})=(\{D_{k}\},\{h_{k}\})$. We note that trees with identical tuples $(\{D_{k}\},\{h_{k}\})$ may still possess different connectivities. The ensemble of trees is then characterized by the joint $(2G-1)$-dimensional PMF of number of the nodes and leaves in all generations, $P_{2G-1}(\mathcal{D}_{1},...,\mathcal{D}_{G};h_{1},...,h_{G-1})\equiv P_{2G-1}(\{\mathcal{D}\},\{h\})$. The tree-ensemble averages, Eqs. (12), are then approximated by

where the summations run over all possible values of $\mathcal{D}_{1},...,\mathcal{D}_{G}$ and $h_{1},...,h_{G-1}$. As a second way of coarse-graining, we consider trees with the same total number of nodes and leaf nodes as identical. As we will show in Sec. III.2, this simplification yields sufficient results in the strong-coupling limit $\kappa\rightarrow\infty$. We parameterize a particular realization of tree network by the tuples $(\mathcal{H}_{k},\mathcal{N}_{k})$ and then carry out averaging similar to Eq. (14), but with 2-dimensional PMF of the total number of leaves and nodes, $P_{2}(\mathcal{H},\mathcal{N})$.

Throughout the paper we focus on small trees, $2<G\leq 4$, with branching supported on a bounded interval. This is consistent with the topology of branched myelinated terminals of sensory neurons Banks, Barker, and Stacey (1982); Banks et al. (1997); Lesniak et al. (2014); Marshall et al. (2016). For such trees the number of configurations with distinct tuples $(\{\mathcal{D}\},\{h\})$ in the same ensemble, is rather small. In particular, for binary trees this enables us to list all non-identical trees in the ensemble and calculate the corresponding joint PMF. Furthermore, dynamical measures, such as the firing rate and CV, can be calculated numerically for the complete small set of trees. Thus, the structure-induced variability of these measures can be calculated according to Eq. (14).

In binary trees each node has at most two offspring. Here we consider two types of binary trees: full and non-full binary trees. The so-called full binary tree, is a tree in which every internal node has two offspring and leaves have none. In contrast, in non-full binary trees, which we term as general binary trees, the number of offspring of any internal node can be either one or two.

In the physiologically-relevant case of strong coupling, the stochastic firing of all nodes is synchronized. In the synchronized state, the dynamics of individual tree network realizations, Eqs. (1–6), can be approximated by that of an effective single node (see Appendix C):

Here $V_{k}(t)$ is the effective membrane potential and $\mathcal{J}_{\text{eff},k}(t)$ is the effective stochastic input current. The index $k$ refers to the actual tree realization. The latter encodes the tree structure. The nodal ionic current, $I_{\text{ion}}[V_{k},m_{k},h_{k}]$, and the gating variables, $m_{k}$ and $h_{k}$, are given by the same equations as for the network model in Sec. II.1. Equations for $\mathcal{J}_{\text{eff},k}(t)$ for regular trees were derived in Kromer et al. (2017). In Appendix C, we extend their approach to random trees and obtain,

Here $J$ and $S$ are the constant current and the noise intensity, respectively. Both specify the strengths of the noisy input current applied to the leaves in the actual tree realization; $\xi(t)$ is Gaussian white noise.

In the following, we use Eqs. (19) and (III.2) as an approximation for the dynamics of the root node, i.e. $v_{k}(t)\approx V_{1,k}(t)$, in the strong coupling limit. Thus, we replace the ensemble of trees by an ensemble of their effective root nodes. The respective total numbers of nodes and leaves, $(\mathcal{H}_{k},\mathcal{N}_{k})$, in the individual tree realization determine the effective current, $\mathcal{R}_{k}J$, and noise intensity, $S\mathcal{S}_{k}$, in Eq. (III.2).

The randomness of tree ensembles may lead to qualitatively different dynamics of individual network realizations. In the present section we consider its consequence for the threshold current setting the onset of repetitive spiking of the tree’s root node in the case of deterministic input currents, i.e $S=0$.

We numerically calculate the threshold current as a function of coupling strength for each of $25$ distinct full binary trees depicted in Fig. 4. This results in one curve for each tree network realization, each one similar to the curve shown in Fig. 2b. Besides its dependence on the coupling strength $\kappa$, the threshold current for a single tree network realization, $J_{\text{th},k}$, is also a function of the number of nodes in the $3$-rd and $4$-th generation, i.e. $J_{\text{th},k}=J_{\text{th}}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})$. It can be expressed in units of $J_{\text{AH}}$ by introducing the dimensionless scaling factor $\mathcal{R}^{-1}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})$, i.e. $J_{\text{th},k}=J_{\text{th}}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})=\mathcal{R}^{-1}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})J_{\text{AH}}$. At $J=J_{\text{AH}}$ a single isolated node enters the repetitive spiking regime by undergoing an Andronov Hopf bifurcation. We will refer to $\mathcal{R}^{-1}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})$ as the normalized threshold current in the following. For strong coupling the threshold current of a single tree network realization approaches its limiting value, $J_{\infty,k}(\mathcal{H},\mathcal{N})$, given by Eq. (21). The latter depends only on the total number of nodes and leaf nodes, i.e. $J_{\infty}(\mathcal{H},\mathcal{N})=\mathcal{R}^{-1}_{\infty}(\mathcal{H}_{k},\mathcal{N}_{k})J_{\text{AH}}=(\mathcal{N}_{k}/\mathcal{H}_{k})J_{\text{AH}}$. Here $\lim_{\kappa\to\infty}\mathcal{R}(\kappa,\mathcal{D}_{3,k},\mathcal{D}_{4,k})=\mathcal{R}_{\infty}(\mathcal{H}_{k},\mathcal{N}_{k})=\mathcal{H}_{k}/\mathcal{N}_{k}$.

Normalized threshold currents are shown in Fig. 8(a,b). We find that it increases monotonically with the coupling strengths and finally saturates for strong coupling. For weak coupling, the actual tree structure matters as it determines the paths action potentials have to travel in order to excite the root node. Here trees with various configurations, but identical numbers of leaves and nodes, fall into tight clusters around distinct values of $\mathcal{R}^{-1}$. With the increase of coupling the clusters blur. For general binary trees, Fig. 8(b), the range of threshold currents is significantly larger than for full binary trees, Fig. 8(a). Finally, for strong coupling curves for individual sample trees saturate. Their limiting values are well approximated by the theoretically predicted strong coupling limit $J_{\infty,k}=(\mathcal{N}_{k}/\mathcal{H}_{k})J_{\text{AH}}$, as illustrated in Fig. 8(c,d).

These results indicate that for weak and moderate coupling the onset of spike generation may be strongly affected by the particular tree structure. However, for physiologically relevant coupling strengths, $\kappa>100$ mS/cm², threshold currents are close to their strong-coupling limits and their values mainly depend on the statistics of the total number of nodes and leaf nodes.

For stochastic input currents, i.e. $S>0$, variability in tree structures interacts with variability caused by stochastic inputs. In order to study the stochastic dynamics we prepare ensembles of excitable trees, i.e. no spike generation at the root node if the input current had no stochastic component. To this end, we set the value of the constant input current such that it is below the sample trees’ threshold currents for strong coupling, but causes non-vanishing firing rates, $>2$ Hz, of all possible tree realizations. For the two types of binary trees, used in the previous section, this is achieved by using a constant input current of $J=38.5~\mu$A/cm². This corresponds to the normalized threshold current, $\mathcal{R}^{-1}=1.325$, shown by the dashed lines in Fig. 8a,b. For the coupling strengths, $\kappa>30$ mS/cm², all curves of the threshold current in Fig. 8a,b lie above the dashed lines, indicating that all trees are indeed excitable.

Then, we simulated Eqs.(1–6) for all possible non-identical tree network realizations and estimated the firing rate and CV of their root nodes, $r(\kappa,\{\mathcal{D}\},\{h\})$ and $\mathcal{C}_{\tau}(\kappa,\{\mathcal{D}\},\{h\})$, respectively.

The obtained measures of spike train statistics are shown in Fig. 9. For individual tree realizations, theses measures show qualitatively similar a dependences on the coupling strength. In more detail, we find that firing rates become maximal at intermediate coupling strengths for both full and general binary trees. CVs attain their minimum values at similar coupling strengths. This indicates most regular firing of the root node at intermediate coupling strengths. Note that this is in qualitative agreement with previous results on regular tree networks presented in Kromer et al. (2017).

Considering the spike train statistics of the entire tree ensemble, we find that structural variability hardly affects firing rates and CVs for weak coupling. In contrast, firing rates and CVs of individual tree realizations strongly differ for physiologically relevant strong coupling. For such coupling, the spike train statistics of individual tree realizations are well-approximated by those of single isolated nodes with effective currents and noise intensities given by Eq. (III.2). This is further illustrated in Fig. 9(c,f), where the ISIs statistics of root nodes of tree realizations is compared to those of effective isolated nodes according to the strong-coupling approximation.

Importantly, the strong-coupling approximation also allows for the prediction of lower and upper bounds of the firing rates and CVs, shown by the dashed lines in Fig. 9a–d. The scaling parameters $\mathcal{R}_{\infty,k}$ and $\mathcal{S}_{\infty,k}$ of the effective current in Eq. (III.2) define the range of the trees’ firing rates and CVs at strong coupling, see Fig. 7. In that sense, they provide bounds for the influence of structural variability on the spike train statistics of the root node. Since all tree realizations are excitable for strong coupling, the highest values of scaling parameters $\mathcal{R}_{\infty,k}$ yield the highest firing rate and lowest CV. In contrast, small values of $\mathcal{R}_{\infty,k}$ yield low effective currents, Eq.(III.2), and drive the network deep into the excitable regime. This causes Poisson-like statistics of spike generation resulting in low firing rates and CVs close to one.

Next we consider the ensemble-averaged ISI statistics. It is obtained by averaging the firing rates and CVs according to Eq.(14). Averaging is simplified for the full binary trees, as the joint PMF, Eq. (16), for $G=4$ depends only on the numbers of nodes in the $3$-rd and $4$-th generations. For the general binary trees with $G=3$ the ensemble averaging is performed using the $5$-dimensional joint PMF, $P_{5}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3};h_{1},h_{2})$ which we estimated numerically. Ensemble-averaged measures of spike train statistics are shown in Fig. 10 for different values of the zero-branching probability $p_{0}$. The ensemble-averaged firing rates and CVs follow curves that are qualitatively similar to those for individual trees. However, $p_{0}$ strongly affects the ensemble-averaged statistics in the strong coupling regime. Low probabilities cause on average taller sample trees with smaller fractions of leaf nodes. In consequence, the effective current and noise intensity in Eq. (III.2) become smaller, which drives the network deep into the excitable regime, and results in low firing rates and large CVs. Larger values of $p_{0}$, refer to an increased fraction of short trees with larger fractions of leaves, which receive inputs. The latter results in higher firing rates and smaller CVs.

In order to quantify the influence of structural variability on the firing statistics of root nodes, we consider the normalized standard deviation, $\mathcal{C}_{r}$, of the distribution of firing rates, see Eq. (13). We find that firing rate distributions of general binary trees show larger variability than those for full binary trees. Besides the previously noticed fact that the structural variability is most pronounced for strong coupling, we find a maximum of the $\mathcal{C}_{r}$ for a finite value of zero-branching probability. For full binary trees the number of possible distinct trees approaches one for $p_{0}\rightarrow 0$, i.e. only the regular tree with branching two, and $p_{0}\rightarrow 1$, i.e. only the regular tree with the minimal number of generations. The pronounced maximum of $\mathcal{C}_{r}$ expresses the trade-off between trees becoming more regular as $p_{0}$ approaches either one of those limits. In general binary trees, however, the limit $p_{0}\rightarrow 0$ still results in an ensemble of 17 possible trees with distinct pairs of number of nodes and leaf nodes, as all combinations of branching one and two are possible. Consequently, the normalized SD at $p_{0}=0$ remains finite.

The variability of the root node’s firing rates indeed depends on the input current to the leaves. For strong coupling a particular coupling structure is imprinted in the scaling of input current according to Eq.(III.2), and so the spread of effective input currents for trees in the ensemble translates into the spread of their firing rates. This can be seen in Fig.7 by comparing locations of tree realizations for two values of input currents to leaves. For $J=38.5~\mu$A/cm², all tree realizations are in the excitable regime and their positions in the heat map cut across a wide range of firing rates. For $J=70~\mu$A/cm², however, most tree realizations are in oscillatory regime, cutting across a narrower range of firing rates. A decrease of the input current shifts trees to the left in Fig. 7, i.e. deeper to the excitable regime. The increase of $J$ moves trees to the right and trees become oscillatory. As a result, the ensemble variability of firing rate is high for small input currents and low for large currents. These results are summarized in Fig. 12, which shows the dependence of the normalized SD of firing rate, $\mathcal{C}_{r}$, on the input current and the zero-branching probability for the strong-coupling limit. Except for small input currents, $J<30~\mu$A/cm², the firing rate variability of general binary trees is larger than that of full binary trees, as they cut across wider ranges of scaling factors of input current in Eq. (III.2).

We finally consider an example of random non-binary trees. In this example, random branching is drawn from a uniform distribution with at most $4$ offspring for each node, and trees with at most $G=4$ generations are allowed. To avoid short trees we set the zero branching probability to zero, $p_{g}(d)=0$, for generations $g=0,1,2$, as in the example of full binary trees. Thus, the branching PMF is given by

The PMF in Eq. (23) yields a large number of non-identical trees, i.e. trees that differ in the their numbers of nodes and leaf nodes per generation. Figure 13a exemplifies $3$ non-identical tree network realizations obtained from the PMF. Instead of counting distinct trees and calculating their corresponding probabilities, we analyzed tree network ensembles obtained from the PM using a brute-force approach. We generated an ensemble of $2000$ tree network realizations and then proceeded with the analyses of collective dynamics of coupled HH nodes as in previous sections.

For deterministic inputs, the threshold current as a function of coupling strength possesses a similar shape as those for binary trees (data not shown). This included the strong coupling regime, where threshold currents approached the value predicted by the strong coupling theory, $J_{\infty,k}=(\mathcal{N}_{k}/\mathcal{H}_{k})J_{\text{AH}}$.

We then set the constant input current at $J=38.5$ $\mu$A/cm², for which all network realizations resulted in excitable trees for coupling $\kappa>100$ mS/cm². Firing statistics of root nodes of the $2000$ tree realizations showed qualitatively similar dependences on the coupling strength as for binary trees of Fig. 9 (not shown). As for binary trees, structural variability of the firing rate of non-binary trees is small for weak and intermediate coupling and large for strong coupling. This is illustrated in Fig. 13c, by the probability distribution of the firing rate of root nodes. The latter broadens as the coupling strength increases.

As for binary trees, the firing statistics of non-binary tree network realizations can be predicted using the strong-coupling theory of Sec.III.2. In order to compare the firing rate statistics obtained from simulations for the $2000$ tree network realizations with that resulting from the strong coupling approximation, we first calculate the pairs of effective current and noise intensity for each tree realization using Eq. (III.2). The locations of those pairs in the current noise intensity space is shown in Fig. 13b. The figure also shows a heat map of the firing rate of a single node $\tilde{r}$ as a function of current and noise intensity. Then, using Eq. (III.2.1) we obtain an approximate of the firing rate for each pair of effective currents, i.e. for each tree network realization. The resulting distribution of firing rate approximates (single) is compared to the root nodes’ firing rates for different coupling strengths as obtained from simulations of network activity in Fig. 13c.