Inverse stochastic resonance in adaptive small-world neural networks

Abstract

Inverse stochastic resonance (ISR) is a counterintuitive phenomenon where noise reduces rather than increases the firing rate of a neuron, sometimes even leading to complete quiescence. ISR was first experimentally verified with cerebellar Purkinje neurons [A. Buchin et al., PLOS Computational Biology 12, e1005000 (2016)]. These experiments showed that ISR enables a locally optimal information transfer between the input and output spike train of neurons. Subsequent studies have further demonstrated the efficiency of information processing and transfer in neural networks with small-world network topology. We have conducted a numerical investigation into the impact of adaptativity on ISR in a small-world network of noisy FitzHugh-Nagumo (FHN) neurons, operating in a bistable regime with a stable fixed point and a limit cycle — a prerequisite for the emergence of ISR. Our results show that the degree of ISR is highly dependent on the value of the FHN model’s timescale separation parameter $\epsilon$. The network structure undergoes dynamic adaptation via mechanisms of either spike-time-dependent plasticity (STDP) with potentiation-/depression-domination parameter $P$, or homeostatic structural plasticity (HSP) with rewiring frequency $F$. We demonstrate that both STDP and HSP amplify the effect of ISR when $\epsilon$ lies within the bistability region of FHN neurons. Specifically, at larger values of $\epsilon$ within the bistability regime, higher rewiring frequencies $F$ are observed to enhance ISR at intermediate (weak) synaptic noise intensities, while values of $P$ consistent with depression-domination (potentiation-domination) consistently enhance (deteriorate) ISR. Moreover, although STDP and HSP control parameters may jointly enhance ISR, $P$ has a greater impact on improving ISR compared to $F$. Our findings inform future ISR enhancement strategies in noisy artificial neural circuits, aiming to optimize local information transfer between input and output spike trains in neuromorphic systems, and prompt venues for experiments in neural networks.

The impact of noise on nonlinear dynamical systems often yields counter-intuitive behaviors, such as the stabilization of otherwise unstable deterministic states, and the inhibition or enhancement of oscillations. Classic examples of stochastic enhancement include phenomena like stochastic resonance (SR) and coherence resonance (CR). In SR, adding (an optimal intensity of) noise to a nonlinear bi-stable system (or, systems with sensory thresholds, such as neurons) enhances its response to weak external periodic signals, making imperceptible signals detectable Longtin (1993). Experimental research has shown that SR maximizes information flow in sensory neurons at an optimal noise level Nozaki et al. (1999). CR, on the other hand, does not require an external periodic signal Pikovsky and Kurths (1997) and occurs when the regularity of noise-induced oscillations in an excitable system is a peaked (non-monotonic) function of the noise amplitude, with oscillators being optimally correlated at a certain non-zero noise intensity. CR has also been studied experimentally Pisarchik and Hramov (2023). More recent studies of stochastic individual neurons and neural network models have identified another form of non-linear response to noise, inverse stochastic resonance (ISR) Gutkin, Jost, and Tuckwell (2009); Tuckwell, Jost, and Gutkin (2009), an effect which later also has been observed experimentally in neurons Buchin et al. (2016). ISR leads to an inverse resonance, i.e., the response curve possesses a trough (minimum), unlike SR and CR which result in a response curve with a peak (maximum); or even the complete silencing of the spiking activity. While the effect of network adaptivity on SR and CR has been investigated extensively, research on the effect of adaptivity inverse stochastic (ISR) in single neuronal models and large-scale networks is lacking and the effects are unknown. The present study fills this gap and examines ISR in adaptive small-world neural networks, with two key mechanisms of adaptivity, spike-timing-dependent plasticity and homeostatic structural plasticity. Our study shows that adaptive mechanisms strengthen ISR, thus offering insights into the role and control of ISR in neural and other complex systems, and offers guidance for future experimental investigations.

I Introduction

Noise-induced resonances are counter-intuitive phenomena where the introduction of stochastic fluctuations, or noise, into a nonlinear dynamical system modulates the system’s activity and response. Stochastic resonance (SR) is a well-documented phenomenon where noise enhances a nonlinear system’s response to weak periodic signals and has been extensively studied and leveraged across various fields and particularly in neuroscience Hänggi (2002); Longtin (1993); Vázquez-Rodríguez et al. (2017); Gluckman et al. (1998); Bulsara et al. (1991). For coherence resonance, the addition of a specific amount of noise in excitable system renders oscillatory responses most coherent Pikovsky and Kurths (1997), and thereby optimizes signal detection and processing. Inverse stochastic resonance (ISR) is a counterpart phenomenon, since stochastic fluctuations or noise results in the reduction or even suppression of the system’s activity. Thus, ISR conceptually represents a paradigm shift in our understanding of noise-induced behaviors in complex systems, as it displays that noise cannot only play an excitatory but also an inhibitory role Tuckwell, Jost, and Gutkin (2009).

ISR was first identified in the context of neuronal dynamics with the observation that certain neurons, subjected to an optimal level of noise, could decrease and sometimes even completely quench their mean firing rate into quiescence Gutkin, Jost, and Tuckwell (2009). The experimental validation of ISR in both biological Buchin et al. (2016) and physical Huh, Shiomi, and Miyagawa (2023) systems underscores its relevance in the real world. By elucidating how noise can modulate neural activity, ISR provides insights into the delicate balance between excitation and inhibition in the brain. The discovery of ISR has significant implications for understanding both the function and dysfunction of the brain, especially in terms of regulatory mechanisms of neural circuits. For example, it is generally believed that neurons convey information via spiking interactions. Consequently, the occurrence of ISR can on one hand be viewed as a constraint on information processing in neural systems. On the other hand, ISR could also be crucial for computational processes that require diminished firing activity without chemical inhibitory neuro-modulation, or for those processes requiring intermittent bursts of activity Paydarfar, Forger, and Clay (2006). The presence of ISR could be beneficial in these specific scenarios.

Moreover, the quenching effect of ISR holds therapeutic potential, where controlled noise might be used to reduce excessive neural activity and prevent conditions marked by hyper-excitability, such as epilepsy Buchin (2015). Beyond neuroscience, the intriguing ISR effect is also relevant to various fields, including physical systems Huh, Shiomi, and Miyagawa (2023), ecological systems Touboul, Staver, and Levin (2018), particularly neuromorphic engineering.

From a dynamical system perspective, ISR typically arises in bistable and monostable stochastic nonlinear dynamical systems, where the underlying mechanisms are identified as noise-induced biased switching between periodic (limit cycles) and stationary (fixed points) attractors of the deterministic dynamics and noise-enhanced stability, respectively Bačić and Franović (2020); Zhu (2021). At an intermediate noise level, the switching rates become notably asymmetric, with the system spending substantially more time in a quasi-stationary state. This results in a distinctive non-monotonic relationship between mean-firing rate and noise, a hallmark of ISR.

Recent years have shown a surge of interest in exploring various aspects of the ISR phenomenon. One investigation Guo (2011) delved into the influence of temporal noise correlations on ISR, revealing that colored noise exerts a more significant suppressive impact on neural activity when compared to Gaussian white noise. Another study Tuckwell and Jost (2011) scrutinized ISR in a more realistic setting including spatial extension, demonstrating that mild noise can impede spiking when signal and noise inputs coincide spatially on the neuron; vice versa, if the signal and noise are unevenly distributed, the noise does not disrupt spiking activity, irrespective of the neuron’s extension on a spatial domain. A further study Uzuntarla, Barreto, and Torres (2017) found that ISR can emerge in static networks as a consequence of a variety of factors, including channel noise, connection strength, synaptic currents with excitatory and inhibitory terms, and topological features of the network, including degree distribution and mean connectivity degree.

