^{1}

^{*}

^{1}

^{2}

^{*}

^{1}

^{2}

Edited by: Arpan Banerjee, National Brain Research Centre (NBRC), India

Reviewed by: Spase Petkoski, INSERM U1106 Institut de Neurosciences des Systèmes, France; John David Griffiths, University of Toronto, Canada

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Frequency-dependent plasticity refers to changes in synaptic strength in response to different stimulation frequencies. Resonance is a factor known to be of importance in such frequency dependence, however, the role of neural noise in the process remains elusive. Considering the brain is an inherently noisy system, understanding its effects may prove beneficial in shaping therapeutic interventions based on non-invasive brain stimulation protocols. The Wilson-Cowan (WC) model is a well-established model to describe the average dynamics of neural populations and has been shown to exhibit bistability in the presence of noise. However, the important question of how the different stable regimes in the WC model can affect synaptic plasticity when cortical populations interact has not yet been addressed. Therefore, we investigated plasticity dynamics in a WC-based model of interacting neural populations coupled with activity-dependent synapses in which a periodic stimulation was applied in the presence of noise of controlled intensity. The results indicate that for a narrow range of the noise variance, synaptic strength can be optimized. In particular, there is a regime of noise intensity for which synaptic strength presents a triple-stable state. Regulating noise intensity affects the probability that the system chooses one of the stable states, thereby controlling plasticity. These results suggest that noise is a highly influential factor in determining the outcome of plasticity induced by stimulation.

One hypothesis for the basis for communication in the brain is that it occurs

Plasticity is the mechanism by which the nervous system adapts to external stimuli over multiple spatial and temporal scales. The phenomenon can manifest in functional or structural changes to neural networks as we experience the world and underlies cognitive processes such as learning and memory. In a previous study, we applied a periodic stimulus at a range of frequencies to a neural network model based on Wilson-Cowan (WC) oscillators connected

