We built a network of 10,000 leaky integrate-and-fire (LIF) neurons, consisting of 8,000 excitatory (E) and 2,000 inhibitory (I) units, with a 10 % connection probability and plastic synapses, to represent a cortical network (see Materials and methods and Table 1). Under these parameters, the network exhibits a Synchronous-Irregular (SI) balanced state (Brunel, 2000), characterized by a power spectrum peaked in the upper β band, with an endogenous frequency f ∼ 30  Hz . We use f to denote the endogenous frequency, i.e., the frequency observed in the network in the absence of stimulation, and ω s to refer to the exogenous frequency, i.e., the stimulation frequency.

To promote entrainment and improve the signal-to-noise ratio (i.e., contrast between endogenous oscillations and tACS), we set the system in a weak-coupling regime. In this configuration, the ratio of synaptic input to stimulation amplitude remains comparable, especially during the early stages of the simulation, before plasticity significantly alters connectivity. Although individual synaptic weights are small in this regime (see Table 1), the net synaptic current amplitude is comparable to stimulation-induced fluctuations: the average maximum synaptic current during population synchronous spiking is approximately 0.5  mV , with a standard deviation of ∼ 0.1  mV (see Supplementary Material for more details). In contrast, the strong-coupling regime emerges after ∼ 600  s of spontaneous network activity in the absence of stimulation. During this period, synaptic plasticity modifies the connectivity such that synaptic input dominates over stimulation amplitude. This regime avoids competition between recurrent synaptic inputs and stimulation-induced fluctuations (Krause et al., 2022; Lefebvre et al., 2017). We explored this condition in the Supplementary Material and found qualitatively similar results. In the rest of this study, we focus on results obtained under the weak-coupling regime.

We subjected this network to periodic stimulation of various amplitudes ( A s ) (Schwab et al., 2019), and frequencies ( ω s ) , for a period of 15 seconds (simulation time, from t = 5 s to t = 20 s ). We then compared changes between the dynamics observed before stimulation (i.e., pre-stimulation) and after stimulation (i.e., post-stimulation) over epochs of 4 seconds. Specifically, we calculated the power spectrum over the pre-stimulation epoch (i.e., t = [ 1   5 ] s ), the stimulation epoch (i.e., t = [ 10   14 ] s ), as well as the post-stimulation epoch (i.e., t = [ 20.5   24.5 ] s ). The time intervals corresponding to each of those epochs have been selected to avoid any transient effects. To investigate the relationship between MTC heterogeneity and the persistence of stimulation-induced aftereffects, we plotted representative dynamics of the network in pre- and post-stimulation epochs in Figure 1 (See Supplementary Material for strong-coupling regime). We randomly selected 60 excitatory neurons and compared both network connectivity and the relative magnitude of synaptic weights between pre- and post-stimulation epochs in Figures 1A1,B1, respectively. Comparing these panels, one can readily notice stimulation-induced changes in synaptic weights and/or connectivity persisting well after stimulation offset. This effect was found to be mediated by variability in MTC. Corresponding synaptic weight matrices are plotted in Figures 1A2,B2, respectively. As shown in Figures 1A3,B3, the endogenous synchronous irregular activity present in the pre-stimulation period has been amplified in the post-stimulation epoch, accompanying a persistent increase in neuronal firing rates (Note that firing rates are lower than the network endogenous frequency as expected from irregular synchronous dynamics (Vogels and Abbott, 2005), i.e., the median of excitatory neurons firing rate is ∼ 0.5   H z for pre-stimulation and ∼ 1.5   H z for post-stimulation. See Figures 1A4,B4). The underlying population’s local field potential (LFP) (see Equation 5 in Materials and methods) also exhibits a significant increase in spectral power, especially salient at the endogenous (i.e., resonant) oscillation frequency and outlasting stimulation duration (see Figures 1A5,B5, also Figures 1A6,B6). We generalized these results in Supplementary Figures S2, S3, by choosing different distances to the threshold (by increasing the distance between resting and threshold potential) for each cell in the network, and by introducing heterogeneity in the threshold values. Note that this parameter change, even though it increased the amplitude of the oscillation, did not qualitatively change these results. The underlying mechanism behind this phenomenon may involve neurons being activated at different phases of stimulation, which induces selective synaptic weight modification and leads to amplified oscillatory activity. Further research is needed to explore the reasons behind this response, which are beyond the scope of this study.

