Effects of transcranial alternating current stimulation on spiking activity in computational models of single neocortical neurons.

Neural oscillations are a key mechanism for information transfer in brain circuits. Rhythmic fluctuations of local field potentials control spike timing through cyclic membrane de- and hyperpolarization. Transcranial alternating current stimulation (tACS) is a non-invasive neuromodulation method which can directly interact with brain oscillatory activity by imposing an oscillating electric field on neurons. Despite its increasing use, the basic mechanisms of tACS are still not fully understood. Here, we investigate in a computational study the effects of tACS on morphologically realistic neurons with ongoing spiking activity. We characterize the membrane polarization as a function of electric field strength and subsequent effects on spiking activity in a set of 25 neurons from different neocortical layers. We find that tACS does not affect the firing rate of investigated neurons for electric field strengths applicable to human studies. However, we find that the applied electric fields entrain the spiking activity of large pyramidal neurons and large basket neurons at < 1 mV/mm field strengths. Our model results are in line with recent experimental studies and can provide a mechanistic framework to understand the effects of oscillating electric fields on single neuron activity. They highlight the importance of neuron morphology and biophysics in responsiveness to electrical stimulation.

Introduction

Variations of local field potentials (LFPs) take their origin from brain oscillations which play a key role for brain communication. These oscillations are known to affect spiking activity by means of de- or hyperpolarization of neural membranes (Draguhn and Buzsáki, 2004). An increasing number of psychiatric diseases (e.g. schizophrenia) are being understood as deviations from healthy neural oscillatory activity (Uhlhaas and Singer, 2006). Thus, tools to directly interact with brain oscillations hold promise for new therapeutic approaches to harmonize pathological brain rhythms. Consequently, there has been a rising interest in non-invasive brain stimulation methods to modulate brain oscillations in humans. One popular method, transcranial alternating current stimulation (tACS), works by passing a weak oscillating electric current through electrodes attached to the scalp, creating an electric field inside the brain (Antal and Paulus, 2013). This electric field interacts with the spiking activity of neurons through transient membrane de- or hyperpolarization at the applied stimulation frequency.

Many experimental studies have investigated the effects of oscillating electric fields either in-vitro (Chan, Hounsgaard, 1988; Deans et al., 2007) or in-vivo animal studies (Fröhlich and McCormick, 2010; Johnson et al., 2020; Krause et al., 2019; Simal S. Ozen et al., 2010). One of the main hypotheses for tACS is the idea of neural entrainment. Entrainment means that neural activity at either the single neuron or population level synchronizes to a certain extent to the applied stimulation waveform (Chan and Nicholson, 1986; S. Simal Ozen et al., 2010). However, entrainment is only one specific form of more general changes in spike-timing, meaning a deviation from regular spiking activity due to the applied stimulation. Thus, despite a growing body of experimental studies, the exact mechanisms involved in the alteration of neural activity due to tACS are still not fully understood.

Computational modeling is a key tool, complementary to experimental work by providing a theoretical framework in which experimental findings can be understood. Modeling further allows to efficiently explore a larger range of stimulation parameters than feasible in experimental studies. Several computational studies have explored the effects of tACS on single neurons or neural populations (Ali et al., 2013; Cakan and Obermayer, 2019). These studies have highlighted the idea of network resonance in neural populations due to the applied oscillating electric field. Common to these studies is the use of simplified neuron models such as the ball-and-stick model (Aspart et al., 2016) or the two-compartment model (Ladenbauer and Obermayer, 2019). While such models are well suited for modeling large neuronal populations, they can give little insight about single neuron responses. This is due to the lack of a detailed dendritic morphology which strongly shapes neural activity (Ostojic et al., 2015). Further, differences in neural morphology and biophysics cannot be captured by simplified models but require detailed individualized neuron models. Thus, to study the effects of tACS on the single neuron level in different cell types, it is necessary to develop morphologically realistic neuron models with detailed biophysics.

Here, we conduct a computational study to investigate the effect of tACS on single-neuron spiking activity in a set of morphologically realistic neocortical neurons, using cable theory to accurately model the underlying biophysics (Rall, 1962). Ongoing spiking activity was created through synaptic input randomly located in the dendritic branches. We explore the effect of electric field strength on neural activity at 10 Hz stimulation frequency and investigate neural entrainment and changes in spike timing and firing rate. Our computational study can thus give important insights into tACS mechanisms.

