A biophysical model of dynamic balancing of excitation and inhibition in fast oscillatory large-scale networks

Over long timescales, neuronal dynamics can be robust to quite large perturbations, such as changes in white matter connectivity and grey matter structure through processes including learning, aging, development and certain disease processes. One possible explanation is that robust dynamics are facilitated by homeostatic mechanisms that can dynamically rebalance brain networks. In this study, we simulate a cortical brain network using the Wilson-Cowan neural mass model with conduction delays and noise, and use inhibitory synaptic plasticity (ISP) to dynamically achieve a spatially local balance between excitation and inhibition. Using MEG data from 55 subjects we find that ISP enables us to simultaneously achieve high correlation with multiple measures of functional connectivity, including amplitude envelope correlation and phase locking. Further, we find that ISP successfully achieves local E/I balance, and can consistently predict the functional connectivity computed from real MEG data, for a much wider range of model parameters than is possible with a model without ISP.


Defining ISP convergence
By definition, plasticity ceases when !" #$ % !% is zero. Hence, the convergence or success of ISP would normally be quantified by comparing the target activity rate to the instantaneous excitatory activity: = However, in this study we will be operating the model within an oscillatory regime. As a result, !" #$ % !% will vary over the period of an oscillation, such that on electrophysiological timescales, the system oscillates between being above and below the target activity level.
Instead of assuming the system has converged when !" #$ % !% is zero, we consider the system to have converged if ,-is zero integrated over a period spanning several oscillations, which implies there is no sustained change in ,-over time. Integrating both sides of Eq. (4) yields ,/0 ,-= − , A1 and setting ,-= 0 provides the requirement for convergence as

= . A2
Note that Eq. (A2) corresponds to a weighted average of excitatory activity. A further subtlety to consider is what happens when the system is operating at a point where it produces sustained nonlinear oscillations. Since inhibitory activity in the model is driven by local excitatory activity, is largest when is also large. The weighted average in (A2) therefore weights times with high excitatory activity more strongly, and an unweighted average of excitatory activity over this same time would be smaller than the expression in (A2). In these cases, the system will thus converge to a state where that the average excitatory activity is smaller than .

Robustness to noise
To demonstrate the robustness of dynamics and functional connectivity to noise, Figure A2 shows simulations using the same delay/coupling parameters as used in Fig. 5, but with different levels of noise. Note that we included the ISP stage of the simulation as well, to verify that ISP is also robust to the choice of noise signal. For these simulations, the noise signals were sampled from Gaussian distributions with zero mean and standard deviations of 0, 0.05, 0.1, and 0.2, which correspond to no noise, and to 5, 10, and 20 times the original amount of noise. The results are qualitatively very similar even with a 10-fold increase in noise, although larger differences begin to emerge after this point. The phase-based connectivity metrics are more sensitive to noise than AEC, but Fig. A2 shows that the magnitude of the phase connectivity metrics is still highly robust to a moderate increase in noise. This provides confirmation that our results primarily reflect the nonlinear dynamics of the model, rather than the noise process.

Synchrony with orthogonalisation
As synchrony and metastability are affected by spatial leakage, it is necessary to include an orthogonalisation step when comparing them to the model. The synchrony and metastability matrices shown in Figure A3 are from the same simulations shown in Fig. 6c and 6d, except that the parcel activity timecourses were orthogonalised prior to estimating the phase timecourses. Experimental estimates of synchrony and metastability should be compared to