Having identified post-stimulation amplification in endogenous oscillations, we next evaluated how this phenomenon depends on stimulation parameters. In Figures 2A,B, we plot the peak LFP power spectrum for various stimulation frequencies, both during and after stimulation offset. Stimulating at frequencies ranging from ω s = 1   H z to ω s = 40   H z ( A s = 1   ( m V ) ) invariably increases LFP power during entrainment, especially for stimulation frequencies near the resonant endogenous frequency. The effect carried over to the post-stimulation epoch: as can be seen in Figure 2B, peak power remained high around the population endogenous frequency despite no stimulation being present, indicative of stimulation-induced engagement of synaptic plasticity. Optimal post-stimulation peak power was observed at a stimulation frequency of ω s ∼ 23   H z , which we note is different from the network endogenous oscillation observed before stimulation onset ( f ∼ 28   H z ) . This indicates that stimulation-induced changes in synaptic coupling might be higher at non-resonant frequencies, which possibly reflects the interaction of the neurons’ MTCs with the stimulation frequency.

Stimulation amplitude is also crucial to elicit - and possibly maintain - persistent entrainment and associated changes in synaptic coupling. We plotted in Figures 2C,D the peak LFP power as a function of stimulation amplitude (i.e., A s ) both during and after stimulation offset. Two stimulation frequencies (i.e., ω s = 25 , and 30   H z ) were considered as they both reside within the range of frequencies for which the effect of post-stimulation LFP power is significant (see Figure 2B). While peak LFP power increases linearly with stimulation amplitude during stimulation epochs (see Figure 2C), a thresholding effect can be observed in the post-stimulation period. Indeed, a minimum stimulation amplitude appears to be required to cause post-stimulation LFP power amplification (Figure 2D). These results indicate that a high stimulation amplitude is required to modulate the neurons’ membrane potential and spiking response, to cause changes in connectivity significant enough to yield observed post-stimulation effects. The difference in LFP spectral power between the two selected stimulation frequencies (i.e., ω s = 25 and 30   H z ) indicates that, despite expected stimulation-induced resonance (here at ω s = f = 30   H z , see Figure 2D), amplification may occur at different, non-resonant stimulation frequencies. We however, emphasize that stimulation-induced change in synaptic coupling may trigger shifts in endogenous oscillatory activity, causing the peak power to fluctuate around a frequency of f ∼ 28   H z ( ± 1   H z std.).

We further investigated whether and how MTC heterogeneity is involved in generating those results. Is the LFP power amplification observed post-stimulation due to a global, non-specific increase in synaptic coupling, or is it instead due to selective, MTC-mediated synaptic plasticity? To answer this question, we first explored the effects of STDP on post-stimulation power amplification. As shown in Figure 2E, in the absence of stimulation (i.e., sham; A s = 0 ) while the network remains exposed to STDP due to its own endogenous activity, no significant shift in LFP power can be observed.

Stimulation-induced amplification in post-stimulation power was found to rely heavily on selective synaptic modifications, i.e., synapse-specific directional changes resulting from periodic entrainment of neurons possessing distinct MTCs (Pariz et al., 2023). To expose the role of such selectivity, we randomly shuffled synaptic weights amongst neurons of the same cell-type while preserving their overall statistics (see Materials and methods). Figure 2F compares the spectral power obtained without stimulation (sham control; A s = 0   ( m V ) ) and post-stimulation (stim. Control; ω s = 25   H z , A s = 1   ( m V ) ) conditions with those obtained by shuffling and/or sampling synaptic weights randomly while preserving their respective distributions, within and between cell types. To do this, we first calculated the synaptic weight distribution amongst all synaptic types (i.e., E → E , E → I , and I → E ; Note that I → I remained unchanged). We next randomly shuffled synaptic weights in the network and examined whether post-stimulation oscillatory amplification could be observed over epochs of 4 s (no stimulation was applied during that period). As shown in Figure 2F, no post-stimulation increase in power could be observed, indicating that while displaying the same overall statistics (i.e., being shuffled, there are no changes in synaptic weights value and the distribution remains unchanged within each cell-type; See Materials and methods), selective plasticity between neurons with distinct MTCs is essential in generating amplified oscillation. We pushed the analysis further and sampled synaptic weights independently, only using the cell-type specific distributions calculated above (i.e., agnostic of the actual values of those weights). With this, the same result could be observed: in the absence of selectivity, post-stimulation oscillatory amplification vanishes.