A significant number of studies investigated ISR in individual neurons and neural networks as a function of various types of noise Wang et al. (2022); Lu et al. (2020); Zhao and Li (2019), spatial extension of the neuron model Tuckwell and Jost (2011); Liu et al. (2024); Zhang, Li, and Xing (2021), electromagnetic induction of due to ions moving in the neuronal membrane Ye, Yang, and Jia (2023), conductance-driven input Tuckwell, Jost, and Gutkin (2009), neuronal morphology Liu et al. (2024), time delays and coupling strength Liu et al. (2024); Zhang, Li, and Xing (2021), electrical synapses versus (inhibitory and excitatory) chemical synapses Uzuntarla, Barreto, and Torres (2017), electrical and chemical autapses (i.e., time-delayed synaptic connections where a neuron forms a synapse with itself) Zhang, Li, and Xing (2021), the average degree of scale-free neural network size Uzuntarla, Barreto, and Torres (2017). Despite the large range of situations explored in these research efforts, the question of how the dynamics in adaptive neural networks affect ISR remains underexplored. As we show in the present study, adaptive dynamics in a network due to dynamic plasticity have a profound effect on ISR.

One study Uzuntarla et al. (2017) examined ISR by analyzing how a single neuron with synaptic dynamics is affected by the independent (uncorrelated) spiking activity of numerous other neurons. In this setup, presynaptic neurons are modeled as independent Poisson spike generators emitting uncorrelated spikes with a certain frequency. For synaptic transmission to the postsynaptic neuron, the authors adopt the dynamic synapse formulation of Tsodyks and Markram Tsodyks, Pawelzik, and Markram (1998). It was shown that dynamic synapses featuring short-term synaptic plasticity may expand or shrink the interval of presynaptic firing rate over which ISR in the single postsynaptic neuron is present. Examining short-term depression and facilitation of the noisy postsynaptic current, other authors Uzuntarla et al. (2017) found that double inverse stochastic resonance (DISR) can occur with two distinct dips at different presynaptic firing rates. However, this research fails to incorporate a large range of other crucial plasticity principles that regulate the adaptability of neural networks in the brain. Clearly, there remains a large gap in understanding in literature how ISR non-adaptive neural networks are influenced by spike-timing-dependent plasticity (STDP) and homeostatic structural plasticity (HSP) — filling this gap is the aim of this study.

It is worth pointing out that in neurobiology, short-term synaptic plasticity (STP) Tsodyks, Pawelzik, and Markram (1998); Uzuntarla et al. (2017) and STDP Gerstner et al. (1996); Markram et al. (1997) are two distinct mechanisms that modulate synaptic strength, but they operate on different timescales, under different principles, and have different functional roles. In terms of timescale, STDP involves long-lasting changes, while STP involves transient changes. In terms of dependence on timing, STDP is highly dependent on the exact timing of pre-and postsynaptic spikes, while STP depends more on the recent history of synaptic activity and involves mechanisms such as neurotransmitter release probability. In terms of function, STDP is primarily associated with learning and memory, encoding long-term changes, whereas STP is involved in modulating synaptic transmission over short periods, affecting real-time signal processing.

Synaptic plasticity in neural networks denotes the ability to adjust the potency of synaptic links over time and/or transform the neural network’s structural configuration according to certain principles. Two primary mechanisms linked to adaptive regulations in neural networks are spike-time-dependent plasticity (STDP) and homeostatic plasticity (HSP). Synaptic modifications induced by STDP hinge on the repeated pairing of pre-and postsynaptic membrane potentials, with the extent and direction of these changes contingent on the precise timing of the neuronal firing. The exact timing of pre-and postsynaptic spikes determines whether synaptic weights undergo long-term depression (LTD) or long-term potentiation (LTP), which correspond to a lasting decrease or increase in synaptic strength, respectively Gerstner et al. (1996); Markram et al. (1997).

Synaptic modifications induced by HSP entail altering neuronal connectivity through the creation, pruning, or rearrangement of synaptic connections. This leads to modifications in the network’s architecture while preserving its operational framework, thereby enhancing the specialized functions of interconnected neuronal groups and enhancing the efficiency of sensory processing Shine et al. (2016). Initial indications of structural plasticity were identified through histological analyses of spine density in response to new sensory experiences or training Greenough and Bailey (1988). Additional studies revealed that the micro-connectome, which describes the connectome at the level of individual synapses, undergoes rewiring Bennett, Kirby, and Finnerty (2018); Van Ooyen and Butz-Ostendorf (2017); Yamakou and Kuehn (2023). Although brain networks conform to distinct topologies like small-world and random networks Hilgetag and Goulas (2016); Valencia et al. (2008) and exhibit dynamic behavior over time, recent research indicates that these networks can enhance information processing efficiency through homeostasis Butz, Steenbuck, and van Ooyen (2014). Motivated by these studies, this paper examines ISR in a time-varying small-world network of FitzHugh-Nagumo neurons evolving via STDP while adhering to its small-worldness via HSP at all times.

The main objectives of this study are the following. First, we examine how ISR is influenced by $\epsilon$, i.e., the timescale separation between the fast membrane potential and the slow recovery current variables of the neuron model; specifically, we investigate how the bistability necessary for the emergence of ISR is affected by $\epsilon$. Second, we study how STDP and HPS change the non-monotonic mean-firing rate response that is characteristic of ISR while varying the STDP control parameter $P$ (which determines whether STDP induces potentiation or depression-domination average synaptic weight) and the HSP control parameter $F$ (which determine how quickly the synapses of the small-world network architecture rewires while maintaining its small-world properties). To do this, we employed systematic and extensive numerical simulations to investigate these issues.

II Model

II.1 Neuron model

We consider a paradigmatic model with well-known biological relevance, the FitzHugh-Nagumo (FHN) neuron model Fitzhugh (1960); FitzHugh (1961); Nagumo, Arimoto, and Yoshizawa (1962):

