Electroencephalography data were recorded with a MR-compatible EEG cap with 31 sintered Ag/AgCl electrodes mounted on an elastic cap (10–20 system) (Easycap, Falk Minow, Munich, Germany). An additional electrode was attached on the subject’s back to record the electrocardiogram (ECG). The electrode impedance was kept below 20 kΩ. The signal was recorded using BrainVision Recorder software (Brain Products GmbH, Gilching, Germany) with an online bandpass filter (0.016–250 Hz). EEG data were sampled at a rate of 5,000 Hz and amplified using a BrainAmp MR plus amplifier (Brain Products, Munich, Germany).

All images were acquired using a 3 Tesla Magnetom Prisma scanner (Siemens Healthcare, Erlangen, Germany) with a 64-channel phased-array head coil. In the first 7–8 min, adjustment scans and a high-resolution T1-weighted anatomical scan (three-dimensional (3D) turbo fast low angle shot; repetition time (TR): 2,250 ms, inversion time (TI): 900 ms, echo time (TE): 3.3 ms, flip angle 9°, isotropic resolution of 1 mm × 1 mm × 1 mm) were recorded. During the further procedure, functional data were acquired using gradient-echo, echo-planar imaging with T2*-weighting (TR: 2 s, TE: 30 ms, flip angle 50°, 35 slices of 3-mm thickness, gap between slices of 20%, in-plane resolution of 3 mm × 3 mm). A total of 210 whole-brain volumes were acquired in each functional run. Before every functional scan, an additional scan with the same parameters but opposite phase encoding direction was recorded to correct for distortion in the functional data.

To improve the individualized volume conduction model for electrical field simulation (see section “electrical field simulation”), we recorded additional anatomical scans in a separate session, including a T2-weighted anatomical scan (3D turbo spin echo; TR: 3500 ms, TE: 282 ms, variable flip angle, with integrated parallel acquisition technique: factor 2, isotropic resolution of 1 mm × 1 mm × 1 mm, no fat saturation), and diffusion-weighted MRI performed using a spin-echo echo planar imaging technique at 1.7 mm isotropic resolution (repetition time = 10,000 ms, echo time = 88 ms, parallel acquisition factor 2; acquisition matrix: 128 × 128, 75 slices), acquiring 64 image volumes with diffusion weighting (along 64 diffusion directions, b = 1,000 s/mm²) and one reference image without diffusion weighting.

Due to technical problems, T2 images are lacking in seven subject datasets. The T2 contrast was used to improve brain matter segmentation [skin, skull, cerebral spinal fluid (CSF), white matter, gray matter] for calculating the volume conduction model needed for tACS electrical field simulation.

After each experiment, the subject completed a questionnaire designed to assess their sleepiness during the experiment and get information of possible sensation effects of the stimulation. Sleepiness was rated for each run with a handwritten mark on a 10-cm long line representing a continuum between very awake and very sleepy [Visual analogue scale (VAS); Funke and Reips, 2012].

The data were analyzed using MATLAB R2015b (The MathWorks Inc, Natick, MA, United States). For analyzing the EEG data, we used the Fieldtrip toolbox (Oostenveld et al., 2010) and EEGLAB (Delorme and Makeig, 2004) in combination with the FMRIB plug-in provided by the University of Oxford Centre for Functional MRI of the Brain (FMRIB) (Niazy et al., 2005). FMRI data were analyzed using FSL (Smith et al., 2004; Jenkinson et al., 2012) for top-up correction, and the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012) for further analysis steps. For electrical field stimulation we used SimNibs (Windhoff et al., 2013; Thielscher et al., 2015).

Gradient and cardioballistic artifacts were attenuated using optimal basis sets (OBS) as implemented in FMRIB Plug-in for EEGLAB (Delorme and Makeig, 2004; Niazy et al., 2005). Further analysis steps were performed using the Fieldtrip toolbox (Oostenveld et al., 2010). Data were downsampled to 1,000 Hz and low-pass filtered at 160 Hz. An independent component analysis (ICA) (runica method, 20 components) was calculated on a band-pass filtered version (1–20 Hz) of concatenated pre and post EEG data (separately per session). Components representing typical spatial topographies of blinks, saccades, residuals of cardioballistic artifacts or residuals of gradient-related artifacts were visually identified (on average three to five components per subject). The mixing matrix of the ICA decomposition model was then applied to the 160 Hz low-pass filtered data. Subsequently, the artifactual components were removed before the data were back-projected to channel space. These data were then re-referenced to an average of all channels and visually inspected, whereas artifact-contaminated parts were marked for later removal.