Experimental evidence has shown the existence of both multistability and metastability in biological neural networks (see Tognoli and Kelso,

On the other hand, multistability is ubiquitous in biological systems in general (for example, Crabb et al.,

Moreover, noise in dynamical systems can have a significant role in determining the stable-state landscape and in allowing transitions (switching) between those states. Therefore, our aim was to characterize synaptic plasticity dynamics in a model of interacting neural populations as a function of stimulation frequency and noise intensity. To that purpose, we implemented a WC model of neural activity that included activity-dependent and homeostatic plasticity mechanisms.

The mathematical details describing the time-evolution of the activity for the E and I populations are given in Equations (1–3). In the equations, the subindex

Where _{EE}, _{EI}, _{IE}, _{II}. The model also receives zero-mean and unit-variance white-noise input, ξ_{i}(_{EE} = 23, _{II} = 0, _{IE} = 35, _{EI} = 15, as these have been shown to generate stable oscillatory behavior when coupled with the background activity levels given below, see Wang et al. (_{0}, _{0} for the E and I populations respectively and fixed as _{0} = 0.5, _{0} = −5, as in Lea-Carnall et al. (_{E/I}, and the weights _{ij} controlling excitatory connections are described in the next section. A schematic of the models are given in

Schematic diagrams for the networks used in analysis. _{EE}, _{IE}, _{II}, and _{EI}. Inter-unit connection strength is governed by connectivity parameters contained in the connectivity matrices _{jk} and _{jk}. Mathematical details for the model dynamics are given in Equations (1–6).

The WC unit has an intrinsic resonance frequency, _{r}, that is controlled _{E} and τ_{I} (Wang et al., _{r} = 4 Hz: τ_{E} = 0.017 s, τ_{I} = 0.013 s, _{r} = 8 Hz: τ_{E} = 0.024 s, τ_{I} = 0.014 s, _{r} = 12 Hz: τ_{E} = 0.011 s, τ_{I} = 0.007 s, and _{r} = 23 Hz: τ_{E} = 0.014 s, τ_{I} = 0.006 s, these were derived computationally.

Inter-unit connectivity structure and strength are contained in the quantities _{ij} and _{ij}. In particular, all connections to the inhibitory populations are set _{ij} = 0.1. For the connections to the excitatory population, an initial value is set for _{ij} = 0.15 at

Equation (4) specifies the learning rule which governs the dynamics of the excitatory connections _{ij}. In order for _{ij} to increase, the product of the rates _{i} and _{j} must be greater than a threshold _{ij} and the synapse is enhanced. Conversely, if _{i} = _{j} = 0 or if the product of the fractional firing rates _{i} and _{j} is less than the threshold _{ij} will decay toward zero. The nonlinear threshold is implemented by the Heaviside function Θ(_{E/I} is calculated for each population for each unit and is subtracted from the population's activity at each point in time, see Remme and Wadman (

Where τ_{h} = 2.5 s, γ = 1, _{∞} = 0.2, _{∞} = 0.2, τ_{SE} = 1 s, τ_{SI} = 2 s, and _{h} and γ were chosen so that the slope of the decay (when the two units were not coincidentally firing) was equal to that of the increase, as in Lea-Carnall et al. (

Each unit within the network received an independent Gaussian white-noise input, ξ (

In the case of the rhythmic input

The Euler-Murayama method was used for the integration of stochastic differential equations with an integration step of 1 ms (Higham, ^{5} time steps in the case of the 2-unit model and 10^{6} time-steps for the 10-unit model. The results given in ^{6} time-steps and the first 3 × 10^{5} were discarded, while in the second, the model ran for 2 × 10^{6} time-steps and 1 × 10^{5} were discarded. Transmission delays were assumed to be instantaneous, although we note that in models of STDP, delays within the network have been shown to influence plasticity outcomes (Madadi Asl et al.,

Initially, each time series was band-pass filtered ± 5 Hz around _{r}, Hilbert transformed, and then the instantaneous phase was extracted for each point in time

The complex phase locking value (

Where Δϕ(_{1}(_{2}(

In

The Kuramoto order parameter (Kuramoto,

Where

All simulations were performed in Matlab (The Mathworks Inc., MATLAB ver. R2019b).

The WC model, a neural mass model, is used to represent each unit (or mass) within our network (Wilson and Cowan, _{ij}, for excitatory connections and _{ij}, for inhibitory connections which are fixed connecting unit

We first consider a model consisting of two WC units, with _{12} and _{12} being the only non-zero synapses. Since there is only one excitatory weight _{12}, from this point we omit the subscript and refer to it as _{d}) for different resonant frequencies of the units (_{r}) and levels of noise intensity (_{r} frequencies. There is a cut-off point for _{r} and driving frequency (_{d}). Within the multistability, there is further structure evidenced by the appearance of a staircase-like pattern, with sub-structures that repeat at approximately every first sub-harmonic of the resonance (_{r}). The multistability appears for all values of _{r} shown, indicating that this behavior could be present in different regions of the human brain known to exhibit a preferred resonance frequency (Galambos et al., _{d} in which multistability occurs extends further for systems with higher resonance frequencies (_{r}). The resonance frequencies chosen here are 4 (delta), 8 (theta), 12 (alpha), and 23 (beta) Hz; broadly chosen to represent the natural dominant frequencies found in the brain. To see the temporal evolution of

Effect of noise and driving frequency on connection strength. The effect of driving frequency (_{d}) and noise (_{r}) of 4 Hz

In order to explore the effect of synaptic noise on the dynamics of the synaptic strength, we analyzed the behavior of the connectivity strength for four fixed combinations of intrinsic resonance (_{r}) and driving frequency (_{d}), for _{r} and _{d} (chosen so as to highlight one of the multistable regions) that the multistable states exist for low levels of additive noise. In all cases, there is a critical value of the noise intensity after which there is a transition to a single broad region of attraction that tends to narrow down for higher levels of noise. Moreover, within the multistability regions shown in

Effect of noise on connection strength with fixed resonance and driving frequency. The effect of noise on multistable plasticity in the WC model. _{r} = 4 Hz, _{d} = 20 Hz; _{r} = 8 Hz, _{d} = 32 Hz; _{r} = 12 Hz, _{d} = 48 Hz; _{r} = 23 Hz, _{d} = 86 Hz.

There are two key aspects to explore further in relation to the results in

Effect of noise on the probability of connection strength with fixed resonance and driving frequency. _{r}/_{d} combinations of _{r} = 4 Hz, _{d} = 20 Hz, _{r} = 8 Hz, _{d} = 32 Hz, _{r} = 12 Hz, _{d} = 48 Hz, and _{r} = 23 Hz, _{d} = 86 Hz. It is apparent that for all cases, lower levels of noise result in higher mean connection strength and also the highest probability of attaining the high

To characterize the role of noise on the relative probability of the system choosing a particular synaptic strength state _{k}|_{r} and _{d} combinations described previously. We observe that for all cases, especially those cases with higher _{r}, which include _{c}), it is more probable that _{c}, the three states collapse into a single state which is close in value to the mid-state

In what follows, a generic excitatory synaptic weight _{ij} will be represented as

Let us take a time-average of the above differential equation over a coarse-grained time interval Δ >> τ_{h}. Thus, we have:

Where we used the fact that time averaging and differentiation commute, as shown below. For simplicity, we also assumed that during the time Δ the product of the activities was above the threshold. Let's notice that _{Δ} indicating the threshold correlation within the time-window Δ. Therefore,

We can now solve the equation within the time-window Δ. The solution is:

Where _{0} is the value of the synaptic strength at the start of the interval Δ. Since we assumed that Δ ≫ τ_{h} we can study the above solution for times of order Δ that satisfy Δ > _{h}. Thus, in that regime the exponentials are very small and can be neglected, leading to the final approximation,

Therefore, we conclude that the coarse-grained average of the synaptic strength is proportional to the correlation _{Δ}, with proportionality constant equal to γ. This result is consistent with Hebbian-like synaptic dynamics stating explicitly that synaptic strength between two units is directly proportional to the correlation between the activity of such units.

We now show that the time averaging and differentiation commute. Without any loss of generality, we assume that the time average can be written as a convolution with an appropriate kernel. That is,

Then, we take the derivative with respect to

Where the first term in Equation (16) is zero due to the kernel finite support and we used the prime as a short-hand notation of the time derivative of the kernel

As shown in Section 3.3, within the synaptic dynamics incorporated in our model a direct proportionality is expected on average between the correlation of activity between the sub-units' excitatory populations and synaptic strength. This is supported by the results shown in _{r} = 12 Hz, _{d} = 48 Hz, and

The effect of connection strength on correlation, synchrony, and phase-space dynamics. _{Δ}, where _{1} plotted against _{2}, and the time series for the fractional firing rates _{1} and _{2} (_{1} is blue, _{2} is red) for the three _{1} (blue) and _{2} (red) for the same sections of the simulation as in

To better understand these relationships, in

The relationship between the three stable states with both the correlation and synchrony suggests a nontrivial relationship in the phase-space dynamics of the excitatory populations. In _{1} and _{2}, in the “low” (left), “mid” (middle), and “high” (right), synaptic strength states respectively marked in blue in _{1} and _{2}. In particular, for most of the trajectory either variable remains almost fixed, while the other changes over a range of values. In the “mid” synaptic state, the trajectory in the phase space of both _{1} and _{2} detaches further from the axes leading to slightly higher correlations. Finally, on the right, we observe behavior consistent with the high correlation and phase locking values associated with the higher synaptic state: both _{1} and _{2} follow a cycle of high eccentricity aligned close to the diagonal in the _{1}, _{2} plane. In

In

Thus far, we have characterized synaptic multistability in a two-unit WC system. To explore how the results apply to larger systems we extended the model to include ten fully-connected units. The network was tuned to exhibit resonance at 12 Hz. Initially, we tested the effect of driving frequency on the final connection strength between each pair of units, _{d} between 20 and 80 Hz, chosen to highlight the region of multistable behavior. It can be seen that as was observed in the two-unit system, in the two central plots there are regions of multistability (highlighted in pink), although the behavior of the network is noisier. Next, we fixed _{d} = 48

Multi-unit simulation. Results for the 10-unit network are shown. _{r}) with a range of inputs (_{d}) between 20 and 80 Hz for the noise values z = 0 (left), _{r} = 12 Hz and _{d} = 48 Hz, chosen so as to highlight the multistable behavior for a range of additive noise (_{r} = 12 Hz and _{d} = 48 Hz, fixed as above. _{r} = 12 Hz and _{d} = 48 Hz as above.

Finally, we fixed the ten-unit network with parameters _{r} = 12 Hz, _{d} = 48 Hz,

Multi-unit simulation. The parameters for the 10-unit network were fixed as _{r} = 12 Hz, _{d} = 48 Hz, and

Each self-contained bubble (left and right) relates to a 1,000 ms time period indicated by the dashed lines. The connectivity matrices contain the final value for

The time-series for each of the 10 units within the selected time-frame are plotted color-coded to match the star diagram to illustrate the high coherence between units within each sub-network. Traces from the same sets appear so highly correlated that they completely overlap. The sub-networks switch and reorganize many times during the simulation and it can be seen that the connectivity architecture has been completely rearranged between these two time points. For instance, during the first period (left bubble), Units 1, 5, 6, and 10 (pink) are connected

Neural ensembles residing within disparate brain regions oscillate across a wide range of frequencies. To orchestrate complex brain functions, these sub-networks may transiently connect with each other across multiple spatial and temporal scales dependent on the task at hand (Fries,

Initially, we tested the effects of driving frequency and additive noise on the connection strength in a two-unit model. We found that for low levels of additive noise (

We next explored in more detail the effect of adding noise to models with a specific combination of resonance and driving frequencies chosen to highlight the multistable regions (

There has been significant computational work recently aiming to elucidate how multistable brain-states aid cognition. For example, Golos

A major aim of this work was to provide predictions that could be tested in plasticity experiments. In particular, a direct interpretation of our results indicate that it may be possible to optimize synaptic strength from the interplay between noise and driving frequency. With this in mind, we calculated the mean connection strength between the units as a function of noise for the same four combinations of resonance and driving frequency described previously and we found that for low values of noise (

To characterize the dynamics of the model for the different stable states of the connection strength, we compared the scaled truncated correlation γ_{Δ} (see Methods), to the coarse-grained connection strength, _{1} vs. _{2} shows a high degree of correlation (see

Finally, it is important to understand whether the behavior described so far for a small two-unit model extends to larger networks. For this, we implemented an all-all connected ten-unit network with fixed resonance of 12 Hz. We first repeated the investigation into the effect of driving frequency on connection strengths between the units and found that multistable regions of the connection strength did exist for specific ranges of driving frequency and additive noise, as was the case for the two-unit model (see

The findings presented here suggest that plasticity outcomes can be enhanced by specific levels of additive noise as a consequence of an interaction between noise and multistability. This prediction could be tested in multiple ways experimentally; the effect of extrinsic noise on plasticity and learning could be tested by controlling the intensity of noise added to a stimulus designed to elicit plasticity in a targeted brain region. As an example of this, Xie et al. (

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

CL-C and MM: conceptualization, data curation, formal analysis, investigation, methodology, visualization, writing—original draft, and writing—review and editing. LT: investigation, methodology, writing—original draft, and writing—review and editing. All authors contributed to the article and approved the submitted version.

CL-C was funded by the Medical Research Council (MRC) Grant MR/PO14445/1.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

The Supplementary Material for this article can be found online at: