Intensity- and frequency-specific effects of transcranial alternating current stimulation are explained by network dynamics

Objective. Transcranial alternating current stimulation (tACS) can be used to non-invasively entrain neural activity and thereby cause changes in local neural oscillatory power. Despite its increased use in cognitive and clinical neuroscience, the fundamental mechanisms of tACS are still not fully understood. Approach. We developed a computational neuronal network model of two-compartment pyramidal neurons (PY) and inhibitory interneurons, which mimic the local cortical circuits. We modeled tACS with electric field strengths that are achievable in human applications. We then simulated intrinsic network activity and measured neural entrainment to investigate how tACS modulates ongoing endogenous oscillations. Main results. The intensity-specific effects of tACS are non-linear. At low intensities (<0.3 mV mm−1), tACS desynchronizes neural firing relative to the endogenous oscillations. At higher intensities (>0.3 mV mm−1), neurons are entrained to the exogenous electric field. We then further explore the stimulation parameter space and find that the entrainment of ongoing cortical oscillations also depends on stimulation frequency by following an Arnold tongue. Moreover, neuronal networks can amplify the tACS-induced entrainment via synaptic coupling and network effects. Our model shows that PY are directly entrained by the exogenous electric field and drive the inhibitory neurons. Significance. The results presented in this study provide a mechanistic framework for understanding the intensity- and frequency-specific effects of oscillating electric fields on neuronal networks. This is crucial for rational parameter selection for tACS in cognitive studies and clinical applications.