The stimulation frequency was determined from the 3-min eyes-closed recording period before the main experiment started (see section “study protocol”). Using BrainVision Analyzer 2.0 (Brain Products GmbH, Gilching, Germany), the EEG data were fast Fourier transformed with a frequency resolution of 0.2 Hz, and the frequency with the most dominant peak between 8 and 12 Hz from the occipital electrodes (O1, Oz, O2, POz) was visually identified. This frequency was used as the stimulation frequency. Since the IAF peak location can vary slightly across conditions, e.g., for eyes-closed and watching the video and also from day to day, we calculated an averaged IAF from all EEG data recorded during the video for further data analysis. Therefore, artifact-free EEG segments of pre and post in both sessions were split into 10 s epochs. Next, a Fourier transformation was calculated [0.1 Hz resolution, 0.5 Hz smoothing, multi taper (dpss)] for the Oz and POz electrodes and subsequently averaged over both channels. The frequency with the highest peak between 8 and 12 Hz was taken as IAF.

To evaluate changes of alpha activity between pre and post and between sessions, we determined the alpha amplitude as follows: a Hilbert transformation of preprocessed EEG data (bandpass filtered ± 2 Hz around IAF, 4th order Butterworth filter with hamming window) was calculated for Oz and POz channels separately. The absolute values of the complex Hilbert transformation, which represented the envelope of the filtered signal, were averaged in artifact-free segments over time and over both channels (Nelli et al., 2017). For each subject we calculated the difference in the alpha amplitude change between the two sessions (ΔΔ ALPHA_AMP) as Δ alpha amplitude (post - pre) calculated for the tACS session minus Δ alpha amplitude (post - pre) calculated for the sham session.

In a first step, fMRI data were top-up (Andersson et al., 2003) corrected using FSL (Smith et al., 2004; Jenkinson et al., 2012) and further analyzed using the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012). After correction for motion, slice timing correction and outlier detection, functional data were spatially normalized into MNI space. The data were then downsampled to 2 × 2 × 2 mm and smoothed with a Gaussian kernel (6 mm FWHM). In a subsequent denoising step, six motion regressors, the mean signals of white matter and CSF, and a various number of noise components (one for each detected outlier) were regressed out of the data. Finally, functional data were band-pass filtered between 0.008 and 0.24 Hz.

To determine the voxelwise cortical regions that are related to the EEG alpha band activity recorded on the scalp, we performed an EEG-guided fMRI analysis. A general linear model (GLM) was calculated for every voxel of the preprocessed fMRI data, using a regressor derived from the EEG alpha amplitude. For each run of each subject, an alpha-related regressor was calculated using the signals from the Oz and POz electrodes. Separately, the continuous (7 min) EEG signals were Hilbert-transformed (bandpass filtered ± 2 Hz centered around IAF, 4th order Butterworth filter with hamming window) and their absolute values were convolved with the hemodynamic response function (HRF). Subsequently, the time courses were downsampled according to the fMRI acquisition rate (0.5 Hz) and, finally, averaged between electrodes.

Results of the EEG-guided fMRI analysis were then used to create binary masks to determine the electrical field strength at alpha related occipital voxels (see section “electrical field simulation”). For each subject, we created a binary mask that included the 1,000 occipital gray matter voxels showing the strongest (negative) alpha-BOLD correlations (MASK_alphaBOLD). We focused on only negative correlations because previous studies have shown that in occipital areas, BOLD signal and occipital alpha amplitude are inversely related (Goldman et al., 2002; Laufs et al., 2003, 2006). Then, the first-level results of all runs (in total four: pre and post of tACS sessions plus pre and post of sham sessions) were averaged separately per subject. The voxels in the occipital-parietal areas with the most negative t-values were determined individually for every subject and served as a binary mask for further analysis steps.

