## 2 Materials and Methods
### Data Analysis
These differences in coherence coefficients were averaged over the ensemble of all subjects and presented as topographic maps. On these maps, the median values of the obtained differences were marked in color. The value of the phase coherence coefficient can quantitatively characterize the non-directional interaction of the systems that generated the analyzed signals. To analyze the structure of directional couplings, we used the phase dynamics modeling method [50]. This approach involves approximation of the increments in the instantaneous phases \(\phi_{1}(t)\) and \(\phi_{2}(t)\) of two analyzed processes:  
\[\varphi_{1,2}(t+\tau)-\varphi_{1,2}(t)=F_{1,2}(\varphi_{1,2}(t),\ \varphi_{2,1}(t-d),\ \mathbf{a}_{1,2})+\varepsilon_{1,2}(t), \tag{2}\]  
where \(\tau\) is the phase increment interval, \(F_{1,2}(t)\) are trigonometric polynomials of the third order, \(\mathbf{a}_{1,2}\) are vectors of their coefficients, \(d\) is a trial delay taking into account the time delay in couplings, and \(\varepsilon_{1,2}(t)\) are model residuals. The parameter \(\tau\) in (2) was equated to one characteristic period of oscillations in the corresponding frequency range, which was calculated as the reciprocal of the frequency corresponding to the center of the analyzed frequency range: 10.0, 4.0, 0.6, and 0.2 s for \(\delta_{1}\)-, \(\delta_{2}\)-, \(\delta\)-, and \(\theta\)-ranges, respectively. The coefficients \(\mathbf{a}_{1,2}\) were estimated from time series \(\phi_{1}(t)\) and \(\phi_{2}(t)\) using the least squares method. The trial delay \(d\) was varied within one characteristic oscillation period. Based on [50], the indices of directional coupling were calculated as follows:  
\[c_{1,2}^{2}(d)=\frac{1}{2\pi^{2}}\int_{0}^{2\pi}\int_{0}^{2\pi}(\partial F_{1,2}(\varphi_{1,2}(t),\ \varphi_{2,1}(t-d),\ \mathbf{a}_{1,2})/\partial\varphi_{2,1})^{2}d\varphi_{1}d \varphi_{2}. \tag{3}\]  
These indices characterize the influence of the phase of the first signal on the phase of the second signal. For each pair of EEG lead signals, the indices (3) were calculated in \(\delta_{1}\)-, \(\delta_{2}\)-, \(\delta\)-, and \(\theta\)-ranges. For further analysis, we used the maximum values of these indices under variation of \(d\). The absolute values of the directional coupling indices depend on the individual characteristics of the signals, and their quantitative interpretation is not trivial [50]. Therefore, we used the indices (3) to analyze the qualitative changes in couplings between EEG leads at different stages of the experiment. For each pair of EEG leads, we compared the distributions of directional coupling indices at stages S1 and R1, and stages S2 and R2 using the records of all subjects. If these distributions were significantly different at the stages of stress and rest based on the Mann-Whitney U-test (\(p<0.05\)) [51], then the distribution of individual differences between the directional coupling indices was calculated. Further, a conclusion was made about the sign of the median of this distribution, indicating an increase or decrease in the coupling during stress-inducing cognitive tasks. The results of this analysis were presented in the form of an EEG lead map.