Introduction
Transcranial alternating current stimulation (tACS) is a non-invasive neuromodulation method that can directly interact with brain oscillations [1].Besides an increasing number of basic research studies, tACS has recently been explored as a clinical tool for the treatment of psychiatric and neurological disorders [1].With its various adjustable parameters, in particular stimulation intensity and frequency, tACS is wellsuited for the development of personalized protocols.However, currently there is an incomplete mechanistic understanding of the intensity-and frequencyspecific effects of tACS at the neuronal and network levels.Specifically, the effects of stimulation intensity and frequency on different neuron types and their interactions have yet to be explored.
However, for weak stimulation intensities, neurons may also be de-entrained to the natural endogenous oscillation [14].The entrainment effects are typically frequency-dependent and increase with higher stimulation amplitudes [15].In the study by Huang et al. [15], tACS was applied to the ferret brain using a frequency range between ±4 Hz and the animal's endogenous alpha frequency.They found that the range of frequencies at which entrainment occurred for neurons expands as the strength of the stimulus increases, also known as Arnold tongues [15].
In humans, neurophysiological evidence of tACS effects is more indirect.Several studies have shown that tACS delivered in the alpha and beta frequency ranges can increase motor-evoked potential amplitudes, thus indicating an increase in corticospinal excitability of the primary motor cortex [16,17,18,19].The effects of tACS on cortical excitability as well as perception and cognition are often specific to a particular frequency range (e.g., theta, alpha, or beta) [20][21][22][23].Further, various studies have demonstrated frequency range-specific increases in EEG power [24,25].These observations are in line with the hypothesis that EEG frequency ranges are functionally distinct.Crucially, some human studies pointed toward frequency-specificity within a frequency range [26][27][28], although these findings are not demonstrated consistently [29].
Despite a wealth of existing literature on the effects of tACS across a wide range of stimulation parameters, theories that can integrate findings across studies are scarce.Computational models are crucial to developing a mechanistic understanding of experimental work.Additionally, modeling makes it possible to explore a wider range of stimulation parameters than is practical for experimental studies.Singlecell models demonstrated that pyramidal neurons (PY) are more sensitive to the exogenous electric field than other cell types due to their complex and elongated morphology [30,31].However, an isolated single PY cannot account for all tACS effects.To study the network-level effect of tACS, many modeling studies have utilized mean-field/neural mass models [5,[32][33][34][35], as well as spiking neuronal networks of single-compartment point neurons [15,36,37].These models demonstrated that an applied oscillating electric field causes network resonance in neuronal populations [5,[32][33][34][35][36][37][38].Thus, stimulation frequencies that are closer to the endogenous frequency are more effective at entraining the network.However, these models have limitations: the mean-field models do not account for individual variation in tACS effects on single neurons, and pointneuron models are unsuitable for direct electric field coupling due to the lack of distance over which the extracellular potential changes.Therefore, the twocompartment models in a spiking neuronal network [33,39] allow for a more realistic representation of the effects of electric fields.Neuron membranes can undergo depolarization and hyperpolarization at opposite ends, accounting for differences in depolarization caused by external electric fields.
Here, we developed a unique approach using a neuronal network with two-compartment model neurons and considered the orientation of neurons in the network to investigate the effect of tACS.We show that tACS entrains neuronal population dynamics via phase-locking.With this approach, we present two novel findings.Recent experimental findings suggest that tACS entrainment is intensity-dependent and is related to a shift in phase [14].Here, we aim to understand these findings using a computational approach.First, we demonstrate for the first time that tACS indeed induces a shift in the network's neurons' preferred phases.Second, we show that tACS causes desynchronization to the endogenous oscillation and resynchronization to the exogenous electric field as stimulation strength increases.Additionally, our modeling results suggest that excitatory and inhibitory synaptic couplings are important factors in tACS-induced network entrainment.Our network model provides a parsimonious account of known experimental results and provides important insights into the effects of intensity and frequency of tACS, which can inform future experimental studies.

Neuron models
We implemented conductance-based twocompartment neuron models consisting of soma and a dendrite (figure 1(A)).Single-compartment neurons are commonly used in network modeling, however, two-compartment models can capture the effects of electric fields [33] while still being computationally efficient.We used NetPyNE [40], a Python package to simulate and model biological neuronal networks based on the NEURON simulation environment [41].The membrane potential of the model neuron can be calculated by the following equation based on the cable theory [42]: where E || (x,t) is the parallel component of the electric field along the neuron compartment, V m is the membrane potential, r i is the cytoplasmic resistance, c m is the cell membrane capacitance, i ionic are the currents passing through membrane channels, x is the location of neuron compartment, and t is the time.The morphologies of each type of neuron are shown in figure 1(A).Large PY are often considered to be more susceptible to electric fields than small interneurons (IN) [30,31,43].Active currents in PY include a fast sodium current (I Na ), fast potassium current (I Kv ), slow non-inactivating potassium current (I km ), leak current (I L ), calcium current (I Ca ), and Calcium-dependent potassium current (I KCa ).The Hodgkin-Huxley-style kinetic equations of these currents are adapted from the model [44].Membrane electrical properties, such as membrane capacitance and ion channel conductance, were modified to generate a regular-spiking firing pattern.The IN model was adapted from the single-compartment model [45] and a dendritic compartment is added for electric field coupling.IN contain I Na and I Kv currents to create fast-spiking activity.The parameters used for the morphologies of the neurons are shown in table S1.

Single neuron model tuning
The reduced neuron model with simple morphologies helps speed up the simulation of the large network model.Neuron models with realistic morpho-logies can produce accurate predictions of the electrophysiological effect.We investigated the entrainment of reduced neuron and realistic neuron models under electric fields at 10 Hz, in the range of 0 mV mm −1 to 3 mV mm −1 .We scaled the length and the diameter of the soma and dendrite of the reduced model to match those of realistic neuron models [30].The electrical properties were adjusted accordingly (table S4).The firing rate of both neuron models was tuned to 10 Hz by systematically adjusting the weight of the synaptic input.Finally, we matched the entrainment of the realistic neuron by fine-tuning the electric field coupling of the reduced model.The reduced model's entrainment closely resembles that of the realistic model (figure S9).Although our reduced model does not approximate the complex morphology of the realistic neurons, given the high computational costs of network modeling, the use of reduced models is justified within the scope of the study, where we only focus on the electric field orientation along the somatodendritic axis of the neurons.

Network architecture and connectivity
The network consists of 1000 neurons with 800 PY and 200 IN [46].PY in the network were aligned to the vertical axis, and IN were distributed in 3D space (1000 µm × 1000 µm × 2430 µm) with random orientations (figure 1(A)) to reflect the typical cytoarchitecture of cortical columns.Neurons are connected to other neurons in-or outside their population.Excitatory and inhibitory synaptic connections within the network were constructed as shown in figure 1(A).Excitatory synapses connect to the dendrites, and inhibitory synaptic connections on PY and IN are located on the soma [47].All neurons in the network received background Poisson input to allow for spontaneous firing to mimic the ongoing in vivo activity of neurons.The rate of the Poisson process was set to 50 Hz.We used the NEURON class NetStim that generates an independent random Poisson spike train.We chose a value of 1 for the fractional randomness (0-1) which indicates that the spike train is completely random.
The connection probability was similar to what was chosen in Izhikevich [2006, 46], where the neurons were sparsely connected with 0.1 connection probability.Here, we slightly increased the connection probability of IN-PY to 0.15 and included selfinhibition of IN to generate the network oscillation [48].Each PY was connected to 100 random neurons, which could be PY or IN, each IN connected to 120 random PY and 20 random IN.Synaptic weight represents the NEURON connection weight of synaptic input to neuron models.It is associated with a change in synaptic conductance, governed by a dual exponential function (NEURON Exp2Syn).

Synaptic connections
The electric current I syn that results from a presynaptic spike was modeled as follows [49,50]: where the effect of transmitters binding to and opening of postsynaptic receptors is a conductance change (g syn ) in the postsynaptic membrane.V (t) denotes the transmembrane potential of the postsynaptic neuron and E syn is the reversal potential of the ion channels that mediates the synaptic current.If the neuronal transmembrane potential is below the reversal potential, presynaptic spike arrival leads to a depolarization of the neuron.The equation is used to describe the time course of the synaptic conductance [49,50]: where t s is the time of a presynaptic spike, τ 1 and τ 2 characterize the rise and decay time of the synaptic conductance in the dual exponential function.The conductance change resulting from α-amino-3hydroxy-5-isoxazolepropionic acid (AMPA), N-methyl-D-aspartate (NMDA), and γaminobutyric acid type A (GABA A ) mediated synaptic current are simulated using a dual exponential function (equation ( 3)).Parameters used for excitatory and inhibitory synaptic conductance are adapted from [51] and are shown in table S2.The synaptic conductance resulting from the background Poisson input is also modeled as a dual exponential decay with E syn = 0 mV, τ 1 = 2 ms and τ 2 = 10 ms [30].

Modeling cortical oscillations
To generate an ongoing, endogenous oscillation in the network, we adjusted the Poisson input to both the PY and IN, as well as the synaptic weights within and between the population.We sought to model a cortical oscillation, such as the alpha-band oscillations (8-13 Hz), which are dominant oscillations in the human brain [52], and have previously been effectively targeted in a number of tACS human studies [53][54][55].Additionally, to determine whether the same entrainment effects persist for various network dynamics, we modeled an endogenous oscillation at low beta frequency.We used the Genetic Algorithm to explore and find the best selected synaptic weights.Then we manually fine-tuned the parameters to achieve the final optimized values.We aimed to tune the network to obtain physiologically constrained mean firing rates across all neuron populations (approximately 10 Hz for PY and 55 Hz for IN) as well as to generate LFP signals with alpha (10 Hz) and low beta oscillations (14 Hz).The firing rate of PY and IN are in line with experimental studies [11,15,56,57].We implemented a reciprocal local connectivity scheme [58].The parameters used for the network connectivity are shown in table S3.
To add the heterogeneity to the network, we used a distinct seed of the Poisson input for each single neuron.

Local field potential (LFP)
The LFP signal was obtained by summing the extracellular potential induced by each segment of each neuron at the electrode location.Extracellular potentials were calculated using the line source approximation (LSA) [59,60] and assuming an Ohmic medium with conductivity sigma = 0.3 mS mm −1 [61].LSA implementation has been shown to better approximate extracellular signals than point source approximation and were previously validated on PY [61].The LFP of the network was computed and filtered using a 2nd order Butterworth bandpass filter (cutoff frequencies = 0.1 Hz and 300 Hz).

Modeling tACS
To couple the tACS electric fields to the neuron model, electric potential (Ve), was calculated for each neuron compartment based on the external electric field induced by tACS.When the electric field is assumed to be uniformly distributed, the electric potential equation can be simplified as follows [62]: where ⇀ E is the electric field vector, ⇀ s is the displacement.E x , E y and E z indicate the Cartesian components of the electric field in a three-dimensional space, and x, y, and z denote the Cartesian coordinates of each neuron compartment.The electric potentials were calculated at each neuron segment and were set as the extracellular potentials in the 'extracellular mechanism' in NEURON environment [63].Because an electric field is considered quasistatic, it can be divided into spatial and temporal components [64,65].
We used a sinusoidal tACS waveform and simulated a range of electric field strengths (0-1 mV mm −1 in increments of 0.1 mV mm −1 ) and frequencies (1-30 Hz in increments of 1 Hz).Using this range, we cover typical intensities and frequencies applied in humans [13] as well as explore the effect of tACS at a higher dosing range.

Phase-locking analysis
We quantified neural entrainment by calculating the phase-locking value (PLV), which measures spike timing synchronization with respect to an ongoing signal.The PLV is calculated as follows [66]: where N is the number of spikes and θ k denotes the phase of the oscillation at the time where the kth spike occurs.A uniform distribution of spike timings over all phases (0 to 2π) results in a PLV of 0, while a value of 1 means perfect synchrony to a specific phase of tACS.Under baseline condition, we measured PLV with respect to the LFP signal at a frequency of interest (filtered into a ±1 Hz range).While for tACS condition, PLV was calculated with respect to the tACS waveform [14].Note that for generating heat maps, cubic interpolation was used to smooth over data points.We quantified the phase preference of each model neuron by calculating the circular mean of the phases of the spike timing corresponding to the oscillation [67].

Evaluation of network effects
To evaluate the network effects of stimulation, the network size was decreased while maintaining the network's structure and neuron proportion, both in terms of number of neurons and number of synapses per neuron.As a result, we created a smaller network of 100 neurons with 20 synapses per neuron.The synaptic weights in the smaller network were adapted to obtain an endogenous oscillation with the same frequency (∼10 Hz) as the original network [68].

Model simulation and data analysis
The simulations were done in NEURON version 8.0 using a fixed-time-step cnexp method (modified Crank-Nicolson method) with a time increment dt = 0.05 ms.The total duration of an individual simulation was 4 min-a 2 min baseline period (no tACS) and a 2 min tACS period.All simulations were done using the Minnesota Supercomputing Institute (MSI).All post-simulation analyses were done in MATLAB 2021b.

Results
We connected 800 two-compartment PY and 200 fast-spiking interneurons (IN) in a 3D volume (figure 1(A)).To mimic the intrinsic network activity, we applied Poisson input to all neurons.We modeled tACS with the electric field vector aligned along the vertical axis of the network.We simulated network activity for 4 min: a 2 min baseline period (no tACS) and a 2 min tACS period.We measured the local field potential (LFP) from a virtual electrode inside the network (figure 1(A)).We first analyzed network activity at the baseline without external electric field.
The LFP signal during the 2 min baseline had a strong spectral peak at ∼10 Hz due to the endogenous network activity (figure 1 S1).
We then quantified the effect of 10 Hz tACS on neural entrainment in the network.Trajectories on the polar plots in figure 2(A) show the effects of increasing electric field strength on spike timing of both PY and IN. Figure 2(B) demonstrates a variability in PLV and phase preference of all neurons.In the baseline condition, all neurons were strongly entrained to the network's 10 Hz LFP component.For an electric field intensity of 0.3 mV mm −1 , the neural entrainment decreased relative to baseline.However, at higher strengths, the electric field reinstated and increased entrainment at a different preferred phase.This indicates that weak external electric field reduces entrainment to the endogenous oscillation, but at higher intensities, tACS overcomes and increases entrainment to the tACS waveform.Interestingly, the phase preference of neurons shifts as the exogenous electric field interacts with the endogenous oscillation.Given that the network response was dependent on the stimulation intensities, we undertook a secondary analysis by stimulating the network model with a broad range of electric field intensities (between 0.1 and 1 mV mm −1 , with 0.1 mV mm −1 increments) and stimulation frequencies (between 1 and 30 Hz, with 1 Hz increments).For each simulation, we computed the firing rate and PLV of all neurons and created heat maps to visualize the evolution of the firing rate and PLV across simulations.Firing rates of both PY and IN neurons were not affected by weak electric fields (<0.9 mV mm −1 ) (figure S2).At higher electric field strength (>0.9 mV mm −1 ), the change in firing rate was less than 5% (figure S2).Furthermore, the maximum electric field intensity of 1 mV mm −1 used in the simulation is subthreshold and insufficient to initiate action potentials.For entrainment effects, we found high-entrainment regions centered on the endogenous alpha frequency (∼10 Hz) and first harmonic (∼20 Hz) for both PY and IN neurons (figure 3(A)).This phenomenon cannot be explained by the small difference in somatic polarization as the stimulation frequency shifts from the frequency of the endogenous frequency.The somatic membrane polarizations across the study's stimulation frequencies (1 Hz-30 Hz) are similar (figure S10).The range of frequencies at which entrainment occurred expands as the strength of the stimulus increases.Interestingly, we found that both neuron types show high baseline entrainment at around alpha and the first harmonic frequency.Starting from the baseline levels of entrainment, neurons were desynchronized relative to baseline as low intensity tACS was applied.At higher electric field strength, neurons were resynchronized to the external electric field, resulting in entrainment within the high-synchronization region.To illustrate this behavior, we discuss four specific cases of a sample PY's entrainment within the stimulation parameter space (figure 3(B)).For 10 Hz tACS, as the electric field strength went up, the PLV first decreased from the baseline value and then increased.Polar histograms (figures 3(B-I)-(III)) demonstrate a shift in the phases of the spike timings.At the same electric field strength (1 mV mm −1 ), 10 Hz tACS entrained the neuron more effectively than 20 Hz tACS (figures 3(B-III) and (IV)).To test if the synchronization regions also appear for different endogenous frequencies, we repeated the simulations for another endogenous cortical oscillation.We tuned our network model to exhibit an endogenous oscillation around the low beta band (14 Hz) (figure S3).We observed similar high-synchronization regions on the PLV heat map (figure S4).This demonstrates that the same qualitative behavior arises at different endogenous frequencies.In the absence of endogenous activity, applying an electric field of 1 mV mm −1 and 10 Hz shows that the membrane potential solely reflects the tACS waveform with no entrainment effect (figure S8).These results suggest that cortical oscillations and tACS interact via phase locking at the single-neuron level within the network and network level.
Since single neurons within the network have a variable PLV and preferred phases under both baseline and tACS conditions across the population (figure 2(B)), we explored whether this spread of the PLV and phase preferences is due to the intrinsic firing rates of the neurons.We investigated the relationship between the neural entrainment of PY and their intrinsic firing rates under 10 Hz tACS, where we observed the largest tACS effects.Figure 4 shows that most PY with higher intrinsic firing rates have higher PLV under both baseline and tACS conditions.Notably, this trend plateaus at higher firing rates (above 10 spikes sec −1 ).Similarly, for phase preferences, PY with higher intrinsic firing rates are more likely to be entrained towards phases that are shifted backward (after peak) in time.Likewise, this trend plateaus after 10 spikes sec −1 .Under tACS, PY exhibit bigger phase shifts as electric field strength increases.Moreover, the logarithmic regression shows a high coefficient of correlation between PLV and firing rates (baseline condition: r = 0.859; tACS condition: 0.778, 0.822, and 0.892 for 0.3, 0.5, and 1 mV mm −1 respectively).In contrast, IN showed a different trend of PLV for varied intrinsic firing rates (figure S5).IN with high firing rates were less entrained than IN with lower firing rates.Because our model neurons are interconnected within the network via excitatory and inhibitory synaptic connections, entrainment is affected not only by the difference between neuron firing rates and stimulation frequencies, but also by synaptic coupling.
To study the mechanism of synaptic coupling on neural entrainment during tACS, we employed a simple network model with two neurons and one synaptic connection (figure 5).We investigated how the excitatory synapses (AMPA and NMDA) and inhibitory synapses (GABA A ) affect the PLV of the postsynaptic neuron.For unbiased comparison, we included a control condition (no synapse) in which the postsynaptic neuron receives no synaptic connection from the presynaptic neuron while maintaining the same firing rates by adjusting the Poisson input.The excitatory or inhibitory synaptic weights were adjusted to keep the firing rates of PY and IN at 10 and 55 spikes sec −1 , respectively.Figure 5 shows that the PLV of postsynaptic neuron was higher with an AMPA synapse than it was without synapse.Furthermore, the phase preference of the postsynaptic neuron was shifted further backward by NMDA synapses than by AMPA synapses (figure S6).GABA A synaptic connections, however, reduced the PLV of the postsynaptic neuron in comparison to the control condition (figure 5).Therefore, excitatory synapses enhanced the neural entrainment and caused phase shifts of the postsynaptic neuron.In contrast, inhibitory synapses reduced entrainment.Compared to the PLV of the isolated PY and IN, tACS effects on single neurons were amplified by synchronized synaptic inputs.As a result, excitatory and inhibitory synaptic coupling plays an important role in the network dynamics and thus contributes to the variability in single-neuron entrainment in the network (figure 4).
To understand how the entrainment of a single neuron enhances through network effects, we investigated the effects of tACS on PLV when the number of cells in the population changes.We used the same 10 Hz tACS protocol for a network of 100 neurons, which resulted in less entrainment than the original 1000-neuron network (figure S7).
Although both PY and IN were entrained by tACS at the frequency close to the endogenous frequency of the cortical oscillations (figure 3), we hypothesize that PY is the driving force of the network entrainment due to its higher sensitivity to the electric field.To determine whether the tACS-induced neural entrainment was driven by PY or IN, we performed a control analysis.Under 10 Hz tACS condition, we changed the electric field scaling factor (α) of each neuron type, which determines how neuron models respond to the external electric field.For both neuron types (PY and IN), we compared three conditions: (1) the neurons were not affected by electric field (α = 0), ( 2) the electric field coupling to neurons were doubled (α = 2), and (3) a default condition (α = 1).Due to IN's low sensitivity to the electric field, we added another condition with α IN = 10 to match the ratio of length between PY and IN.We computed PLV for all seven conditions (figure 6).By doubling the α of PY, the entrainment of both PY and IN increased strongly (figure 6(A)).When PY were not affected by electric field (α = 0), the PLV of all neurons was lowered to less than 0.05, indicating that the network was not entrained by tACS.On the other hand, by changing the α of IN, the entrainment of both neuron types changed minimally, suggesting that IN neurons only marginally contribute to the direct entrainment of the network.This shows that PY were directly entrained by tACS, whereas IN indirectly contribute to the entrainment effects of exogenous electric field through synaptic inputs from PY.

Discussion
We have developed a computational neuronal network that consists of coupled two-compartment excitatory and inhibitory neuron models.Using this network model, we investigated how the intensity and frequency of tACS influence its interaction with ongoing cortical oscillations.We observed a non-linear relationship with first de-entrainment followed by re-entrainment at a different phase.The phase shift reflects the transition from endogenous entrainment to exogenous entrainment and was seen in both PY and IN.Besides intensity, we show that entrainment of ongoing cortical oscillations also depends on frequency.When the tACS frequency matches the endogenous frequency more closely, the phase-locking between the ongoing endogenous oscillations and the exogenous electric field increases.Additionally, excitatory synaptic coupling and network effects amplify the tACS-induced entrainment.We found that PY were entrained directly by the electric field and further drive the indirect entrainment of IN.
We observed that tACS entrains single neurons at intensities that are feasible for human studies, without causing considerable change in firing rates (figure S2).This in line with previous in vivo nonhuman primate studies [11,12].As far as we are aware, we employed network modeling to demonstrate, for the first time, that the entrainment effects were non-linear.We observed de-entrainment of the network at low electric field strength and its subsequent re-entrainment at higher electric field.At intensities below 0.3 mV mm −1 , PLV decreased compared to baseline, suggesting desynchronization of neural firing in relation to endogenous oscillations.At intensities above 0.3 mV mm −1 , an increase in PLV was observed, which reflects synchronization of neural firing to the exogenous alternating current.In accordance, a recent meta-analysis suggested that AC stimulation can positively bias spike timing at intensities above 0.3-0.4mV mm −1 [13].Further, our findings match a recent in vivo study in non-human primates where it was shown that tACS competes with the ongoing neural activity in controlling spike timing of single neurons [14].Depending on the stimulation intensity, this competition between tACS and ongoing neural activity determines the increase or decrease in the rhythmicity of neuron firing and shift of the preferred phase of spiking within an oscillatory cycle.Our results also demonstrated the same phenomenon.Entrainment to the ongoing oscillations was disturbed by exogenous electric field at low intensity.The network was then synchronized to the exogenous electric field at a different preferred phase at higher intensities (figure 2).This phase shift can be interpreted as the transition of phase-locking from the ongoing oscillation to the external electric field.Further, tACS induced phase shifts have recently been demonstrated in humans at a network level [69].Although these findings are preliminary, it shows that tACS-induced entrainment is non-linear and can involve changes in phase preference.
We explored the parameter space of tACS and found that it interacted with ongoing cortical oscillations through phase-locking mechanisms, as demonstrated by resonance zones in the stimulation parameter space centered around the network endogenous frequency (figures 3(A) and S4B).Replicating this phenomenon, referred to as an Arnold tonguea synchronization region involving coupled oscillators and an external force [70]-we corroborated findings from previous modeling studies [5,15,[36][37][38].It suggests that tACS can entrain network oscillations more effectively when stimulating with the endogenous frequency of the network oscillations [4,5,9].A recent study has shown that tACS entrains alpha oscillations by following the Arnold tongue in both ferret experiments and a computational model of the thalamo-cortical system [15].Currently no direct evidence, such as single unit recordings, of the Arnold tongue phenomenon in humans has been reported.However, some indirect evidence suggests that effects of tACS are more pronounced when using individual peak frequencies, compared to standard values [26] or control frequencies slightly above or below the peak frequency [27,28].It should, however, be noted that such strong frequency dependency is not universally observed [29].
Variability in PLV and phase preferences among neurons was demonstrated during network entrainment in our model (figure 2).This observation prompted further investigation, focusing on the difference in intrinsic firing rates among individual neurons (figure 4).A computational modeling study [30] suggested that at the same stimulation frequency, an isolated single neuron that fires at higher rates relative to the stimulation frequency has lower PLV than the same neuron that fires at lower rates.This might explain why the PY trend plateaus after 10 spikes sec −1 during 10 Hz tACS and the PLV of IN shows an overall downward trend.Since our model neurons are coupled in the network via excitatory and inhibitory synaptic connections, the variability in PLV cannot be simply explained by the intrinsic firing rates of single neurons.Therefore, neural entrainment is affected by network dynamics.
Neuronal networks can amplify the tACS induced entrainment.Studies using in vitro and in vivo models provided evidence for the idea that active networks may be more sensitive to electric fields than isolated single cells [5,9,15].Our model also demonstrates that neuronal networks can amplify the tACS induced entrainment via excitationinhibition balance.Specifically, our model provided mechanistic insight: excitatory synaptic connections enhanced neural entrainment via AMPA and NMDA receptors, while inhibitory synaptic connections reduced entrainment at higher electric field strengths (figure 5).Recurrent feedback between excitatory and inhibitory neurons causes an increase in inhibition after excitation is increased.Such balance allows for generating stable periods of activity and modulating network oscillations [71,72].Previous research has suggested that tACS may target specific types of neurons: PY have been shown to be more susceptible to electric field while small interneurons are less susceptible [6,30,63].Other studies suggest that network interactions make tACS have especially strong effects on interneurons, which drives the remaining neurons [15].Despite being highly entrained, our modeling results indicate that interneurons were primarily entrained indirectly by the excitatory neurons (figure 6).This means that without connections from PY, IN do not highly synchronize to the exogenous current.Finally, the model implies that tACS effects are network dependent.Neural entrainment becomes weaker as network size and number of synapses decrease (figure S7).Entrainment relies on slight shifts in spike timings at single cell level that are amplified by network effects, such as number of neurons and synaptic connections affected by the weak exogenous electric field [5,32,73].
Our model, while capable of reproducing existing experimental findings and highlighting crucial tACS mechanisms, has some limitations.First, the network model does not include varied orientations and locations of PY.Although we only study tACS effects on a local network level, it can affect broad area in the cortex [69].Second, we did not consider various lengths of two-compartment PYs and INs.To generate more accurate estimates of the tACS effects while remaining computationally efficient, we tuned our cortical neuron models to match the PLV of those with realistic morphologies.However, different morphologies of neurons can potentially lead to a larger variability in their response to tACS [30,31,39].Third, we did not investigate tACS-induced modulation of neural plasticity mechanisms that could explain stimulation-related aftereffects [24,55,74].Recent study has shown that tACS-induced phase shifts are related to plasticity [69].Future work could incorporate spike-timing-dependent plasticity into the network and investigate the NMDA receptormediated plasticity during tACS.Lastly, the neuron density in our model is lower than experimental known results [75].Running long simulations with realistic cell density requires significantly longer computation time, which is not feasible with the current technical setup.Future work can explore this in more detail.
Our study focused on the mechanistic understanding of interactions between ongoing oscillations in a local network and exogenous electric field, yet it would be fascinating to study how tACS affects the generation of cortical oscillations by including exogenous driving inputs from other regions [15,47].In the network, we observed variability in the entrainment of neurons that fired at different rates (figure 4(A)).Future research could systematically examine how neuron firing rates affect tACS effects in both isolated single neuron models and coupled neuron populations.We studied the effects of network size on entrainment and demonstrated network amplification effects (figure S7).An intriguing area of research will be to compute the impact of other network properties including time constants and connectivity on tACS-induced entrainment.

Conclusions
This study used computational modeling to provide insight into the mechanisms of tACS effects on neuron populations.The findings highlight the importance of tailoring stimulation parameters to specific intensities and frequencies based on the endogenous neural oscillation.Different stimulation intensities would result in diverse neural entrain-ment effects.On the other hand, using a stimulation frequency that matches the intrinsic neural oscillation, phase-locking can occur at lower electric field strengths.Our model allows the exploration of tACS parameters and could pave the way for more effective use of tACS in clinical settings.

Figure 1 .
Figure 1.Network structure and endogenous alpha oscillation.(A) Top panel: a network of PY (n = 800, blue) and IN (n = 200, yellow) populations with excitatory and inhibitory synaptic connections within and between populations.Poisson input was applied to all neurons to generate ongoing network oscillation.Bottom panel: a schematic of the network in 3D space.PY are oriented along the vertical axis, while IN have random orientations to mimic the distribution of cortical neurons.LFP was measured in the extracellular space via a virtual recording electrode, shown in red.An alternating electric field (sinusoidal waveform) is applied to both PY and IN along the vertical axis.(B) Endogenous oscillation of the neuronal network.Top and middle panel: raster plot of neuron spiking activity and corresponding LFP signals for 1 simulated second.Bottom panel: power spectral density (PSD) of the LFP signals shows a peak in the alpha range (10 Hz).

Figure 2 .
Figure 2. 10 Hz tACS reduces entrainment and reinduces entrainment with phase shifts, as the stimulation intensity increases.Polar plots indicate the preferred phase (angle) and PLV (radial distance from the origin) of neurons in PY (blue) and IN (yellow).0 degree corresponds to the peak of the oscillation and 90 degree corresponds to the falling edge of the oscillation.(A) Each dot shows the average across neurons in each population.Trajectories for both PY and IN show reinstated entrainment (B) PLV and phase preferences of all neurons for four conditions in (A).Each dot represents an individual neuron.

Figure 3 .
Figure 3. Network response to tACS with frequencies between 1 and 30 Hz.The lower section of each heat map represents PLV at baseline condition when no stimulation was applied.The upper section of each heat map shows intensities between 0.1 and 1 mV mm −1 .The color at each point shows the average PLV across neurons.(A) Entrainment maps of neurons to tACS for each neuron type.Both PY and IN neurons show high PLV regions centered on the peak frequency (∼10 Hz).(B) Four points on PLV map are selected to show different degrees of entrainment effects of a sample PY in the population: (I) a point with high PLV under baseline condition.The sample PY is entrained to the endogenous 10 Hz LFP component.Rayleigh test.p-value = 1 × 10 −63 .(II-IV) Points with different stimulation parameters (II: 0.3 mV mm −1 and 10 Hz; III: 1 mV mm −1 and 10 Hz; IV: 1 mV mm −1 and 20 Hz.Rayleigh test.p-value = 1 × 10 −14 , 1 × 10 −123 , and 1 × 10 −26 , respectively).

Figure
FigureThe effects of tACS differ between PY in the network.10 Hz tACS was applied at three electric field strengths (0.3, 0.5, and 1 mV mm −1 ).(A) Phase-locking values across all PY are shown for baseline and three electric field strengths.Individual dots indicate values for each PY in the network.PY that fire at higher intrinsic firing rates (x axis) show greater PLV.(B) Preferred phases across all PY.At baseline, PY show phase preferences before 0 degree (peak).PY with higher intrinsic firing rates tend to be entrained after the peak.Under tACS, the same trend remained.With increasing electric field strength, PY show larger backward phase shifts.

Figure 5 .
Figure 5. Synaptic connections affect entrainment of the postsynaptic neuron.Two neurons were connected through a monosynaptic connection.10 Hz tACS was applied along the somato-dendritic axis of both PY (blue) and IN (yellow) neurons.PLV plots and corresponding polar plots (the radial axis shows PLV) and show the effects of synapse on the entrainment of postsynaptic neuron (with synapse) relative to control (no synapse).(A), (B) AMPA synaptic connection enhanced entrainment of postsynaptic neuron.(C) GABAA reduced entrainment of postsynaptic neuron.

Figure 6 .
Figure 6.The entrainment of the network is affected by the sensitivity of neurons to external electric field.α value of 0 indicates that the model neurons are not affected by tACS; 1 indicates that the electric field sensitivity of neurons is the same as the original network; and 2 indicates that the sensitivity is doubled.(A) The sensitivity of all PY to the tACS (αPY) was altered in the network.The averaged PLV of all PY and all IN are shown (standard deviation is plotted as the shaded area).(B) The sensitivity of all IN to the tACS (αIN) was altered in the network.Due to IN's low sensitivity to the electric field αIN was further increased by 10 times.10 Hz tACS was applied to the network in all conditions.