$$\left\{\begin{array}[]{lcl}\displaystyle{\frac{dV_{i}}{dt}}&=&V_{i}\big{(}a-V_{i}\big{)}\big{(}V_{i}-1\big{)}-W_{i}+I_{i}^{\text{syn}}(t)+\eta_{i}(t),\\[8.53581pt] \displaystyle{\frac{dW_{i}}{dt}}&=&\epsilon\big{(}bV_{i}-cW_{i}\big{)}.\end{array}\right.$$ (3)

where $V_{i}=V_{i}(t)\in\mathbb{R}$ and $W_{i}=W_{i}(t)\in\mathbb{R}$ represent the fast membrane potential and slow recovery current variables of the neuron $i=1,\ldots,N$, respectively; $0<\epsilon\ll 1$ defines the timescale separation between $V_{i}$ and $W_{i}$; $a$, $b>0$ and $c>0$ are parameters changing the dynamic behavior of the neuron. The terms $\eta_{i}$ ($i=1,...,N$) are independent Gaussian noises with zero mean, standard deviation $\sigma$ (noise intensity), and correlation function $\langle\eta_{i}(t),\eta_{j}(t^{\prime})\rangle=\sigma^{2}\delta_{ij}(t-t^{\prime})$.

Note that, when we study a single neuron ($N=1$), we drop the subscripts $i$ and let the synaptic input $I^{\text{syn}}=0$.

II.1.1 Numerical Integration of SDE.

To integrate the stochastic differential equations (SDEs) Eqs. (3) in time, we used the Euler–Maruyama algorithm Higham (2001) with a small time step $\operatorname{d\!}t=0.0025$ and an integration time of $T=7.0\times 10^{3}$ units, which is sufficiently long to overcome transient behavior.

II.2 Network model

To include synaptic interactions in a neural network, we introduce the synaptic input current $I_{i}^{\text{syn}}(t)$ in Eq. (3), which models excitatory uni-directional chemical synapses between the neurons along their synaptic connections. The synaptic input current $I_{i}^{\text{syn}}(t)$ for the $i$th neuron at time $t$ is given by

$$I_{i}^{\text{syn}}(t)=-\frac{1}{k_{i}(t)}\sum_{i\neq j=1}^{N}\ell_{ij}g_{ij}s_{j}(t)\big{[}V_{i}(t)-V_{\text{syn}}\big{]}.$$ (4)

This term sums the synaptic input currents from all pre-synaptic neurons adjacent to neuron $i$. Such an interaction occurs if the neuron $j$ is pre-synaptic to the neuron $i$, i.e., if the connectivity matrix $L(=\{\ell_{ij}\})$ is $\ell_{ij}=1$ ; otherwise, $\ell_{ij}=0$. Specifically, we consider a small-world (SW) network Bassett and Bullmore (2006); Bassett et al. (2006); Liao, Vasilakos, and He (2017); Muldoon, Bridgeford, and Bassett (2016) constructed using a Watts-Strogatz network algorithm Watts and Strogatz (1998); Strogatz (2001), where the network’s Laplacian matrix $L(=\{\ell_{ij}\})$ is a zero-row-sum matrix. The sum over all synaptic inputs is normalized by the in-degree of the $i$th neuron (i.e., the number of synaptic inputs to the neuron $i$), $k_{i}(t)=\sum_{j\neq i}\ell_{ij}=k_{i}$. The matrix $g_{ij}$ represents the weight of the connection from the pre-synaptic neuron $j$ to the post-synaptic neuron $i$. Note that the connectivity matrix and the synaptic weights will adapt over time when we introduce plasticity mechanisms (see Sec. II.3).

An input current is modulated by the fraction of open synaptic ion channels, $s_{j}$, in a pre-synaptic neuron $j$. Finally, the membrane potential $V_{i}$ of the incident neuron $i$ is compared to the reverse potential $V_{\text{syn}}$.

The fraction of open synaptic ion channels at time $t$ of the $j$th neuron is represented by $s_{j}(t)$ in Eq. (4) and evolves in time according to Yu et al. (2015):

$$\frac{ds_{j}}{dt}=\frac{2(1-s_{j})}{1+\displaystyle{\exp\Bigg{[}-\frac{V_{j}(t)}{V_{\text{shp}}}\Bigg{]}}}-s_{j},$$ (5)

where $V_{j}(t)$ is the action potential of the pre-synaptic neuron $j$ at time $t$; $V_{\text{shp}}=0.05$ determines the threshold of the membrane potential above which the post-synaptic neuron $i$ is affected by the pre-synaptic neuron $j$.

Our study focuses on the inhibition of spiking activity triggered solely by noise through the effect of ISR. Nevertheless, it is well known that inhibitory synapses can also inhibit spiking activity in neural networks independently of ISR, operating through distinct mechanisms that can manifest both with and without the bi-stability required for ISR. Thus, to avoid the inhibition of spiking activity that inhibitory synapses might induce and focus only on the inhibition of spiking activity induced by the ISR effect, we shall fix the reversal potential $V_{\text{syn}}$ in Eq. (4) at $V_{\text{syn}}=2.0$ so that the network in Eq. (3) is entirely excitatory. Of course, one could study ISR in neural networks with inhibitory synapses, but for the reason given above, we are not interested in that case.

II.3 Adaptive network model

We intend to investigate the behavior of ISR in an adaptive network of $N$ FHN neurons exhibiting two forms of plasticity — spike-time-dependent plasticity (STDP) and homeostatic structural plasticity (HSP). Thus, the connectivity and synaptic weights are functions of time, i.e., $l_{ij}=l_{ij}(t)$ and $g_{ij}=g_{ij}(t)$, which are updated according to the rules of STDP and HSP explained in the following.

II.3.1 Spike-time-dependent plasticity (STDP)

The synaptic strength $g_{ij}(t)$ for each synapse is updated according to a nearest-spike pair-based STDP mechanism Morrison, Aertsen, and Diesmann (2007). The update rule according to Xie, Gong, and Wang (2018) is then implemented as follows:

$$\left\{\begin{array}[]{lcl}g_{ij}(t+\Delta t)=g_{ij}(t)+\Delta g_{ij},\\[8.53581pt] \Delta g_{ij}=g_{ij}(t)M(\Delta t),\\[8.53581pt] M(\Delta t)=\left\{\begin{array}[]{ll}\displaystyle{A\exp{(-\lvert\Delta t\rvert/\tau_{a})}\>\>\text{if}~{}\Delta t>0,}\\[2.84526pt] \displaystyle{-B\exp{(-\lvert\Delta t\rvert/\tau_{b})}\>\>\text{if}~{}\Delta t<0,}\\[2.84526pt] 0\>\>\text{if}~{}\Delta t=0.\end{array}\right.\end{array}\right.$$ (12)

This rule updates the synaptic coupling strength $g_{ij}(t)$ multiplicatively via the synaptic modification function $M$, where $\Delta t=t_{i}-t_{j}$, $t_{i}$ and $t_{j}$ represent the spiking time of neurons $i$ and $j$. The amount of synaptic modification $M$ is controlled by the adjusting potentiation and depression rate parameters represented by $A$ and $B$, respectively. The potentiation and depression temporal windows of the synaptic modification are controlled by $\tau_{a}$ and $\tau_{b}$, respectively.

Studies conducted experimentally have shown that the timeframe during which synaptic weakening occurs aligns closely with that of synaptic strengthening Bi and Poo (1998); Feldman and Brecht (2005); Song, Miller, and Abbott (2000). For our model, synaptic potentiation reliably occurs when the post-synaptic spike occurs within a 2.0 time unit window following the pre-synaptic spike, while depression is induced conversely. Thus, the temporal window parameters for potentiation and depression are fixed at the same value, i.e., $\tau_{a}=\tau_{b}$ = 2.0. The same studies have also shown that the ratio of the adjusting depression and potentiation rate parameters determines whether STDP exhibits long-term potentiation (LTP) or long-term depression (LTD). It was shown that STDP is depression-dominated if $P:=B/A>1$ and if $P:=B/A<1$, it is potentiation-dominated Bi and Poo (1998); Feldman and Brecht (2005); Song, Miller, and Abbott (2000).

In the present study, we wish to investigate ISR when STDP is both depression- and potentiation-dominated. We, therefore, keep the depression rate parameter fixed at $B=0.5$ so that the potentiation rate parameter is always given by $A=0.5/P$. Thus, we may consider $P$ as the single control parameter so that STDP is depression-dominated when $P>1$ Yamakou and Kuehn (2023); Yamakou, Desroches, and Rodrigues (2023), and potentiation-dominated when $P<1$. We vary $P$ in the interval given by $[5.0\times 10^{-6},5.0]$.

Furthermore, we wish to prevent (i) unbounded growth; (ii) negative coupling strengths, as it may give rise to inhibitory synapses, which we wish to avoid (for the reason given earlier); and (iii) the complete elimination of synapses (i.e., $g_{ij}=0$). To achieve this, we require that $g_{ij}$ remains bounded, i.e., $g_{ij}\in[g_{\text{min}},g_{\text{max}}]=[0.5\times 10^{-3},0.1\times 10^{-2}]$. Here, the maximum synaptic weight $g_{\text{max}}$ represents the value above which the bi-stability between the stable fixed point and limit cycle (indispensable for the onset of ISR) disappears and leaves only an unstable fixed point and a stable limit cycle. We choose a small but non-zero $g_{\text{min}}$ to prevent the complete deletion of synapses while allowing room for synaptic modifications that do not exceed $g_{\text{max}}$. Moreover, the initial weight of all excitable synapses is normally distributed in the interval $[g_{\text{min}},g_{\text{max}}]$, with mean $g_{0}=0.75\times 10^{-3}$ and standard deviation $\sigma_{0}=0.15\times 10^{-3}$.

II.3.2 Homeostatic structural plasticity (HSP)

To mimic HSP in the neural network dynamics given by Eqs. (3) and  (4), we generate a time-varying small-world network with rewiring probability $\beta\in(0,1)$ that adheres to its small-worldness at all times during the integration interval. To achieve this, we implement the following process during the rewiring of synapses Yamakou and Kuehn (2023): To build an initial small-world network, we used the Watts-Strogatz algorithm Watts and Strogatz (1998); Strogatz (2001) with rewiring probability of $\beta=0.25$ and average degree of $\langle k\rangle=4$. A synapse between two distant neurons is rewired to one of the neuron’s nearest neighbor with probability $(1-\beta)F\operatorname{d\!}t$. If the synapse is already between two nearest neighbors, it is replaced by a synapse to a randomly chosen distant neuron with probability $\beta F\operatorname{d\!}t$. We consider a node $i$ to be distant to node $j$ if $\lvert i-j\rvert>\langle k\rangle$, where $\langle k\rangle$ is the average degree of the original ring network used in the Watts-Strogatz algorithm Watts and Strogatz (1998); Strogatz (2001) to generate the initial small-world topology.

With a small integration time step of $\operatorname{d\!}t=0.0025$ and a rewiring probability of $\beta=0.25$, the parameter $F$ determines whether or not the probabilities given by $(1-\beta)F\operatorname{d\!}t$ and $\beta F\operatorname{d\!}t$ are large enough for the original and subsequent small-world networks to rewire as time advances in steps of $\operatorname{d\!}t$. Thus, the parameter $F$ becomes a proxy for the frequency at which the neurons in a small-world network swap their synapses while preserving the network’s small-worldness. In this paper, we call $F$ the characteristic frequency (which we will measure in per second (Hz) to make the probabilities dimensionless) of the time-varying network topology.

If $F=0$, then $(1-\beta)F\operatorname{d\!}t=0$ and $\beta F\operatorname{d\!}t=0$, and none of the synapses will be rewired. Consequently, the small-world network is time-invariant (static network). As soon as $F>0$, $(1-\beta)F\operatorname{d\!}t>0$ and $\beta F\operatorname{d\!}t>0$, and there is a non-zero probability that the network rewires and becomes a time-varying network. However, if $F$ is small, the topology will only slowly change over time. As $F$ increases, the network rewires faster because the probabilities $(1-\beta)F\operatorname{d\!}t$ and $\beta F\operatorname{d\!}t$ also increase. For example, if the probability that a synapse between two distant neurons is rewired to a nearest neighbor of one of the neurons is unity, i.e., $(1-\beta)F\operatorname{d\!}t=1$, then we compute the maximum rewiring frequency of the network as $F=1/((1-\beta)\operatorname{d\!}t)$, which is $\approx 533$ Hz  for the value of $\beta=0.25$ and $\operatorname{d\!}t=0.0025$ used in our computations. At the same time, the probability that a synapse replaces the synapse between two nearest neighbors to a randomly chosen distant neuron is computed as $\beta F\operatorname{d\!}t=0.25\times 533\times 0.0025\approx 0.33$.

Finally, note that we lose the time dependence in the average degree connectivity (average number of synaptic inputs per neuron) $\langle k\rangle$, because neurons would be able to change their neighbors via the rewiring rules, but we require that they do not change the number of neighbors. That is, the number of neighbors is always fixed, but the individual neighbors may be swapped over time with a certain probability.

III Numerical observations/measurements

Ensemble averages.

To quantify the dynamic behavior of Eq. (3) and to ensure robust statistical results of the neural network, for a fixed parameter, we calculate value ensemble averages over $R=300$ (independent) realizations with random initial conditions and initial small world network structure. For each realization of the neural network, initial conditions $(V_{i}(0),W_{i}(0))$ of the $i$th neuron ($i=1,...,N$) were drawn randomly from a uniform distribution within the range covering the basins of attraction of the stable fixed point and the limit cycle, i.e., $V_{i}(0)\in(-0.5,1)$, $W_{i}(0)\in(-0.05,0.2)$.

Measurements.

Measurements are always made by excluding a transient time of $T_{0}=1.0\times 10^{3}$ units. For each realization, we then calculated the number of spikes $n_{\text{spike},\ell}$ that occur during the remaining $T-T_{0}$ time units. Spikes were recorded when the membrane potential variable $V_{i}(t)$ crosses the threshold $V_{\mathrm{th}}=0.25$ from below.

An ensemble average firing rate $\langle r_{i}\rangle$ was calculated for each individidual neuron $i$ in the network (Eq. (3)) (or the single neuron in Eq. (3)) as follows:

$$r_{i}=\frac{1}{R(T-T_{0})}\sum\limits_{\ell=1}^{R}n_{\text{spike},\ell}.$$ (13)

The collective mean firing rate $\langle r\rangle$ of the neural network of $N=70$ neurons was then calculated as

$$\langle r\rangle=\frac{1}{N}\sum\limits_{i=1}^{N}r_{i}.$$ (14)

Spike-time-dependent plasticity (STDP).

To investigate how the average coupling strength $\langle g_{ij}\rangle$ of the network changes with the STDP parameter, which affects ISR, we averaged the synaptic weights over the entire population and time:

$$\langle g_{ij}\rangle=\displaystyle{\Bigg{\langle}\frac{1}{N^{2}}\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}g_{ij}(t)\Bigg{\rangle}_{t}},$$ (15)

where $\big{\langle}\cdot\big{\rangle}_{t}$ represents the average over the time interval $[T_{0},T]$.

In the following section on Results, we use the mean firing rate $\langle r\rangle$ of Eq. (14) to study the effect of (i) the bi-stability parameter $\epsilon$, (ii) the noise intensity $\sigma$, (iii) the STDP rule (controlled by the parameter $P$ defined earlier), and (iv) the HSP rule (controlled by the characteristic frequency parameter $F$) on the occurrence ISR. For an example of the control flow used in the simulations, see Appendix in Sec. VI. This flow of control is easily adapted to produce the rest of the simulations presented in this paper.

IV Results and discussion

IV.1 Bifurcation analysis for a single noiseless neuron

We provide a brief bifurcation analysis for Eq. (3) without noise ($\sigma=0$). An important goal of this analysis is to pinpoint the conditions on the parameters that enable the co-existence of a stable fixed point and a stable limit cycle, resulting in a bistable regime, the presence of which is essential to observe ISR.

For the FHN neuron in Eq. (3), there is a unique fixed point $(V_{0},W_{0})=(0,0)$ if and only if $(a-1)^{2}/4<b/c$. This fixed point is stable when $-a/\epsilon<c$ and $a>-b/c$, that is, when $a<0$ is sufficiently close to zero ($-1\ll a<0$), and in the limit $\epsilon\to 0$ only for $a\geq 0$.

The $V$-nullcline given by the graph of $W=-V^{3}+(a+1)V^{2}-aV$ loses normal hyperbolicity at its maximum $V_{+}$ and minimum $V_{-}$ points (i.e., the fold points), each located at $V_{\pm}=(a+1)/3\pm\sqrt{(a+1)^{2}/9-a/3}$. Here, it is worth noting that $V_{-}<V_{0}=0$ if and only if $a<0$.

When $c^{2}<b/\epsilon$ and $3\epsilon c\leq a^{2}-a+1$, we observe a Hopf bifurcation at $V_{\text{HB}}=(a+1)/3-\sqrt{(a+1)^{2}/9-(a+\epsilon c)/3}$. Since $\epsilon c>0$, it is easy to see that $V_{-}<V_{\text{HB}}$, and consequently, the Hopf bifurcation $V_{\text{HB}}$ is to the right of the minimum $V_{-}$ of cubic $V$-nullcline, hence on its ascending branch.

Whenever a fixed point $(V_{0},W_{0})=(0,0)$ is on the left descending branch of the $V-$nullcline, that is, to the left of the minimum $V_{-}$, it is stable, and this stability persists a little into the ascending branch for specific choices of the parameter values $a$, $b$, $c$, and $\epsilon$. For $a<0$, $V_{-}<V_{0}=0$, a Hopf bifurcation occurs when $\epsilon=-a/c=:\epsilon_{\text{HB}}$. Thus, the fixed point $(V_{0},W_{0})=(0,0)$ loses stability via Hopf bifurcation as $\epsilon$ decreases. The stable fixed point $(V_{0},W_{0})$ and an unforced stable limit cycle $\big{[}\bar{v}(t),\bar{w}(t)\big{]}$ co-exist as long as $-a/c<\epsilon$. In other words, we have bi-stability between the fixed point $(V_{0},W_{0})$ and a limit cycle if and only if $V_{-}<V_{0}=0<-a/c=\epsilon_{\text{HB}}$.

Fig. 1 illustrates the dynamic behavior of a single deterministic FHN neuron, where we chose parameters (see caption) such that bistable behavior is present, i.e., $V_{-}=-0.025244<V_{0}=0<-a/c=\epsilon_{\text{HB}}=0.025>0$. The bifurcation diagram in panel (a) displays a stable fixed point at $(V_{0},W_{0})=(0,0)$ and stable limit cycle for the parameter range $\epsilon\in[\epsilon_{\text{HB}},\epsilon_{\text{FB}})=[0.025,0.027865)$, where $\epsilon_{\text{FB}}$ represents the fold bifurcation point at which the stable and unstable limit cycles coalescence and annihilate, leaving behind only the stable fixed point $(V_{0},W_{0})=(0,0)$. The phase portrait in panel (b) shows the $V$- and $W$-nullclines in black and green, respectively. The stable fixed point $(V_{0},W_{0})=(0,0)$ (blue dot) lies to the right of the minimum of the $V-$nullcline and is surrounded by an unstable limit cycle (red) and a stable limit cycle (blue). Throughout the remainder of this paper, we fix the parameter values for the single and coupled FHN neurons such that each neuron resides within the bistable regime; i.e., we set $a=-0.05$, $b=1.0$, and $c=2.0$.

Figure 1: (a) Bifurcation diagram for $N=1$ neuron with $I_{\text{syn}}=0$, Eq. (3). Membrane potential $V$ vs. timescale parameter $\epsilon$ of a single noiseless FHN neuron. The bistable regime with a stable fixed point at $V=0$ (blue line) and stable limit cycle (blue dots) resides in the interval $\epsilon\in[0.025,0.027865)$. An unstable limit cycle (red dots) separates the two stable attractors. (b) Phase portrait of a single FHN neuron with $\varepsilon=0.02785$ displays the bistability between the fixed point at the origin (the unique intersection of the cubic $V$- and $W$-nullclines) and limit cycle (blue). An unstable limit cycle (red) separates the two attractive states. Other parameters: $a=-0.05$, $b=1.0$, $c=2.0$, $\sigma=0.0$.

IV.2 ISR in a single FHN neuron

We illustrate in Fig. 2 how different noise intensities impact the spiking behavior of a solitary FHN neuron of Eq. (3). Initial conditions lie in the basin of attraction of the stable limit cycle (i.e., $V_{i}=1.0$ and $W_{i}=0.2$), and we examine a range of timescale parameter values $\epsilon$ within the bi-stability interval $[0.025,0.027865)$.

Figure 2: (a) Time series of a single FHN neuron of the membrane potential $V$ with $\epsilon=0.0266$ for three different values of the noise amplitude $\sigma$ indicated. (b) ISR in a single FHN neuron for moderate noise level is characterized by the non-monotonicity of mean firing rate $\langle r\rangle$ curves. The dependence of mean firing rate $\langle r\rangle$ on noise amplitude $\sigma$ for different timescale separation parameters $\epsilon$. $a=-0.05$, $b=1.0$, $c=2.0$.

When $\epsilon=0.0266$, the time series of the membrane potential $V$ in Fig. 2(a) shows that there is an intermediate noise intensity that induces quiescence in the spiking activity, see the black time series with $\sigma=1.6\times 10^{-6}$. In Fig. 2(b), the behavior of the mean firing rate $\langle r\rangle$ is shown with respect to varying both the noise intensity $\sigma$ and the timescale parameter $\epsilon$. The non-monotonic behavior of the mean firing rate $\langle r\rangle$ as $\sigma$ and $\epsilon$ vary is characteristic of ISR.

A stronger ISR effect is associated with a deeper minimum in the non-monotonic $\langle r\rangle$ curve. Inspecting Fig. 2(b), it becomes evident that the closer the timescale parameter $\epsilon$ is to its Hopf bifurcation value $\epsilon_{\text{HB}}=-a/c=0.025$, the weaker is the ISR effect, see, e.g., the magenta curve with $\epsilon=0.02501$. As $\epsilon$ increases within the bistable interval, the non-monotonic $\langle r\rangle$ curve becomes deeper, indicating a stronger ISR effect; see, e.g., the purple curve in Fig. 2(b) with $\epsilon=0.0272325$.

We also notice that when $\epsilon$ is outside the bistable interval $[0.025,0.027865)$, e.g., for $\epsilon=0.0278650$, the effect of ISR disappears with the disappearance of the bi-stability between the fixed point $(V_{0},W_{0})=(0,0)$ and the limit cycle. For $\epsilon\geq 0.0278650$ (see Fig. 2(b)), only the stable fixed point $(V_{0},W_{0})=(0,0)$ remains, as the stable and unstable limit cycles have collided and annihilated each other. Consequently, increasing the noise intensity $\sigma$ results in a monotonic, rather than non-monotonic, increase in the mean firing rate $\langle r\rangle$.

The strengthening of the ISR effect with increasing $\epsilon\in[0.025,0.027865)$ can be explained in terms of the basin of attractions. Near the Hopf bifurcation threshold $\epsilon_{\text{HB}}=0.025$, the basin of attraction for the fixed point $(V_{0},W_{0})=(0,0)$ is significantly smaller in comparison to that of the stable limit cycle. Consequently, trajectories tend to swiftly depart from the basin of the fixed point while lingering longer within the basin of the limit cycle. This yields, as the noise intensity increases, a shallow $\langle r\rangle$ curve, as depicted by the magenta curve in Fig. 2(b) for $\epsilon=0.02501$.

As $\epsilon\in[0.025,0.027865)$ increases, the basin of attraction of the stable fixed point expands while that of the stable limit cycle contracts (eventually vanishing at $\epsilon=0.027865$). Consequently, the exit times from the basin of the fixed point lengthen compared to those from the limit cycle, leading to non-monotonic $\langle r\rangle$ curves characterized by deeper minima as $\epsilon\in[0.025,0.027865)$ increases.

IV.3 ISR in the adaptive network

IV.3.1 Effect of $\epsilon$ on ISR with STDP only

Since the timescale parameter $\epsilon$ controls the size of the basin of attraction of both the stable fixed point and the limit cycle in the isolated neuron, it is natural to investigate how $\epsilon$ affects the degree of ISR in an adaptive network. To do this, we computed the mean firing rate $\langle r\rangle$ in the STDP-driven SW network without HSP (i.e., $F=0$ Hz) while varying $\epsilon$, see Fig. 3. We found that the ISR phenomenon occurs in a wider range of $\epsilon$ when compared to a single neuron: (i) for $\epsilon<0.025$, where the single FHN neuron possesses a globally stable limit cycle (Fig. 1), the adaptive network exhibits ISR that is more pronounced than that of the single neuron (Fig. 2(b)); (ii) for $\epsilon>0.027865$, where the single FHN neuron has a globally stable equilibrium point, typical ISR curves can be observed for the adaptive network even for $\epsilon$ away from the fold bifurcation point of the limit cycles.

The dynamic behavior of neurons in vivo has to be characterized as a collective phenomenon rather than in isolation, and our simulation describes indeed the setting of a neural network. Thus, our results (see Fig. 3) strongly suggest that ISR should be more pronounced and more readily observed in experiments involving networks of neurons when compared to isolated neurons Buchin et al. (2016).

Figure 3: Collective firing behavior of a SW network with STDP (and without HSP) for different timescale parameter values $\epsilon$ and noise amplitudes $\sigma$. The emergence of ISR in this network is characterized by the non-monotonicity of the mean firing rate $\langle r\rangle$ curves that occur at intermediate noise levels. Stronger ISR occurs at the largest value of $\epsilon$ in the bistable regime $\epsilon\in[0.024,0.0290]$. $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $P=5.0\times 10^{-6}$, $A=B/P$, $F=0.0$, $N=70$.

Our observations for these simulations are qualitatively similar to Fig. 2(b), i.e., the minimum mean firing rate $\langle r\rangle_{\text{min}}$ decreases as $\epsilon$ increases, thus indicating an enhanced inhibition of the spiking activity by ISR. However, this effect gets stronger when $\epsilon\geq 0.027$ and $\langle r\rangle_{\text{min}}$ becomes very low around $\epsilon\approx 0.028$ — close to the upper boundary of the bi-stability region (see Fig. 4 (a)). This behavior can be explained by the rapid coalescence of the stable and the unstable limit cycles at the fold bifurcation as $\epsilon$ increases, see Fig. 1. The ensuing rapid decrease in the size of the basin of attraction of the limit cycle also induces a rapid decrease in the optimal noise intensity $\bar{\sigma}$ required to achieve $\langle r\rangle_{\text{min}}$ (see, e.g., Fig. 4(b)). The plot of $\langle r\rangle$ in the $(\epsilon,\sigma)$-plane (Fig. 4(c)) clearly reveals this effect at the boundary, in agreement with Fig. 4(a) and (b).

Figure 4: Absence and appearance of ISR in an STDP-driven small-world network without HSP. (a) The minima of the mean firing rate $\langle r\rangle_{\text{min}}$ vs. $\epsilon$. (b) Corresponding noise intensity $\bar{\sigma}$ to $\langle r\rangle_{\text{min}}$ response to the timescale parameter $\epsilon$. (c) Mean firing rate $\langle r\rangle$ in dependence of $\epsilon$ and $\sigma$ plane displays ISR for $\epsilon\in[0.024,0.029]$. Parameters are: $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $P=5.0\times 10^{-6}$, $F=0.0$, $A=B/P$, $N=70$.

In Fig. 5 depict examples of dynamic behavior of the network for $\epsilon=0.0285$ with three different noise intensities. For small and large noises, the neurons in the adaptive networks spike frequently. We observe a quasi-total inhibition of the spiking activity due to ISR for the intermediate noise.

Figure 5: Dynamic behavior of neurons in the STDP-driven SW network without HSP at three different noise amplitudes for $\epsilon=0.0285$: (a) $\sigma=1.0\times 10^{-8}$; (b) $\sigma=2.5\times 10^{-6}$; (c) $\sigma=1.25\times 10^{-4}$. $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $P=5.0\times 10^{-6}$, $A=B/P$, $F=0.0$, $N=70$.

IV.3.2 Effect of HSP on ISR

We study the effects of only HSP on ISR, see Fig. 6. Recall that the HSP algorithm rewires the synaptic connections between neurons while preserving the small-worldness of the network topology (see Sec. III). The rewiring frequency $F$ determines how fast the SW network reshapes its synaptic connections. We investigate a wide range of from $F$ $0$ Hz to $500$ Hz, which covers the limiting regimes where the network is static (i.e., rewires with probability zero) up to where the rewiring of synapses between two distant neurons occurs with probability close to unity (i.e., rewires at each time step and hence, changes the topology very quickly). Fig. 6 illustrates how $F$ affects the mean field frequency $\langle r\rangle$ and ISR, as we vary the rewiring frequency $F$ while keeping all other parameters fixed. These simulations are carried out over a range of $\epsilon$-values within the bistable regime of the network.

The overall trend is that faster rewiring implies a lower ISR curve; however, this behavior may be more or less pronounced as we vary $\epsilon$. Specifically, when $\epsilon$ is close to the lower boundary of the bi-stability region (see Fig. 6(a) where $\epsilon=0.02501$), the enhancement of ISR by increasing $F$ is inconspicuous compared with when $\epsilon$ is chosen higher up in the bistable interval. With larger values of $\epsilon$ as in Figs. 6(b)-(d), the enhancement of ISR becomes stronger for higher rewiring frequency $F$. Furthermore, we observe (especially in Fig. 6(b)) that the optimal noise intensity, i.e., where the mean firing rate is most inhibited and ISR is most pronounced, is slightly shifted to the left towards smaller $\sigma$ values. Thus, we observe that HSP results in two effects: (i) larger $F$ lowers the ISR curves, and (ii) larger $F$ slightly shifts the minima of the ISR curves to the left, making them occur at slightly smaller values of the noise intensity $\sigma$.

Conceptually, we may explain these effects as follows. First, note that neurons in the network have independent noise sources and different initial conditions. Therefore, they can spike independently of each other and at different times. Second, fast rewiring (i.e., larger $F$) makes it even more difficult for the neurons to synchronize their spiking activity, as neurons constantly swap neighbors via HSP. When $F$ is large, connecting neurons do not have the time to synchronize their different spiking times before they become disconnected again via HSP. Hence, neurons in the network spiking at different rates would quickly and repeatedly connect and disconnect from each other as time evolves. The overall consequence of the quick connections and disconnections between neurons with varying spiking rates is the creation of a second source of synaptic noise to each neuron in the network. This second source of synaptic noise (induced and measured by $F$) combines with the synaptic noise (measured by $\sigma$) to increase the overall external stochastic forcing of the neurons involved in fast synaptic rewiring. Ultimately, this leads to a downward shift of the ISR curves, with minima occurring at lower values of the synaptic noise intensity $\sigma$ which explains effect (i).

Next, we note that larger $\epsilon$ expands (shrinks) the basin of attraction of the fixed point (limit cycle). Thus, a smaller noise intensity is required to kick trajectories out of the smaller basin of attraction of the limit cycle. Furthermore, the additional synaptic noise source (induced by the fast rewiring of the synaptic connections) makes it even easier for trajectories to escape from the basin of attraction of the limit cycle to that of the fixed point, as it assists the small synaptic noise intensity $\sigma$ in the escape process. This explains effect (ii).

Figure 6: Effect of HSP (controlled by $F$) on collective firing behavior of a SW network at four different values of $\epsilon$ within the bistable regime. Higher $F$ lowers the $\langle r\rangle$ curves, thereby enhancing ISR. However, this effect is significant for larger values of $\epsilon$ and weak noise amplitudes. $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $P=5.0\times 10^{-6}$, $A=B/P$, $N=70$.

IV.3.3 Effect of STDP on ISR

We now set the rewiring frequency to $F=0$ and switch off HSP to study the effect of STDP on ISR only. Recall that STDP can be controlled by the parameter $P$, which represents the ratio between the adjusting depression $B$ and potentiation $A$ rate parameters. By fixing $B=0.5$ so that $A$ is given by $A=0.5/P$, we vary $P$ so that STDP is either depression-dominated (i.e., $P>1$) or potentiated-dominated (i.e., $P<1$). In this section of the paper, we are interested in the effects of STDP on ISR when $P$ varies in $[5.0\times 10^{-6},5.0]$. Here, we also vary $\epsilon$ and consider these effects in various parts of the bistable regime. The results are shown in Figs. 7-9, respectively.

In Figs. 7(a) and (b) with $\epsilon=0.02501$, ISR is not significantly affected by changing values of $P$. This less pronounced effect of ISR as $P$ changes is an immediate consequence of the small basin of attraction of the fixed point when the timescale parameter $\epsilon$ is near $\epsilon_{\text{HB}}=0.02501$. Nevertheless, we can still see from the inset of Fig. 7(a) that the larger value of $P$ (see, e.g., the green curve with $P=5.0=P_{\text{opt}}$) induces a slightly deeper minimum of the $\langle r\rangle$ curve when compared to smaller values of $P$, i.e., ISR is slightly enhanced.

Figure 7: Effect of STDP (controlled by $P$) on collective firing of the SW network at $\epsilon=0.02501$ (panels (a) and (b)). (c) Associated population- and time-averaged synaptic weights of the network $\langle g_{ij}\rangle$ vs. noise $\sigma$ for two non-optimal values and the optimal value of $P$. Parameters are: $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $F=0.0$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $A=B/P$, $N=70$.

In Figs. 8(a) and (b), setting an intermediate value of $\epsilon=0.0275375$ within the bistable interval, we observe that changing $P$ has a more significant effect on ISR when compared to $\epsilon=0.02501$. Larger values of $P$ enhance ISR up to a certain threshold, upon which increasing $P$ further enhances ISR no more. More precisely, there is a significant downward shift in the ISR curve as $P$ changes from $5.0\times 10^{-5}$ to $5.0\times 10^{-4}$, after which ISR cannot be enhanced further.

Figure 8: Effect of STDP (controlled by $P$) on collective firing behavior of the SW network at $\epsilon=0.0275375$ in (a) and (b). Associated population- and time-averaged synaptic weights of the network $\langle g_{ij}\rangle$ vs. noise $\sigma$ for a non-optimal value and the optimal value of $P$ in (c). $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $F=0.0$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $A=B/P$, $N=70$.

In Figs. 9(a) and (b), we set $\epsilon=0.0290$ close to the upper bound of the networks’ bistable interval. Contrasting Figs. 7(a) and (b) and Figs. 8(a) and (b), we now observe two different behaviors: (i) increasing $P$ has a more significant effect on ISR, especially at weaker noise intensities ($\sigma\in[10^{-7},10^{-6}]$); and (ii) the lowest ISR curve is achieved for an intermediate value of $P=5.0\times 10^{-4}$ (rather than the largest value of $P=5.0$).

Figure 9: Effect of STDP (controlled by $P$) on collective firing behavior of the SW network at $\epsilon=0.0290$ in (a) and (b). Associated population- and time-averaged synaptic weights of the network $\langle g_{ij}\rangle$ vs. noise $\sigma$ two non-optimal and the optimal values of $P$ in (c). $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $F=0.0$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $A=B/P$, $N=70$.

To gain further insight into the behaviors shown in Figs. 7(a) and (b)-9(a) and (b), we computed the corresponding variations of the average synaptic weight $\langle g_{ij}\rangle$ of the network as a function of the noise intensity $\sigma$. The results are shown in Figs. 7(c)-9 (c). One sees that when $P$ is at its optimal value [i.e., $P_{\text{opt}}=5.0$ in Figs. 7(a) and 8(a), and $P_{\text{opt}}=5.0\times 10^{-4}$ in Fig. 9(a)], the values of the average synaptic weight $\langle g_{ij}\rangle$ in the intervals of the noise intensity in which lowest ISR curves are achieved [i.e., $\sigma\in(1.25\times 10^{-5},4.0\times 10^{-5})$, $\sigma\in(3.0\times 10^{-6},1.25\times 10^{-5})$, and $\sigma\in(5.0\times 10^{-7},6.0\times 10^{-6})$ in Figs. 7(a)- 9(a), respectively] are lower than those computed at the non-optimal values of $P$. The behavior can be explained by the fact that larger (smaller) values of average synaptic weight $\langle g_{ij}\rangle$ induced by the largest (smallest) value of $P$ in Figs. 7(a)-8(a) and an intermediate value of $P$ in Fig.9(a)] establishes a stronger (weaker) coupling between spiking and quiescent neurons, thereby enhancing (inhibiting) the recruitment of quiescent neurons into the spike state leading to the deterioration (improvement) of ISR.

IV.3.4 Combined effects of STDP and HSP on ISR

Our previous investigations indicated that for intermediate values near $\epsilon=0.0275375$ within the network’s bistable regime and intermediate values of the noise intensity [i.e., $\sigma\in(3.0\times 10^{-6},1.25\times 10^{-5})$], the effect of the HSP parameter $F$ and the STDP parameter $P$ on ISR becomes significant compared to the lower and higher values of $\epsilon$. In particular, it is seen in Figs. 6(b) and 8, where $\epsilon=0.02753755$, that increasing $F$ nd $P$, respectively, lowers the $\langle r\rangle$ curves in the noise interval $\sigma\in(3.0\times 10^{-6},1.25\times 10^{-5})$, indicating an enhancement of ISR in each case.

A natural question arises: Can increasing the HSP and STDP control parameters jointly enhance ISR beyond the level achieved when only one of these parameters is increased? And, if so, which of these parameters has the greater impact on enhancing ISR? To answer this question, in Fig. 10(a), we examine the joint effect of STDP and HSP on ISR for $\epsilon=0.0275375$. It can be seen that the deepest $\langle r\rangle$ curve is achieved when $F$ and $P$ are at their largest values — compare the black $\langle r\rangle$ curve, with $P=5.0$ and $F=500$Hz, to the rest of the other $\langle r\rangle$ curves. Clearly, increasing $F$ and $P$ can enhance ISR beyond the level of enhancement induced when just one of these parameters is increased to a larger value.

Furthermore, in Fig. 10(a), the separation between the minimum $\langle r\rangle_{\text{min}}$ of the red curve (with $P=5.0\times 10^{-6}$ and $F=0$Hz) and the minimum $\langle r\rangle_{\text{min}}$ of the pink curve (with $P=5.0$ and $F=0$Hz) is 3.61. While the separation between the minimum $\langle r\rangle_{\text{min}}$ of the red curve and the minimum $\langle r\rangle_{\text{min}}$ of the blue curve (with $P=5.0\times 10^{-6}$ and $F=500$Hz) is 2.09. Hence, increasing the STDP control parameter $P$ has a stronger effect on ISR than increasing the HSP control parameter $F$.

To obtain deeper insight into the behaviors depicted in Figs.10(a), we calculated the corresponding changes in the average synaptic weight $\langle g_{ij}\rangle$ of the network as a function of noise intensity $\sigma$. The results are shown in Fig. 10(b). We observe that the noise interval, $\sigma\in(3.0\times 10^{-6},1.25\times 10^{-5})$, where the mean firing rate $\langle r\rangle$ curves in Fig. 10(a) reach their minimum, corresponds to the same noise interval in which the average synaptic weight $\langle g_{ij}\rangle$ curves attain their lowest value. Specifically, within the noise range $\sigma\in(3.0\times 10^{-6},1.25\times 10^{-5})$, a decrease in the mean firing rate $\langle r\rangle$ curve corresponds to a dip in the average synaptic weight $\langle g_{ij}\rangle$ curve — the deeper the $\langle r\rangle$ curve, the stronger the dip of the $\langle g_{ij}\rangle$ curve. This behavior can again be attributed to the fact that a weaker average synaptic weight $\langle g_{ij}\rangle$ inhibits the recruitment of quiescent neurons into the spiking state. In addition, because the basin of attraction of the stable fixed point is significantly larger (i.e., $\epsilon_{\text{HB}}<\epsilon=0.0275375$), once neurons in the spiking state are pushed into the resting state, they tend to remain at rest (at least for a very long time), thereby enhancing ISR.

At this point, it is worth mentioning that all the numerical simulations presented in this paper with a small-world network, incorporating STDP and/or HSP have also been implemented with an STDP-driven random neural network that adheres to its randomness at all times via HSP. All the numerical simulations (not shown) for the case of the random network with a different HSP rule yielded qualitatively similar results to those observed in the small-world neural network. To generate a time-varying random network topology (also generated with the Watts-Strogatz algorithm Strogatz (2001) for $\beta=1$) while keeping the statistical network structure constant (i.e., $\beta$), we implement the following process during the rewiring of synapses: During each integration time step $dt$, if there is a synapse between neuron $i$ and $j$, it will be rewired such that neuron $i$ ($j$) connects to any other neuron except for neuron $j$ ($i$) with a probability of $\big{(}1-\frac{\langle k\rangle}{N-1}\big{)}Fdt$.

Figure 10: Combined effect of STDP and HSP on collective firing behavior of the SW network in (a). Associated population- and time-averaged synaptic weights of the network $\langle g_{ij}\rangle$ vs. noise $\sigma$ for various combinations of $P$ and $F$ in (b). $a=-0.05$, $b=1.0$, $c=2.0$, $V_{\text{syn}}=2.0$, $V_{\text{shp}}=0.05$, $\beta=0.25$, $\langle k\rangle=4$, $\tau_{a}=\tau_{b}=2.0$, $B=0.5$, $A=B/P$, $N=70$.

V Summary and conclusions

In summary, we have conducted a numerical investigation into the phenomenon of inverse stochastic resonance (ISR) in a single FitzHugh-Nagumo (FHN) neuron and an adaptive small-world network of FHN neurons, under the influence of spike-time-dependent plasticity (STDP) and homeostatic structural plasticity (HSP). Through a combination of bifurcation analysis and numerical simulations, we identified the specific parameter values and intervals for which both the individual neuron and the small-world network exhibit a bistability regime, characterized by the co-existence of a stable fixed point and a limit cycle, a pre-requisite for the emergence of ISR.

The degree of ISR was shown to be highly dependent on the value of the timescale separation parameter $\epsilon$ of the model within the bistable interval. Using the mean firing rate to gauge the degree of ISR, we found that as $\epsilon$ approaches the Hopf bifurcation threshold—the lower bound of the bistability interval—ISR becomes less pronounced. Conversely, ISR becomes more pronounced as $\epsilon$ moves further away from this threshold.

Our computations demonstrated that, within the bistable interval, the effects of STDP and HSP on the degree of ISR vary depending on both the value of $\epsilon$ and the interval of synaptic noise intensity $\sigma$. The effects of STDP and HSP are less significant when $\epsilon$ is close to the lower bound of the bistability interval and become more significant for larger values of $\epsilon$ within this interval and intermediate values of the synaptic noise intensity. As $\epsilon$ gets closer to the upper bound of this bistability interval, ISR becomes stronger, especially at weaker synaptic noise intensities. The reason why ISR is enhanced as $\epsilon$ is increased (within the bistability interval) is attributed to the fact that the basin of attraction of the stable fixed point (limit cycle) grows (shrinks) as $\epsilon$ increases, thereby inhibiting the spiking activity of the neurons, at least for a very long time, especially at weak noise intensities.

More specifically, our results indicated that at intermediate values of the timescale parameter $\epsilon$, increasing the rewiring frequency $F$ for HSP may noticeably enhance the degree of ISR at intermediate synaptic noise intensities. When the timescale parameter $\epsilon$ is even closer to the upper bound of the bistability interval, augmenting $F$ enhances ISR, especially at weak noise intensities. Our rationale behind this behavior is that the fast rewiring of synapses constitutes a secondary source of noise (controlled by $F$), which — in addition to the synaptic noise — helps inhibit neural spiking and enhance ISR. Furthermore, our results indicated that at intermediate values of timescale parameter $\epsilon$, increasing the STDP control parameter $P$ (which determines whether the neural network exhibits a potentiation- or depression-dominated in terms of the population- and time-averaged synaptic weights of the network), noticeably enhances the effect of ISR within intermediate noise intensities. When the timescale parameter $\epsilon$ is closer to the upper bound of the bistability interval, the intermediate value of $P$ noticeably enhances ISR the most for weak noise intensities. Furthermore, our simulations indicated that the STDP control parameter $P$ has a greater ISR enhancement capability than the HSP control parameter $F$.

It is important to note that the results presented in this study may be influenced by the specific rewiring strategies employed to preserve the small-world characteristics of the time-varying small-world networks. Nevertheless, we expect that other rewiring rules that capture the nature of HSP (i.e., swapping synaptic connections while maintaining the small-world network topology) would lead to similar results. Our findings strongly suggest that the inherent synaptic noise, the timescale difference between the membrane potential and recovery current variables of neurons, and the prevalent STDP and HSP can collectively contribute to improving ISR.

An experimental study has demonstrated that ISR can lead to a locally optimal information transfer between the input and output spike trains of the Purkinje neurons Buchin et al. (2016). Recent experiments have shown that signaling molecules, such as acetylcholine, can modulate STDP Brzosko, Mierau, and Paulsen (2019). Advances in experimental neuroscience have facilitated the manipulation of synaptic control in the brain through drugs that affect neurotransmitters Pardridge (2012) or by using optical fibers to stimulate genetically engineered neurons selectively Packer, Roska, and Häusser (2013). Our findings yield practical implications in locally enhancing optimal information transfer between input and output spike trains via ISR in both experimental settings and artificial neural circuits. Thus, our study prompts exciting venues for further experimental work on ISR in neural networks. Finally, our findings suggest the possibility of developing ISR enhancement strategies for noisy artificial neural circuits, aimed at optimizing local information transfer between input and output spike trains in neuromorphic systems.

Author Declarations

Conflict of Interest

The authors have no conflicts to disclose.

Data Availability Statement

The simulation data supporting this study’s findings are available from the corresponding author upon reasonable request.

VI appendix

/* $X_{i}(t)$ $=\{v_{i}(t),\>w_{i}(t)\}$: variables of coupled SDEs in Eq.(3) ;

$t$: time;

$T$: total integration time;

$T_{0}$: Transient time;

$N$: network size;

$B$: adjusting depression rate parameter;

$\beta$: rewiring probability in Watts-Strogatz algorithm;

$R$: number of realizations;

$F$: HSP control parameter (rewiring frequency of synapses);

$P$: STDP control parameter;

$P_{\text{min}}$: min of $P$;

$P_{\text{max}}$: max of $P$;

$\epsilon$: timescale separation parameter;

$\sigma$: synaptic noise intensity;

$\sigma_{\text{min}}$: min of $\sigma$;

$\sigma_{\text{max}}$: max of $\sigma$;

$\ell_{ij}(t)$: adjacency matrix;

$n_{\text{spike}},i,m$ : number spikes of the $i^{th}$ neuron at the $m^{th}$ realization;

$\langle r\rangle_{i,m}$: meaning firing rate of the $i^{th}$ neuron at the $m^{th}$ realization;

$\langle r\rangle$: average of the meaning firing rate over $R$;

$g_{ij}(t)$: synaptic weights;

$\langle g_{ij}\rangle_{\text{m}}$: average of synaptic weights over $t$ and $N$ of the $m^{th}$ realization;

$\langle g_{ij}\rangle$: average of mean synaptic weights over $R$;

*/

Input: $T$, $T_{0}$, $N$, $R$, $F$, $P$, $\beta$

Output: $\langle r\rangle,\langle g_{ij}\rangle$

$P\leftarrow P_{\text{min}}$ ;

  // Initialize the adjusting rate parameter

1 while $P\leq P_{\text{max}}$  do

       $\sigma\leftarrow\sigma_{\text{min}}$ ;

        // Initialize the timescale parameter

2       while $\sigma\leq\sigma_{\text{max}}$  do

3             for $m\in 1,2,\dots,R$ do

                   Init $X_{i}(t)\>,\>\ell_{ij}(t)$ ;

                    // Random initial conditions of SDEs and initial SW network adjacency matrix

4                   for $t\in 0,\dots,T$ do

                         Integrate network of SDEs in Eqs. (3) ;

                          // Using the Euler-Maruyama method

                         Record the current voltage spike times $t^{n}_{i}$ from $V_{i}(t)$;

                          // Times $t$ at which $v_{i}(t)\geq v_{\mathrm{th}}=0.25$

5                         if $\Delta t_{ij}:=t_{i}-t_{j}>0$ then

                               $\Delta M\leftarrow\displaystyle{\frac{B}{P}}\exp{(-\lvert\Delta t_{ij}\rvert/\tau_{p})}$ ;

                                // $t_{i}$ , $t_{j}$: nearest-spike times of post ($i$) & pre ($j$) neuron

6                              

7                        

8                        if $\Delta t_{ij}<0$ then

9                               $\Delta M\leftarrow-B\exp{(-\lvert\Delta t_{ij}\rvert/\tau_{d})}$

10                        

11                        if $\Delta t_{ij}=0$ then

12                               $\Delta M\leftarrow 0$

13                        

                        $g_{ij}(t)\leftarrow g_{ij}(t)+g_{ij}(t)\Delta M$;

                          // update synaptic weights

                         $\ell_{ij}(t)\leftarrow\widetilde{\ell_{ij}}(t)$ ;

                          // Update the adjacency matrix with $\widetilde{\ell_{ij}}(t)$ obtained by randomly rewiring $\ell_{ij}(t)$ with frequency $F$ according rewiring algorithm that preserve the small-worldness of the SW network generated with rewiring probability $\beta$.

14                        

15                        if $t\geq T_{0}$ then

16                               $r_{i,m}\leftarrow\displaystyle{n_{\text{spike},i,m}}$  , $\langle g_{ij}\rangle_{\text{m}}\leftarrow\displaystyle{\bigg{\langle}N^{-2}\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}g_{ij}(t)\bigg{\rangle}_{t}}$;

17                              

18                        

19                  Add $r_{i,m}$ to $r_{i}$;

20                   Add $\langle g_{ij}\rangle_{\text{m}}$ to $\overline{\langle g_{ij}\rangle}$;

21                  

            $\langle r_{i}\rangle\leftarrow r_{i}[R(T-T_{0})]^{-1}$ ;

              // Compute averages over the $R$ realizations

22             Add $\langle r_{i}\rangle$ to $\overline{\langle r\rangle}$;

             $\langle r\rangle\leftarrow\overline{\langle r\rangle}/N$ ;

              // Compute averages over the $R$ realizations

             $\langle g_{ij}\rangle\leftarrow\overline{\langle g_{ij}\rangle}/R$ ;

              // Compute averages over the $R$ realizations

             $\sigma\leftarrow\sigma+\Delta\sigma$ ;

              // Increment the noise intensity

23            

      $P\leftarrow P+\Delta P$ ;

        // Increment the STDP control parameter

24      

Algorithm 1 Flow of control in the simulations of Fig.7. Adapted for the simulations in the rest of the figures.