To summarize DMN connectivity strength, mPFC-PCC functional connectivity was calculated for each run of each subject using a ROI-based correlation approach. ROIs were defined as 10 mm spheres centered around MNI coordinates adopted from CONN’s ICA analyses of the HCP dataset (497 subjects) [mPFC (1,55,–3) and the PCC (1,–61,38)]. For each ROI, the BOLD time course was averaged over all voxels. Connectivity was calculated as the Spearman correlation between the two ROI time courses and subsequently Fisher-z-transformed (Chang et al., 2013; Clancy et al., 2022). Changes in network connectivity within a session were calculated as the difference between the Fisher-z-transformed correlation factor from that session’s post and pre (yielding Δ CONNECTIVITY_tACS and Δ CONNECTIVITY_sham). The difference in DMN connectivity strength between conditions (ΔΔ CONNECTIVITY) was calculated as Δ CONNECTIVITY_tACS minus Δ CONNECTIVITY_sham.

For each run, subjects reported their experienced degree of sleepiness by marking the appropriate point on a continuous line with extremes representing very sleepy and very awake. The marked position was quantified in relation to the full length of the line giving a sleepiness value between 0 and 1 for every run. In the first step we determined the change of sleepiness during each session by subtracting the pre sleepiness value from the post sleepiness value. In the next step these results were subtracted from each other to determine the difference in sleepiness change between the tACS and the sham session (ΔΔ SLEEPINESS).

SimNIBS (Windhoff et al., 2013; Thielscher et al., 2015) was used to model the electrical field induced by tACS for each subject. First, for each subject, a head model was calculated using T1- and T2-weighted (if available) images. The head models contained five compartments (scalp, skull, gray matter, white matter and CSF) and used a mesh resolution of 0.5 nodes per mm². In the second step, diffusion tensor imaging (DTI) data were integrated into the head model. Simulations were performed using the following experimental parameters in SimNIBS: 2 cm round electrodes in PO7 and PO8 positions (determined via a template cap provided by SIMNIBS), 1 mA peak-to-peak current and an estimated thickness of 0.2 cm conductivity cream between electrodes and scalp. For the anatomical components, standard conductivities were used: scalp = 0.465 S/m, skull = 0.010 S/m, CSF = 1.654 S/m, gray matter = 0.275 S/m, white matter = 0.126 S/m (Wagner et al., 2004). The results of the electrical field simulations were transformed for each subject into MNI space.

A mean electrical field strength was calculated for each subject individually by averaging all simulated electrical field values covered by a predefined mask. In the next step these averaged values were used as a predictor in a general linear regression model. Since we were aiming to find out if the EFIELD strength predicts the aftereffect better when it is determined for regions related to alpha amplitude (EFIELD_alphaBOLD) as compared to the maximum values that reach the occipital cortex (EFIELD_strong), we used two alternative masks per subject. To determine EFIELD_strong the mask included 1,000 gray matter voxels for which the highest electrical field values were estimated (MASK_strong). Alternatively, to determine EFIELD_alphaBOLD, the mask included 1,000 occipital gray matter voxels with the strongest negative correlation between alpha amplitude and BOLD signal (MASK_alphaBOLD, see section “alpha-BOLD correlation”).

To evaluate differences within and between sessions for alpha amplitude, sleepiness and DMN connectivity, we performed a two-way Analysis of Variance (ANOVA) with the two within-subject factors session (tACS and sham) and run (pre and post) for each of these measures.

To test for a difference between sessions (tACS and sham) while considering sleepiness as a covariate, we calculated an Analysis of Covariance (ANCOVA) for alpha amplitude and DMN connectivity strength. For this analysis we subtracted the pre recording values from the post recording values and calculated the ANCOVA with the factor session (tACS and sham) and included ΔSLEEPINESS as a covariate.