Methods

Overview

tACS generates an oscillating electric field in the brain through stimulation electrodes attached to the scalp. Through the choice of an electrode montage (Fig. 1A), one can target a specific brain area such as the motor cortex (Fig. 1B). Stimulation frequency and amplitude are parameters that can be experimentally controlled (Fig. 1C). When a neuron is subjected to an electric field, the branches at one end of the neuron become hyperpolarized while the branches at the other end become depolarized. This de- and hyperpolarization follows the temporal pattern of the stimulation (Fig. 1D). These rhythmic changes in polarization affect the excitability of the neuron and timing of action potentials. Here, we develop computational models to study the effect of tACS on spike timing in detailed neuron models.

Fig. 1.

tACS mechanism. (A) Electrode montage shown on a 3D head model (electrode size–35 cm2, current = 2 mA, C3-FP2 montage) targeting the motor cortex using the open source software package SimNIBS (Thielscher et al., 2015). (B) Resulting electric field distribution. Left. Outward (blue) and inward currents induced by the montage C3-FP2. Middle and right. Distribution of the electric field magnitude. The results were analyzed using Gmsh (Geuzaine, Christophe and Remacle, 2009). (C) Oscillations of the electric field are characterized by the period and the amplitude. (D) The membrane voltage of a pyramidal cell changes with the electric field over time. Shown is the normalized membrane polarization for three different time points of the oscillation cycle.

Neuron models

We used multi-compartmental conductance-based neuron models with realistic morphologies from all cortical layers generated in the NEURON environment (Carnevale and M.L. Hines, 1997). NEURON environment is a simulation software well suited for modeling neurons with complex anatomical as well as biophysical properties. It is based on cable theory which discretizes the neuron morphology into small compartments on which the neural dynamics are computed. In these multi-compartmental models, the membrane voltage changes can be investigated using the cable equation along branches as a function of time and space - to calculate the potential, the branches are divided in small compartments of length dx, and the following equation can be numerically solved: Where t is the time, V (x, t) is the membrane voltage, E|| (x, t) is the induced electric field along the neuron branch, c is the membrane capacitance, ri is the membrane resistance per unit length and i corresponds to the ionic currents.

The original data set can be downloaded from ModelDB (Aberra et al., 2018, accession number 066023; McDougal et al., 2017) and includes adapted rat models from the library of Blue Brain models exhibiting realistic biophysical and geometric properties. 3D reconstructed axonal and dendritic ramifications were thus included and their firing behaviors were experimentally validated (Markram et al., 2015). This set is composed of five neurons for each neocortical layer (25 in total): Layer 1 neurogliaform cell (L1 NGC), Layer 2/3 pyramidal cell (L2/3 PC), Layer 4 large basket cell (L4 LBC), Layer 5 pyramidal cell (L5 PC) and Layer 6 pyramidal cell (L6 PC) (see Supplementary Fig. 1 for the morphologies). Different ion channels and myelination are modelled as outlined in detail in (Aberra et al., 2018; Markram et al., 2015): 13 different ionic channels are present in the soma, apical dendrites, basal dendrites and axon initial segment and cover various ion types: transient and persistent potassium and calcium, SK calcium activated potassium, A-type potassium (Kv3.1), D-type potassium (Kv1), M-current (Kv7), stochastic potassium, H-current, high-voltage-activated calcium, and low-voltage-activated calcium. Originally, this data set was used to investigate neural responses to intracortical microstimulation (ICMS) and uniform DC electric field.

Synaptic activity modeling

The original data set are quiescent neurons exhibiting no spontaneous firing activity. In order to investigate how spiking activity is modulated by the tACS intensity, we modified these cells to incorporate intrinsic synaptic activity. To this end, an excitatory synapse was added to each of the 25 cells at a random location on the apical dendrite for pyramidal cells or the basal dendrite for inhibitory cells. Each model corresponds to a single independent neuron and therefore is not connected to other cells. To simulate the intrinsic activity, a stochastic Poisson process was added. We modeled the synaptic activity within a single postsynaptic compartment on the cell body by using a two-exponential function as follows (Destexhe et al., 1998; Rothman, Jason S and Silver, 2014). The synaptic current Isyn(t) that results from the chemical synapse is given by the following expressions:

The parameter E corresponds to the synaptic reversal potential, g(t) is the time-dependent synaptic conductivity, is the maximum synaptic conductance, V (t) is the membrane potential and τ1 and τ2 are respectively the rise time and the decay time of the modeled synapses. The synapse is activated at t = 0, i.e., when a synaptic input is delivered. The parameters values are E to 0 mV, τ1 = 2ms and τ2 = 10 ms (Cavarretta et al., 2014). Spiking activity was generated through a synapse with a Poisson distribution with an identical seed. The Poisson distribution models an artificial pre-synaptic input which triggers a post-synaptic current through time as bi-exponential curve (eq. 2) in our realistic neuron model. When the somatic membrane voltage reaches a threshold value of 0 mV, we considered it as a spike. Since neurons exhibit different morphologies, we had to tune the synaptic parameters – including the synaptic weight and location – for each individual neuron to have a baseline firing rate of 5 Hz to run a fair comparison. We varied the distance between the soma and the synaptic input (median distance from the soma for L1 NGC: 17.62 μm, L2/3 PC: 27.97 μm, L4 LBC: 35.08 μm L5 PC: 46.62 μm, L6 PC: 13.69 μm).

The data set used in this paper is available at GitHub: https://github.com/OpitzLab/TACS_effects_single_neurons

Modeling of tACS

We used a 10 Hz sine-wave tACS stimulation waveform and simulated a range of electric field strengths between 0.1 mV/mm and 3 mV/mm in 0-to-peak amplitude. The amplitude was increased by 0.1 mV/mm steps between 0.1 mV/mm and 2 mV/mm, and by 0.2 mV/mm steps between 2 mV/mm and 3 mV/mm. We chose this range to cover intensities typically induced in humans as well as to explore a higher dosing range which might be achieved in future tACS applications. The electric field is spatially uniform and aligned with the y-axis which corresponds to the somatodendritic axis of the pyramidal cells (Fig. 2A). The total duration of the stimulation was 8 min, divided in two equal periods – a 4 min tACS-free period and a 4 min period with tACS (Fig. 2B). We used the extracellular mechanism of NEURON to simulate the interaction between the electric field and the cell activity (Anastassiou et al., 2010; Nagarajan et al., 1993). This mechanism is used to compute extracellular potentials in the presence of extracellular electrodes or of an extracellular stimulation. The extracellular potential next to the membrane is calculated for each neuron compartment based on the external electric field induced by tACS using the quasipotential , where is the electric field vector and is the differential vector along the neuron branch (Wang et al., 2018). The external electric field was assumed to be uniformly distributed at the small scale of the neuron. The membrane voltage at each compartment is approximated with the voltage at its geometrical center. If the electric field is uniform and parallel to y-axis, the relation between the extracellular membrane voltage and the amplitude of the electric field is linear and is given by is the extracellular potential of the node j, E is the electric field strength and d is the distance between the cathode and the node j (Wang et al., 2018).

Fig. 2.

tACS simulation setup. (A) Examples of the complex realistic neocortical neuron morphologies used in the simulations. Top left. L5 pyramidal cell. Top right. L6 pyramidal cell. Bottom. L4 large basket cell. The orange arrow indicates the direction of the applied electric field along the y-axis. (B) A 10 Hz oscillating electric field was applied to each neocortical cell for 4 min. For comparison another 4 min without tACS was simulated. (C) Spiking activity was generated through a randomly located synapse with a Poisson distribution. Top. Baseline neural spiking activity. Middle. Spiking activity for an electric field strength of 0.5 mV/mm. Bottom. Spiking activity for an amplitude of 1 mV/mm.

Calculating phase-locking of spiking activity to the applied electric field

As a first step, the firing rate was computed as the number of action potentials per second in the time window of 4 min (separately tACS and tACS-free periods). Multiple methods have been developed to measure neural oscillatory synchronization (Lowet et al., 2016; Picinbono, 1997). To quantify the neural entrainment, we computed the phase-locking value (PLV) which estimates spike timing synchronization relative to an oscillation, in our case the tACS waveform. This metric is commonly used allowing a direct comparison with other studies. Although, at low number of observations, PLV can have biased values, it has been shown that with large sample size, PLV is an unbiased metric (Vinck et al., 2010). However, in our simulations, all cells have a high number of spikes (>1200 spikes for each simulation). The PLV is calculated as follows (Lowet et al., 2016): Where N is the number of action potentials and θ is the phase of the tACS stimulation at which the kth action potential occurs. A PLV of 0 means that no synchronization is existent – e.g. in case the spike times are uniformly distributed across the tACS cycle – while a value of 1 means perfect synchronization.

