## 2 Materials and Methods
### Data Analysis
#### 2.3.4 Methodology: Part B  
We divided our exploratory DL analysis into four phases, each designed with a particular objective in mind:  
1. To determine whether EEG during TEAS differed from EEG at baseline, and whether such differences were dependent on stimulation frequency;
2. To investigate how EEG differed before, during and following TEAS, and whether such differences were dependent on stimulation frequency;
3. To determine whether differences among the EEG bands vary with both time (baseline, stimulation or post-stimulation) and stimulation frequency;
4. The objective here was the same as in Phase 3, but with the addition of comparing how the MLP-NN and CNN-LSTM models performed. Methodology: Part B, Phase 1  
For our analyses, here, we only processed data recorded at 500 Hz, although this did include some quite noisy recordings, resulting in 1806 or 1807 files for each filtered EEG band (delta to gamma). In this phase, the data from all five EEG bands were considered together rather than separately. The data were then standardised using the StandardScaler() command in the Python scikit-learn library [45]. Our objective was to determine whether EEG during TEAS differed from EEG at baseline, and whether such differences were dependent on stimulation frequency. The initial analysis was conducted using a Sequential Application Programming Interface (API) in the 'Keras' DL framework [46], and a hybrid CNN-LSTM model similar to that illustrated in Figure S5 in the online Supplementary Materials, as shown in Table 4.  
The outputs were evaluated using standard confusion matrix metrics, of which sensitivity and specificity are probably the best known. Here, we focused on accuracy, _kappa_ and the area under the'receiver operating characteristic' curve' (ROC curve), known as the AUC or ROC-AUC [47]. If the mean overall accuracy was <0.33, the results were considered non-significant. Sixteen models were created, to examine changes over time (Slots 2, 3, 4 and 5) relative to the baseline (Slot 1), for each of the four stimulation frequencies. Methodology: Part B, Phase 2  
The feature extraction methodology utilised here is based on a combination of the methods used in our previously published papers [48; 49; 41], with data standardised using the StandardScaler() command in the Python scikit-learn library. Standardising the input data can help to reduce the effect of outliers or extreme values and improve the performance of the machine learning model. Some algorithms are particularly sensitive to the scale of the input data (see Supplementary Materials), and thus, may perform better when the data are standardised. Our objective was to investigate how EEG differed before, during and following TEAS, and whether such differences were dependent on stimulation frequency. As in Phase 1, the analysis was conducted using a Sequential API in Keras and a hybrid CNN-LSTM model, but this time, omitting the Dropout layers and adding a further LSTM layer (Table 5). Omitting the Dropout layers and adding a further LSTM layer can potentially improve a model's performance on the training data (but can also lead to overfitting if the model is not properly regularised). It is generally considered good practice to try a variety of different model architectures and regularisation techniques to find the best balance between model complexity and generalisation to new data. Twenty models were created, to examine changes over time (baseline (Slot 1), stimulation (Slots 2-5) and Post-stimulation (Slots 6-8)) in each of the five EEG bands (delta, theta, alpha, beta and gamma) for the four stimulation frequencies. Methodology: Part B, Phase 3  
Here, we did not use Keras and the DL hybrid CNN-LSTM algorithm, but a (relatively) shallow learning method based on Multilayer Perceptron (MLP), a scaled conjugate back-propagation MLP Neural Network (NN) with 10 hidden layers, available as standard in the MATLAB Deep Learning toolbox [50].