Our results indicate that despite the significance of oscillatory amplification and its manifest reliance on MTC heterogeneity, all reported post-stimulation after-effects were found to be transient, as reported in several studies (Zaehle et al., 2010; Kasten et al., 2016; Vossen et al., 2015) and dissipate over time after stimulation is turned off. Upon stimulation offset, prevailing endogenous synchronous irregular activity engages STDP to bring synaptic connectivity back to baseline (see Figure 2G).

We examined the evolution of synaptic weights between all types of synapses in Figure 3 with respect to differences in MTCs, i.e., Δ τ m = τ m pre − τ m post . In Figure 3, we plot synaptic weights evolution for different stimulation frequencies i.e., ω s = 15   H z (Figures 3A1–A4), 25   H z (Figures 3B1–B4), and 35   H z (Figures 3C1–C4). These frequencies were selected to help the comparison between the dynamics and resulting plasticity at stimulation frequencies that either amplify the post-stimulation power (i.e., ω s = 25   H z ) and frequencies that do not ( ω s =   15 ,   35   H z ; see Figure 2B). Although synaptic changes are noticeable in all of these cases, their relative magnitude was found to be highly frequency-specific. For instance, synaptic weights between excitatory and inhibitory neurons (i.e., E → I (Figure 3A2,B2,C2), display a broader range of synaptic modifications at ω s = 25   H z (Figure 3B2) compared to other frequencies (Figure 3A2,C2). This indicates that the stimulation frequency ( ∼ 25   H z ) , solicits MTC heterogeneity more strongly, leading to selective synaptic changes spanning a greater range of Δ τ m and stronger power amplification (see Figure 3A4,B4,C4). This is in contrast to Figure 3A2,C2 where synaptic weight changes were more selective for negative Δ τ m . The same effect could be observed for synapses between different cell types: selective modification observed amongst E → E and I → E synapses displayed a similar trend. Synaptic weight changes observed scaled with MTC mismatch as previously reported (Pariz et al., 2023), and further persisted over time after stimulation offset. In the last column, (A4), (B4), and (C4), for comparison purposes, we plotted the pre- and post-stimulation power resulting from each stimulation frequency used ( ω s =   15 ,   35   H z ) .

Heterogeneity amongst and between different cell types, either excitatory or inhibitory, has different consequences on the post-stimulation power. To quantify this, we explored in Figure 4 the effects of cell-type MTC heterogeneity on post-stimulation LFP power. As shown in Figure 4A, MTC heterogeneity among excitatory neurons (i.e., σ τ m E , along the horizontal axis) enhances post-stimulation power (i.e., S max ), whereas increasing MTC heterogeneity among inhibitory neurons (i.e., σ τ m I , along the vertical axis) abolishes the effects (See Figure 4A). The greater diversity observed among cortical inhibitory interneurons compared to excitatory neurons (Soltesz, 2006), may hinder stimulation effects and possibly prevent power amplification. However, it should be noted that the frequency of stimulation is another factor that determines the stimulation effects. We measured this in Figure 4B, where we varied the level of MTC heterogeneity of E and I neurons (i.e., both E and I neurons were assumed to express the same variation in MTC heterogeneity σ τ m ) and the frequency of stimulation. A similar increase in MTC variability of E and I neurons contributes to the induction of post-stimulation effects over a wider stimulation frequency range, i.e., [ 20   30 ]   H z . Having the same heterogeneity among inhibitory and excitatory neurons amplifies response power and therefore creates the necessary conditions for optimal synaptic weight changes, which ultimately leads to the amplification of oscillation power.

These results highlight the importance of considering plasticity among and between neuron subtypes. To investigate which synapses are more significantly involved in mediating the post-stimulation aftereffects, we applied periodic electrical stimulation on the same population at different degrees of MTC heterogeneity while selectively turning ON and OFF STDP amongst different cell types. This enabled the identification of synapses whose plasticity is more significantly solicited during stimulation. In Figures 4C1–C4, we show that plasticity between excitatory to excitatory neurons (that is, E → E ) and between inhibitory to excitatory neurons (that is, I → E ) is more involved in the amplification of the LFP power. The effect was also found to scale with the level of MTC heterogeneity across cell types (excitatory and inhibitory neurons) as of Figures 4C1–C4 where the post-stimulation power amplified as we increased the σ τ m E ,   I . Introducing plasticity among inhibitory neurons, under the same conditions as previously considered, is found to suppress the amplitude of post-stimulation aftereffects (see and compare Supplementary Figures S1A, 2AB). These results suggest that blocking synaptic plasticity, whenever applicable, among synaptic subtypes may lead to a significant increase in post-stimulation power. Further investigations are required to determine the implications of MTC and cell-type specific synaptic blocking on tACS-induced aftereffects.

We modelled a population of excitatory and inhibitory Leaky-Integrate and Fire (LIF) neurons (Brette, 2015; Tuckwell, 2006). The differential equation for the evolution of the subthreshold membrane potential of each neuron is τ m d v d t = V rest − v + I ζ + I syn + I s , (1) where τ m is the MTC, v is the membrane potential, V rest is the resting membrane potential, and I ζ represents an external current modelled here as white noise with a mean value of μ and a standard deviation σ . The term I syn represents the synaptic current, while I s is the stimulation-induced current, which is here assumed to be a sinusoidal input (representing transcranial alternating current stimulation (tACS), (Krause et al., 2019; Krause et al., 2022; Schwab et al., 2021)), i.e., I s = A s ⁡ sin ( 2 π ω s t ) , where A s and ω s are the amplitude of the periodic signal, and the angular frequency respectively. When a neuron crosses the threshold value v thr = − 54   ( m V ) , it spikes and its membrane potential resets to resting value V rest = − 60   ( m V ) and remains there for τ ref = 2   m s representing the neuronal refractory period. Although having larger refractory periods alters the neurons’ firing, the results remain consistent (not shown). The parameters are in the physiological range (Gerstner et al., 2014) and summarized in Table 1. The total simulation time, unless otherwise stated, is 30 s, including pre-stimulation (sham epoch): t ∈ [ 0   5 ) s , stimulation epoch: t ∈ [ 5   20 ) s , and post-stimulation epoch: t ∈ [ 20   30 ] s (Extended stimulation periods did not show any significant difference. Data is not shown). The total synaptic current for neuron i is given by I syn i = ∑ j = 1 N E g i j E S i j t v i − E syn j + ∑ j = 1 N I g i j I S i j t v i − E syn j where g i j E , I are synaptic weights matrices associated with connections between either excitatory (E) and inhibitory (I) presynaptic neurons towards a postsynaptic neuron i . The sum is taken over N E excitatory and N I inhibitory presynaptic neurons over two nearest spike times. The reversal potential, E syn , for E and I neurons are 0   ( m V ) and − 80   ( m V ) , respectively. The synaptic response function S i j ( t ) for connections from neuron j to neuron i is modeled as S i j t = Λ   exp − t − t s p j − t d i j τ r − exp − t − t s p j − t d i j τ d Λ = 1 / τ r τ d τ r τ d − τ r − τ r τ d τ d τ d − τ r where t s p j is the spiking time of j t h neuron, and t d i j is the axonal delay between presynaptic neuron, j , and postsynaptic neuron, i . The τ r and τ d , are rise and decay synaptic time constants, respectively, associated with GABA a and AMPA receptors (see Table 1); (Gerstner et al., 2014).

Individual synaptic weights are randomly chosen from a normal distribution with mean and standard deviation as given in Table 1 (sham control). In shuffled (see Figure 2F), first we let the simulation run for 20 s , and instantaneously shuffled the synaptic weights at the beginning of the post-stimulation epoch. We shuffled synaptic weights within each synapse category (i.e., the synaptic weights among E and I neurons). In the sampled case (see Figure 2F), we took the following procedure: We let the population in stim. Control evolve for 20s (5s pre-stimulation, and 15s stimulation epochs). We then calculated the distribution of the synaptic weights at the end of the stimulation epoch. Then we used these distributions to randomly sample synaptic weights within each synapse category ( E → E , E → I , and I → E ) using this fitted distribution. To fit the distribution, we used cftool package in MATLAB. We seek any function that fits the data with R − s q u r e > 0.95 % . A representation of this distribution is being shown in S4. These tests demonstrate that while the overall distribution of synaptic weight may remain intact through shuffling or sampling, selective modification is essential for inducing post-stimulation aftereffects.

Plasticity in our population amongst connected neurons is modelled using Hebbian pair-based spike-timing dependent plasticity (Choe et al., 2013; Gütig et al., 2003; Sjöström et al., 2010). To avoid biased synaptic changes (i.e., preferential LTP/LTD.) we chose a symmetric STDP Hebbian learning rule (Bi and Mm, 2001; Gütig et al., 2003; Sjöström et al., 2010). The synaptic weight dynamics in our model follows the below equations: Δ g = A + 1 − g / g max exp − Δ T / γ + , if  Δ T ≥ 0 . − A − g / g 0 exp Δ T / γ − , if  Δ T < 0 . g = g + Δ g (2)

The γ + and γ − are STDP decay time constants. Δ T = t s p post − t s p pre − t d is the time difference between the spiking time of post- and presynaptic neurons, and t d is delay between presynaptic and postsynaptic neurons. Whenever Δ T is positive (negative), the synaptic weight between presynaptic to postsynaptic neurons gets potentiated (depressed). The constant g max denotes the maximum achievable synaptic weight, while g 0 denotes the initial synaptic weight, taken from a narrow Gaussian distribution across all synaptic connections before learning (see Table 1).

Baseline synaptic coupling and threshold were selected to set the network in a weak-coupling regime, sub-threshold regime, in which an isolated presynaptic spike does not guarantee postsynaptic firing. This regime achieved by choosing the synaptic weight from a narrow distribution (see Table. 1) and at the early stage of simulation. Despite weak synaptic coupling, the afferent synchronous synaptic input each neuron receives from the rest of the network is comparable to the stimulation amplitude, i.e., the average of maximum synaptic input (at the onset of every synchronous spiking activity) and its standard deviation is ∼ 0.5   ( m V ) and ∼ 0.1   ( m V ) , respectively. In the strong-coupling regime, which the network reaches after 600 s of simulation time (in the absence of stimulation), the average maximum synaptic input and its standard deviation reach approximately ∼ 1.5   ( m V ) and ∼ 0.45   ( m V ) , respectively. Throughout this report, we used Equation 2 for synaptic modification, and our choice of STDP parameters are A + = 2 A − = 4 × 1 0 − 4 , g max = 2 g 0 and γ ± = 10   m s .

We modelled a randomly connected sparse network of 10,000, LIF neurons (see Equation 1) with a 4:1 ratio of E (8000) and I (2000) neurons with a fixed connection probability of 0.1 (Vogels and Abbott, 2005; Bryson et al., 2021; Campagnola et al., 2022). To balance physiological relevance and computational tractability for the network sizes we used the LIF neurons model (Burkitt, 2006). The synaptic weights and other neurons’ parameters have been selected within the reported physiological range (Campagnola et al., 2022) to be in line with previous studies on LIF cortical network models (see (Burkitt, 2006; Kobay and ashi, 2009; Brunel and Wang, 2003) and references therein), and are further summarized in Table 1. To study the effect of MTC heterogeneity, we randomly sampled neuronal MTCs ( τ m ) from Gaussian distribution with μ τ m = 10   m s , and σ τ m = 3   m s unless otherwise specified. The resulting population exhibits a synchronous irregular activity (SI) (i.e., see Figure 1A3).

To perform spectral analysis of the network’s mean activity, we first calculated the local field potential (LFP), V ̄ as the weighted ensemble average of the membrane potential i.e., (Herrmann et al., 2016; Hutt et al., 2018; Lefebvre et al., 2017; Mazzoni et al., 2015; Bazhenov et al., 2001), V ̄ = 0.8 N E ∑ i = 1 N E V E i + 0.2 N I ∑ i = 1 N I V I i , (3) where the relative proportion of excitatory (0.8) versus inhibitory interneurons (0.2) cells is taken into consideration. The power spectral density of V ̄ was averaged over 10 independent trials in which the same stimulation protocol is applied, but using different baseline connectivity, synaptic weights, and noise realizations. For the purpose of Figures 1, 2, 4, we further took the average of the power spectral density with a moving average window (with MATLAB smooth function) with σ = 1.5   H z , that provided us with a smoothed power-frequency curves (for instance, see Figure 2B6).