Another tool to quantify neural entrainment is the occurrence of spiking relative to the stimulation phase. The phase histogram counts the phase of tACS at the time when spiking occurs. We computed the tACS phase angle during spike occurrences for each neuron and plotted them for different electric field strengths. In the absence of an electric field, neurons do not fire at a preferred phase with respect to a hypothetical 10 Hz oscillation and should thus follow a uniform phase distribution. If tACS synchronizes the firing activity, one should observe a preferred direction in the phase histogram.

Analysis of the neuron morphology

To identify morphological factors which can explain some of the observed differences in spiking activity due to tACS, we performed a quantitative analysis of the neuron morphology. For this we characterized the neuron morphologies through the effective length (units: μm) along the direction of the electric field (y-axis). The effective length L corresponds to the distance between the two furthest points along the y-axis. A high effective length means that a neuron has a morphology specifically oriented in that axis.

Analysis of the neuron biophysics

Pyramidal cells and interneurons differ in many biophysical properties. To control for the effects of biophysics on neural entrainment, we have performed two control simulations for a baseline firing rate of 5 Hz. The first control simulation is to set all the cells dynamics to the same L5 excitatory dynamics. The second control simulation is to set all the cells dynamics to the same L4 inhibitory dynamics (same ionic channel types and distributions).

Data analysis and simulation

Data analysis was performed using custom MATLAB 2019b (Math-works, Natick, MA) scripts. Neuron modeling and tACS simulation were implemented in the NEURON environment v7.4 (Carnevale and M.L. Hines, 1997). Analysis was performed on a desktop workstation (Intel i7 9700k, 32 GB RAM, Nvidia RTX 2080). The NEURON simulation workflow was executed on the supercomputing cluster of the University of Minnesota. In total, for each of the 25 neuron morphologies we performed 130 simulations representing a total of 3250 simulations. Total computation time was approximately 256,000 CPU-hours.

Results

Cortical cell membrane polarization

In order to determine the membrane polarization with respect to the applied electric field strength, we used the same cells but without any synaptic input (silent neuron). For each cell, we recorded the somatic polarization for a given electric field strength. It has been shown in previous experimental studies that the somatic membrane potential increases linearly with the applied electric field strength (Bikson et al., 2004; Radman et al., 2009). The somatic voltage linearly responses to the subthreshold electric field intensities investigated (0.1 – 3 mV/mm). The slope of each response – called polarization length (λ) – determines the subthreshold somatic polarization per unit electric field applied. The unit is in millimeter . The higher λ, the easier it is to depolarize the neuron. It has been shown that the polarization length can differ strongly across different cell types (Radman et al., 2009). The response of cortical neurons to tACS depends on the specific cell type. We thus first determined the amount of membrane polarization with respect to the applied electric field strength for each modeled neuron. We found a linear relationship (mean R2 = 0.9818 over the set of all neurons) between the soma polarization and the applied electric field strength (see Supplementary Fig. 2A).

Non-pyramidal cells had the lowest polarization lengths (see Supplementary Fig. 2B): L1 NGC (range 0.0076 mm – 0.0187 mm), L4 LBC (range 0.0336 mm – 0.0375 mm). L5 PC show the highest values (range 0.0992 mm – 0.1323 mm) followed by L6 PC (range 0.0944 mm – 0.1204 mm) and L2/3 PC (range 0.0221 mm – 0.0726 mm). We found a significant difference in the polarization length for interneurons across layers and L5 pyramidal cells (t(13) = −13.1963, p = 6.6320*10−09), for interneurons across layers and L6 pyramidal cells (t(13) = −11.5400, p = 3.3350*10−08), and between interneurons across layers and L2/3 pyramidal cells (t(13) = −3.1292, p = 0.008). For a similar analysis see (Aberra et al., 2018). The polarization length mainly depends on the neuron morphology, but the cell dynamics also play a role to some degree (see Supplementary Figs. 3 & 4).

Effect of tACS on the firing rate

