## 2 Materials and Methods
### Machine Learning Methods  
Two dimensionality reduction algorithms (PCA and t-SNE) and three classification algorithms (SVM, RF, and Catboost implementation of XGBoost) were applied to the THz spectral data. PCA (sklearn.decomposition.PCA) was applied along with t-SNE (sklearn.manifold.TSNE), without standard scaling. The proposed ML pipeline (Figure 4) includes the following steps: pre-processing; ML model construction and validation; and informative feature selection. The ML model construction and validation used the same algorithms, but different optimized parameters. The optimal values for each pair of groups are calculated based on accuracy metrics and given in Table S2, Supplementary Materials. Figure 3: The averaged THz spectra in the time domain (**a**) and the frequency domain (**b**) for all studied groups. During the pre-processing, each THz spectrum was smoothed by the Savitsky-Golay filter [68]. The impact of the filter parameters (filtering window width--window_size, and degree of polynomial smoothing--polyoder) on the classifier performance was estimated by the cross-validation procedure. The parameter names are given according to the library scipy.signal.savgol_filter in Python. Next, the THz data were normalized by the maximum intensity value (vector normalization) for an appropriate spectra comparison. The effect of data centralization (statistical procedure of subtracting mean and dividing by variance) on the dimensionality reduction methods was tested. The unsupervised ML methods (PCA and t-SNE) were used to verify the homogeneity of the CMI groups (Figure 5). Figure 4: Proposed ML pipeline steps. Figure 5: Illustration of application of unsupervised ML methods (t-SNE and PCA) for the separability analysis of the CMI groups. The number of principal components (PC) for PCA was limited to ten, with the highest explained variance for the visual comparison. The optimized t-SNE parameters were perplexity (15) and number of iterations (300). The latter did not have any impact on the separability of the data. The perplexity significantly changes the distribution of points after projection, so the choice of an appropriate value was crucial. Catboost was chosen because of its excellent computational performance. The optimal values of the hyperparameters for the RF and Catboost ML models were estimated by the RandomizedSearchCV algorithm from the sklearn.model.selection library. Using cross-validation on a user-defined parameter search grid, the algorithm determined the prediction data model with the highest performance metrics for a definite ML algorithm. The choice of hyperparameters for optimization usually depends on human expertise. A dense parameter grid will lead to unreasonable computational costs, and even evolutionary optimization algorithms will require a huge amount of processing time. Randomized search is used to speed up the process. This step was performed to avoid instability of the feature selection algorithms on the models with default parameters. The linear kernel SVM performance can be optimized by a regularization parameter (C). A lower C avoids overfitting, and a higher C allows the fine-tuning of the data model. Usually, SVMs are used with default parameters. In RF implementation, the following parameters of the search grid were used: number of trees in the random forest (50, 100, and 1000), maximum number of levels in each tree (from 10 to 120 with step 12), minimum sample number to split a node (2, 6, and 10), minimum sample number that can be stored in a leaf node (1, 3, and 4), and using bootstrap for data points sampling. The number of cross-validation runs (RandomizedSearchCV) was set to 10. The Catboost parameters grid: learning rate (0.03, 0.1), depth (2, 4, 6, 10), and l2_leaf_reg (1, 3, 5, 7, 9) were optimized by the build-in catboost.randomized_search procedure. A complete list of the optimized values of these parameters is shown in Table S2, Supplementary Materials.