ISP convergence in brain simulations
To investigate the convergence of ISP, we examine the standard deviation of ,-, after providing enough time for the system to converge if possible. Figure A5a shows the standard deviation of ,-computed over the final 500 seconds of the simulation (discarding the first 1000 seconds), averaged over brain regions. We consider three regimes -a highly synchronous regime with low delay and strong coupling, a metastable regime that gives rise to realistic functional connectivity, and a low synchrony regime with weak coupling and long delays. The mean activity and plasticity timecourses for example parameters in each regime are shown in Figs. A5c-A5h.
We first examine the discrepancy in activity that gives rise to ongoing plasticity. Fig. A5b shows the mean I-weighted E activity (Eq. A1) in each brain region, sorted by value. For the low synchrony and metastable parameter regimes, all brain regions have activity close to the target value. However, for the high synchrony parameters, although a small number of brain regions can reach the target level of activity, some exhibit higher levels of activity, while others exhibit lower levels of activity. The corresponding activity time series are shown in ongoing decrease in inhibition from 500s onwards), which indicates that although inhibition in that brain region is decreasing, the amount of activity is not increasing.
Our interpretation of these results is that the high synchrony regime exhibits dynamics that are strongly dependent on the long-range connectivity and the network, resulting in the asymmetries in activity being robust to moderate changes in synaptic strength in parts of the network. In contrast, for the weakly coupled network it is possible to optimize the mean level of activity in each brain region independently, due to the weaker effect of the rest of the network. This enables each region to independently achieve a local balance between excitation and inhibition. The metastable regime shown in Fig. A5g lies somewhere in between. On short timescales (on the order of 0.5-2 seconds), there are asymmetries in activity, but these are not stable and continuously change. The network can be considered to transiently visit different asymmetric oscillatory states. However, on long timescales that average over these states, the metastable regime can achieve the target level of activity in all brain regions, as shown in Fig. A5b. We note that because convergence in the metastable case requires averaging over transient states, it is especially important to ensure that the plasticity timescale ,/0 is slow enough that the synaptic strengths are approximately constant for the duration of each transient state, which is considerably longer than the oscillatory period of the units (~0.1 seconds).  Next, we investigated the convergence properties of ISP when changing the target activity level . Figs. A6a-A6c show the variability in ,-over the final 500 seconds of the simulations, for the low, medium, and high activity targets. For the low target activity = 0.1, ISP can achieve the target firing rate regardless of the long-range delay and coupling.
For the high activity, ISP target = 0.3, variability in ,-is small for the parameter regime that best matches experimental data. However, ISP no longer converges in the weakly interacting regime, with long delays and low couplings. Fig. A6d shows the timecourse of ,for the parameters marked by the red dot in Fig. A6c. There are large ongoing oscillations in synaptic strength. The amplitude of these oscillations does not change when the plasticity rate is decreased (at 500s and 1000s), which indicates that it is not the result of coupling between neural activity and plasticity. Figs. A6e and A6f show a detailed segment of the fluctuation in ,-and the corresponding excitatory activity. For these parameters, high inhibition is associated with low activity, low inhibition with high activity, and there is a bistable regime where both the low activity and high activity dynamics are possible with the same level of inhibition. Initially, the mean activity is below the ISP target in all populations, which causes inhibition to decrease. However, a bifurcation is encountered at which point the amplitude of the oscillations in the network suddenly increases, and the mean activity becomes much higher than the ISP target. In response, inhibition begins to increase again. However, the high amplitude oscillation does not lose stability immediately and when the network transitions back to a low activity state, the mean activity is well below the ISP target. Plasticity thus moves the system periodically over the bistable regime, reversing direction whenever one of the activity patterns becomes unstable, because the ISP target lies between the high and low activity regimes.
In summary, we have identified two distinct cases where ISP fails to converge -a case where strong global interactions prevent all of the units from achieving the same level of mean activity, and a case where a periodic oscillation in synaptic strength occurs if the ISP target lies between bistable high and low activity states. For the parameters that we tested in this study, convergence was not a problem in the parameter regimes that give rise to realistic functional connectivity, but there are clearly restrictions on the level of synchronization and on the ISP target if the synaptic strengths are to converge.
It is a possibility that the periodic changes in synaptic strength are only a local minimum solution, and that a different initialization for the synaptic strengths could result in convergence of ISP. Due to the computationally intensive nature of the ISP simulations, a practical approach for testing this would be to use an offline optimization mechanism to first check for the existence of a stable solution for plasticity with the desired activity level, as this is a necessary requirement for ISP convergence. Even this approach is challenging -the parameter space is high-dimensional with as many synaptic strengths to optimize as brain regions, and we have observed that in the optimal parameter regime, the synaptic strengths are strongly correlated but not precisely correlated with long-range excitation (Fig 7b), which means that the problem cannot be reduced to a two-parameter optimization of the line shown in in Fig 7b. The cost function would also be extremely expensive compared to those used for similar previous work (Deco et al., 2014) because the transient states in the optimal parameter regime require that simulations be run for a sufficiently long time to enable averaging over the states, and delays need to be included. Assessing the stability of a plasticity solution, if one is found, is also challenging because whether plasticity causes perturbations in synaptic strength to grow or decay depends on the highly nonlinear fast dynamics in the model. Regardless, offline optimization algorithms will also need to account for the possibility that there may be no stable solutions. In contrast, the online learning rule we have investigated in this study is simpler to investigate and more biologically relevant, with the trade-off that we do not exhaustively explore the space of synaptic strengths.