To investigate how the rhythmic membrane de- and hyperpolarization affects spiking activity, we first determined the effect of tACS on the firing rate, by comparing the firing rate during tACS to the firing rate without tACS. This was done by calculating the ratio of firing rate during tACS relative to the no-tACS baseline. A value of 1 indicates no change and values smaller or larger than 1, respectively indicate a decrease or increase in the firing rate during tACS. We investigated how the firing rate is affected by increasing electric field strength for various cell types. All the synaptic inputs were tuned for cells have a similar spontaneous firing rate around 5 spikes/second: 5.5150 ± 0.2634 spikes/second (mean ± standard deviation) for L1 cells, 5.3583 ± 0.2850 spikes/second for L2/3 cells, 5.3342 ± 0.1658 spikes/second for L4 cells, 5.2967 ± 0.1739 spikes/second for L5 cells and 5.3625 ± 0.1661 spikes/second for L6 cells (see Supplementary Table 1 for exact value for each cell). We found that in a low electric field strength regime (< 1 mV/mm), the firing rate changed less than 1% irrespective of the cell type (Fig. 3A; see also Supplementary Fig. 5 + Table 1). For the higher electric field intensities, L2/3 and L5 pyramidal cells showed increased firing rate of up to 5%. Overall, we did not find a major change (< 5%) in the firing rate regardless of the cell type, especially for the electric field intensities achievable in Humans (< 1 mV/mm).

Fig. 3.

tACS effects on firing rate and phase locking-value (PLV). (A) Normalized firing rate (change compared to baseline) as a function of electric field strength. Firing rates of different neurons (different colors indicate different cells) are not affected by weak electric fields. At higher intensities an increase below 5% firing rate can occur. (B) Synchronization of neural firing as measured through the PLV increases for all cells with increasing field strength. Pyramidal cells show a larger increase in PLV than other cell types. This increase is present even in the low electric field strength regime. Only the mean of each cell type is displayed.

Effect of tACS on neural entrainment

To quantify the effect tACS has on neural entrainment, we calculated the PLV for all cells as a function of electric field strength. We found that the neural synchronization increases with electric field strength (Fig. 3B, see also Supplementary Fig. 6 for all cells). On average, pyramidal neurons from Layer 5 exhibit the strongest entrainment followed by neurons from L4 (mean slope of 0.0616 and 0.0408, respectively). Pyramidal neurons from Layer 2/3 and Layer 6 exhibit a similar entrainment while interneurons from Layer 1 show no entrainment regardless the electric field intensity (mean slope of 0.0179, 0.0154 and 0.0037, respectively). We explored two additional baseline firing rates: 10 Hz (+/− 1 Hz) which corresponds to the stimulation frequency and 15 Hz (+/− 1 Hz) which is higher than the stimulation frequency. We found consistent results for these two additional firing rates (Supplementary Figs. 11–12).

To further investigate the preferential timing of neural firing with respect to the tACS oscillation, we computed phase histograms (Fig. 4). We found that with increasing electric field strength, phase histograms deviated more from uniformity. Phase distributions differed from circular uniformity (see Rayleigh test statistics for all neurons in Supplementary Table 2 and Supplementary Figs. 13–17) for electric field strengths as low as 0.5 mV/mm. For an electric field strength of 1 mV/mm, only two out of ten non-pyramidal cells exhibit a significant phase preference, and two out of fifteen pyramidal cells exhibit one. Additionally, the preferred phase angle differed between the investigated neurons. For example, the preferred phase angle for the Layer 5 pyramidal neurons was 300–350°, which corresponds to the moment just before the tACS peak. This is in line with experimental literature (Johnson et al., 2020; Krause et al., 2019) showing that preferred phase does not necessarily have to occur at 0-degree (maximum electric field strength) but can also occur at other phase angles.

Fig. 4.

Polar histograms of tACS phase during spike onset for L5 pyramidal cells. Top. Each plot shows the entrainment for different electric field intensities for each cell. Each color corresponds to an electric field intensity - the darker the color is, the higher the intensity is (intensities displayed: 0, 0.5, 1, 1.5, 2 and 3 mV/mm). Radial axis corresponds to the PLV (The farther away from the center, the higher PLV). 0-degrees correspond to the peak of the oscillation and 90-degrees correspond to the falling edge of the oscillation. Bottom left. Morphology of a L5 pyramidal cell (scale bar = 100 μm). Bottom right. Polar histograms across different electric field intensities for the displayed L5 pyramidal cell. Radial axis corresponds to the number of spikes fired at that specific phase. The thick red line corresponds to the mean direction. The stronger the electric field gets, the more the cell converges to a certain preferred phase (324° for this cell).