To disentangle the effect of tACS on the alpha amplitude and the mPFC-PCC connectivity strength from those effects that appeared because of changes in sleepiness, multiple linear regression models were calculated. The models included the covariate ΔΔ SLEEPINESS and the two tACS-related predictors, electrical field strength (EFIELD) and frequency mismatch (MISMATCH). Frequency mismatch was calculated for each subject as the absolute difference between the tACS frequency and the mean IAF during the tACS session.

To test the generalizability of our models, we used a cross-validation technique. Leave-one-out cross-validation (LOOCV) is a validation method that is used on small datasets like our cohort of 18 subjects. For the validation, the model is fitted n times (n = number of all observations) on n-1 data points to predict the data point that was left out in each case. The root-mean-square error (RMSE) obtained between predicted and observed data points quantifies the generalizability of the original model to a new dataset.

Figure 2 shows an overview of the analysis steps.

Although all subjects were stimulated with 1 mA on the scalp, the strength of the EFIELD that finally reached the brain varied strongly between subjects due to the individual anatomy and resulting conductivity (Thielscher et al., 2015; Antonenko et al., 2021). This consideration makes the individual EFIELD strength a crucial factor in explaining inter-subject variability in the tACS-related amplitude increase. The result of our model shows that the individual EFIELD strength is significantly positively related to an increase in alpha amplitude from pre to post. This relation was even stronger when sleepiness as an intrinsic alpha amplitude modulating factor was regressed out before EFIELD and alpha amplitude were correlated. This is in line with previous findings, which showed that a stronger EFIELD induces a stronger tACS aftereffect (Johnson et al., 2020).

However, it needs to be established how to utilize the mean EFIELD strength based on the individual simulation maps. For example, we averaged the values of the simulated EFIELD strength in two alternative regions of the occipital cortex: MASK_alphaBOLD and MASK_strong. The accuracy of our model improved greatly when the individual EFIELD strength was determined in brain areas that are related to alpha activity (EFIELD_alphaBOLD), compared with averaging the strongest EFIELD values for each subject (EFIELD_strong). Consistent with these results, Kasten et al. (2019) showed in their model the importance of considering the overlap of regions with the strongest induced EFIELD and regions estimated to be alpha generators. In order to identify alpha generating areas, Kasten performed a source estimation based on high density MEG data, while we determined alpha-related occipital areas via an alpha amplitude-guided fMRI analysis. Although they used different analyses approaches, both studies show the importance of considering the EFIELD strength in brain areas that are related to the alpha amplitude. Thus, for describing the efficiency of tACS on the alpha amplitude, it might be beneficial to not only determine the maximum EFIELD values for each subject, but rather include the values of the EFIELD strength that exist in the target region, in our case alpha-related areas.

To calculate the mean EFIELD from the simulated maps, Kasten et al. (2019) incorporated an area consisting of 10,000 voxels, including gray and white matter voxels, while we decided to incorporate a much smaller area (1,000 voxels encompassing only gray matter), because our approach was planned to study alpha-specific areas. However, recalculating our models using EFIELD strengths averaged over 10,000 gray matter voxels led to similar results (see Supplementary Table S1). Since both area sizes (1,000 and 10,000 voxels) are somewhat arbitrary, the optimal size for the EFIELD strength calculations remains an open issue.

Transcranial alternating current stimulation-related aftereffects on the alpha amplitude might be undetected or misinterpreted in cases where the natural alpha modulation arising from continuous changes in brain state, e.g., due to increased sleepiness, covers the tACS-induced amplitude change. To disentangle alpha amplitude changes due to mental state changes from those induced by tACS, we assessed self-rated sleepiness as an independent variable that represents an alpha-related brain state. Our results imply that including self-rated sleepiness in a linear regression model can reveal tACS-related amplitude effects that would have otherwise remained undetected. Nevertheless, calculating an ANCOVA and testing the difference between conditions (tACS, sham) while considering sleepiness as a covariate did not reveal significant differences between conditions on the alpha amplitude and also not on the mPFC-PCC connectivity strength. These results support the importance of additionally considering the EFIELD for proving the effectiveness of tACS as discussed above. Although EFIELD showed a significant correlation with ΔΔ ALPHA_AMP, even this correlation was much stronger after regressing out the alpha amplitude effects based on changes in sleepiness. Thus, sleepiness could be a crucial factor in observing a tACS effect targeting the alpha amplitude in experiments, which without stimulation presumably have a strong change of the alpha amplitude based on mental brain changes (e.g., exhausting task or dark environment). It is also possible that sleepiness is particularly important in studies in which subjects are in a supine position, e.g., in the MRI scanner.