Influence of neuron morphology on tACS effects

A neuron subjected to an electric field will exhibit a strong polarization at its distal ends. Thus, the more elongated the somatodendritic axis is, the more polarized the cell membrane will be, resulting in stronger tACS effects on spike timing. In order to quantify the morphology for each neuron, we computed the effective length L which is the distance from the two distal ends of a neuron in the y-axis direction corresponding to the direction of the electric field. We first investigated the relation between the polarization length and the effective length L and we found a high correlation (r(23) = 0.8520, p = 6.58e-08, see Figure 5B). This means that a cell that is elongated in the direction of the electric field (i.e. a large effective length) will exhibit a large net polarization and thus larger effects on its spiking activity (Radman et al., 2009). Therefore, neurons with a higher effective length L are more likely to be entrained by an electric field. We then investigated the relation between the PLV and the effective length L for all 25 neurons at electric field strength of 1 mV/mm. We found a correlation coefficient of r(23) = 0.2452, p = 0.237 with PLV (see Figure 5C). The effective length L is thus easy-to-derive morphological features that can explain the sensitivity of large pyramidal cells to tACS.

Fig. 5.

Response to tACS depends on the neuron morphology. (A) Neocortical cells have a large diversity of morphologies and pyramidal cells exhibit a pronounced somatodendritic axis resulting in a high effective length Ly in the direction of the electric field. Displayed are Layer 1 and Layer 2/3 cells. Each horizontal black line represents 100 μm. (B) Relation between the polarization length and the effective length Ly. The higher the effective length is, the easier the cell is polarized. The dashed line represents a linear regression across all the points. The polarization length is highly correlated with the effective length (r(23) = 0.8520, p = 6.58e-08). (C) Relation between the effective length and the PLV. The dashed line represents a linear regression across all the points computed for an electric field of 1 mV/mm. We found a correlation coefficient of r(23) = 0.2452, p = 0.237 with PLV.

Influence of neuron biophysics on tACS effects

To study the effect of cell biophysics, we performed two additional simulations to control for the effects of biophysics on neural entrainment. The first control simulation is to set all the cells dynamics to the same L5 excitatory dynamics. The second control simulation is to set all the cells dynamics to the same L4 inhibitory dynamics (same ionic channel types and distributions).

Fig. 6 shows that in both excitatory and inhibitory dynamics, the higher the effective length is, the more the PLV increases meaning that neurons exhibiting a high effective length are more likely to be more entrained. However, neurons with inhibitory dynamics exhibit a higher slope meaning that inhibitory dynamics makes them more sensitive to electric field and entrainment. Moreover, the linear regression exhibits a high coefficient of correlation. These results highlight the idea that both morphology and dynamics play important role in entrainment.

Fig. 6.

Neuronal dynamics affect suprathreshold activity. (A) Excitatory dynamics do not highly affect neural entrainment compared to default dynamics. The dashed line represents a linear regression across all the points (r(23) = 0.2088, p = 0.316). (B) By contrast, inhibitory dynamics highly affect neural entrainment of cells – mainly of L5 pyramidal cells. This can be explained by the inhibitory dynamics on the subthreshold activity which increase the polarization length of pyramidal cells (see Supplementary Fig. 3.B). The dashed line represents a linear regression across all the points (r(23) = 0.6257, p = 0.000822). Shown is the PLV for an electric field intensity of 1 mV/mm.

All the cells with excitatory dynamics exhibit a constant firing rate – or exhibit a small increase up to 5% (Supplementary Fig. 7). L5 pyramidal cells exhibit the strongest entrainment followed by pyramidal cells from other Layers. Cells from Layer 1 and Layer 4 show no sign of entrainment as it can be observed on Supplementary Fig. 8.

In the second control analysis (inhibitory neuron dynamics), pyramidal cells exhibit a higher increase in the firing rate – up to 30% for L2/3 and L5 cells. L6 cells exhibit an increase in the firing rate up to 13% (see Supplementary Fig. 9). We also observed that pyramidal cells exhibit a significantly higher entrainment (Supplementary Fig. 10) while inhibitory neurons from Layer 1 show no entrainment.

Discussion

In this study we developed computational models of morphologically realistic neocortical neurons with ongoing spiking activity to study the effects of tACS electric fields on neural activity. We demonstrate 1) a linear relationship between membrane cell polarization and the applied subthreshold electric field strength, 2) varying polarization lengths across different cell types resulting from their different morphologies, 3) no major changes in the firing rate of investigated neurons for electric field strengths applicable to human experiments, 4) spike timing changes and neural entrainment at low electric field strengths, 5) neural entrainment is both affected by the neuronal morphology and the cell biophysics.

Our results are in line with previous experimental and computational work highlighting the importance of morphologic features of neurons on somatic membrane polarizations to electric field stimulation (Berzhanskaya et al., 2013; Chan, Hounsgaard, 1988; Radman et al., 2009; Yi et al., 2017). Pyramidal neurons are characterized through asymmetrically placed somas with respect to their dendritic ramifications. In contrast, most interneurons have a dendritic tree uniformly distributed in space (Amitai and Connors, 1995). These differences in neuronal morphologies result in differences in polarization across cell types (Bikson et al., 2004; Chan, Hounsgaard, 1988; Francis et al., 2003; Jefferys, 1981) which can be quantified through the polarization length. Interestingly, polarization lengths computed in our study overlap very well with experimentally determined values corroborating the utility of our modeling approach (Radman et al., 2009).

We further investigated the effects of tACS on firing rate and spike timing. We did not observe a considerable change in firing rates for the investigated neurons in an electric field regime applicable to human tACS studies. This is in line with recent in-vivo recordings in non-human primates (Johnson et al., 2020; Krause et al., 2019). Only for electric field strengths > 2 mV/mm, did noticeable increases in the firing rate occur. However, we observed a significant modulation of spike timing across all cell types. Entrainment of spiking was present even for low electric field strengths and increasing for higher amplitudes. Our findings support experimental results in non-human primates showing that electric fields strengths of 0.5 mV/mm are sufficient to cause entrainment of spiking activity (Johnson et al., 2020). These effects were most prominent in L5 pyramidal cells which showed stronger phase-locking of spikes compared to interneurons. The effectiveness of tACS to entrain spiking activity was strongly determined by morphological factors such as the effective length of the neuron along the electric field vector. Additionally, we found that neuronal biophysics such as ionic channels distribution play an important role and neural entrainment cannot be entirely explained with morphologies. These findings highlight the importance of taking into account neuron morphology with respect to the applied electric field orientation to optimize tACS physiological outcomes. While our model was able to replicate known experimental findings and highlight important mechanisms of tACS, it has some limitations. First, we only chose one excitatory synapse in proximity to the soma (within a 50 μm range). The synapse location at one proximal dendrite cannot fully take into account nonlinearities and effects arising from dendritic integration (Spruston, 2008). tACS effects on the single neuron level can be amplified through network effects arising from synchronized synaptic input. Depending on the size of neurons populations and on the synaptic connections between them, tACS effects can be magnified by dynamic network activity and these effects need to be investigated in future studies. Future efforts including multiple sources of excitatory and inhibitory synaptic input – NMDA, GABAA or GABAB for example - in our neuron models could overcome some of these limitations and enable predicting tACS effects occurring in larger neural networks. Even though our set of neurons studied here includes various morphologies across all neocortical layers, our findings do not cover all cortical neurons. Moreover, cells from the same type can exhibit large variability in their response to tACS highlighting the need for further large-scale studies. The complex relationship between electric field strength and PLV could be investigated in more detail in future studies. Further work could also explore different frequencies, stimulation waveforms or orientations. Animal studies have shown that the firing rate of neurons increases with a steeper temporal change (Fröhlich and McCormick, 2010) and several human studies have investigated the effects of stimulation waveforms on brain oscillations (Dowsett et al., 2020; Dowsett and Herrmann, 2016) showing that different waveforms differentially modulate brain oscillations.

In conclusion, our findings advance our mechanistic understanding of tACS effects on single neurons. Our results highlight the importance of neuronal morphology on the susceptibility towards external electric fields. Morphological differences across neurons result in differences in membrane polarization and consequently effects of neuronal entrainment to tACS. Changes in spike timing can already occur at low electric field strengths achievable in human studies, most prominent in pyramidal neurons. In future work, our modeling framework will allow the efficient exploration of tACS parameters thus enabling the optimization of stimulation protocols to maximize their physiological effectiveness.