The absolute mismatch between the stimulation frequency and the IAF (mean IAF during the tACS session) had no significant contribution in explaining the variability of tACS-induced aftereffects in our model, neither for alpha amplitude changes nor for mPFC-PCC connectivity changes. This means that, in our study, the offset between stimulation frequency and IAF had no influence on the tACS aftereffect.

Previously, only a few tACS studies took the possible influence of frequency mismatch on the tACS aftereffect into account, with contradictory results (Vossen et al., 2015; Stecher et al., 2017; Stecher and Herrmann, 2018; Kasten et al., 2019). Stecher, for instance, reported a relationship between frequency mismatch and aftereffect on the alpha amplitude in their study (Stecher and Herrmann, 2018). Kasten showed a significant contribution of frequency mismatch in a model explaining inter-subject variance in tACS-induced alpha amplitude increase (Kasten et al., 2019). They considered negative and positive values for the frequency mismatch (IAF above or below stimulation frequency), and they found a stronger tACS aftereffect with stimulation frequencies slightly above IAF. Vossen et al. (2015), on the contrary, suggested that a stimulation frequency below the IAF had a stronger aftereffect on the alpha amplitude.

From a theoretical perspective, a perfect match between the tACS frequency and IAF should indeed matter: A common hypothesis about the tACS mechanism is that intrinsic oscillations are entrained during tACS application, leading to an increase in alpha amplitude. This entrainment requires less tACS intensity if the frequencies are close to each other (Pikovsky et al., 2001). However, it should be noted that in most of the studies, including ours, the mismatch was calculated between the stimulation frequency and a mean of IAF determined before and after stimulation (EEG recording is difficult during stimulation due to the large artifacts induced by tACS). Thus, we cannot really know the true offset and difference between IAF and stimulation frequency during tACS application. In this context, it is also necessary to mention that the subject’s IAF is generally not stable during any experimental procedure (see also Supplementary Figure S1) because the dominant peak frequency varies depending on the mental state. A perfect match between the IAF and tACS at coinciding time points would require a continuous adjustment of the stimulation frequency in a closed loop system. However, even when an intermittent closed loop stimulation protocol is employed, the results might not be optimal. In a recent study, the parietal alpha rhythm was targeted using a closed loop system. The hypothesis was that the closer match of intrinsic and stimulation frequency should lead to an increased detection of visual luminance changes depending on a stronger effect on the alpha peak power in comparison with a protocol using a fixed stimulation frequency (Stecher et al., 2021). The results showed that only the fixed stimulation protocol led to a persistent increase in post-stimulation alpha power compared with sham. Our results also suggest that since the effectiveness of the stimulation did not depend on a perfect match between intrinsic alpha frequency and tACS frequency, such a technical effort might be unnecessary. It should also be noted that in the previously mentioned studies the frequency mismatch was larger than in the present study, which might explain why this factor did not contribute to the variance in our model. Furthermore, Kasten et al. (2019) also considered interactions between frequency mismatch and the strength of the EFIELD. From a conceptual perspective, it makes sense that the frequency mismatch can only have an effect on the alpha amplitude if the EFIELD is strong enough and vice versa. This interaction was not considered in the present study.

In summary, the results of the present study support the importance of individually calculated EFIELD strength along with self-rated sleepiness as two crucial factors that can explain and predict a large amount of inter-subject variability of the alpha amplitude change after tACS applied over occipito-parietal cortex. We found that determining the EFIELD strength in individual, alpha-related brain areas has an advantage over just averaging the individual maximum current strength values. In future studies, individualized dose-control could probably eliminate the variance in EFIELD intensities at a cortical target site. Supposing that the current delivered to the brain directly determines its behavioral consequences, this method may allow for reducing the known variability of tACS effects.