ANALYSIS OF PHASE WAVES IN THE ECOG DATA∗

Subdural ECoG data recorded from the matrix of electrodes during syllable pronunciation are analyzed by the method of circular-linear regression. Phase waves in 1D electrode arrays and in the whole 2D set of electrodes are detected, and their spatial organization and temporal evolution are studied. Phase portraits of wave vectors indicate the presence of sources, sinks, and saddle points. The analysis of temporal evolution of phase portraits shows that they changed more at the beginning of syllable pronunciation. Furthermore, wave sources were more stable in their localization during the pronunciation. Overall, in spite of large variability of phase portraits, they represent some characterization of the dynamics of electric potential in the cerebral cortex. Mathematics Subject Classification. 92-10. Received July 30, 2021. Accepted August 26, 2021.


Introduction
Speaking appears to be an easy and natural task. When performed smoothly and successfully, it requires precise coordination of many articulators (lips, jaw, tongue, larynx). In addition, it involves brain control mechanisms characterized by complex spatiotemporal dynamics in order to retrieve the right words from a mental lexicon and to pronounce 2 or 3 words per second [20].
Methods used for the investigation of these brain mechanisms include different techniques such as grid or strip ECoG (ElectroCorticography), depth LFP (Local Field Potentials), MEG (MagnetoEncephalography) and scalp EEG (ElectroEncephalography) [6]. In addition, recent technological and theoretical advances in several fields have enabled the recording, analysis, mapping and network characterization of cerebral activity with an unprecedented level of power and details. Today, the relation between brain structure and function during language could be described in terms of complexity, dynamics and networks [1,20].
Results obtained using these recording methods suggest that language, like other behaviors, sleep states or even levels of arousal are supported by distinct patterns of brain rhythms. These rhythms vary from very low to very high frequency and they interact with each other (for a review see [5]). According to this work, at least 10 distinct mechanisms contribute to brain waves "including variations in the concentration of extracellular neurotransmitters and ions, as well as changes in cellular excitability" (p. 1). Furthermore, at the network level reference [5] indicates that "the synchronization of the neuron's electrical activity gives rise to rhythmic voltage

1D arrays of electrodes
This section is devoted to a technique of joint analysis of the phases of four electrodes of a single strip. The goal of this study is to check whether their phases are close to an arithmetic progression. The test developed here is elementary and, to our best knowledge, original. The method of circular-linear regression (see [7,8] and Sect. 2.2) can also be used and it will be considered for 2D sets of electrodes. It provides similar information in the form of goodness of fit R, and also additional information about the direction of propagation. However, the method to analyze the regularity of phases, presented here, has the advantage of being faster to compute and easier to scale. The term arithmetic progression of four angles will stand for a sequence of angles of type (φ 1 , φ 1 + δ 1 , φ 1 + 2δ 1 , φ 1 + 3δ 1 ). One may remark that the corresponding sequence of complex numbers e iφ1 , e i(φ1+δ1) , e i(φ1+2δ1) , e i(φ1+3δ1) , of modulus equal 1, is a geometric progression. Furthermore, the signals at any four regularly spaced points are of the form (sin (x 1 − ct), sin (x 1 + d − ct), sin (x 1 + 2d − ct), sin (x 1 + 3d − ct)) in the regime of a plane TW. This form implies that the phases of the four signals form an arithmetic progression of angles, where one denotes: Consider a sequence of four angles represented by complex numbers (z 1 , z 2 , z 3 , z 4 ) of modulus 1 . The regularity of these angles will be measured by the following function: Reg (z 1 , z 2 , z 3 , z 4 ) = . This function is invariant with respect to multiplication of all four numbers by a common factor, that is, by simultaneous rotation of all four input angles. Figure 1 shows some examples of the function Reg. Its maximum (equal 1) is achieved for an arithmetic progression of angles. Its minimum, equal −0.5, is achieved for the sequence (e The distribution of the values of the function Reg for 10 6 sequences of four points chosen randomly at the unit circle with a uniform distribution has been simulated numerically. Figure 2 presents the resulting histogram. The probabilities of values above certain thresholds are given in equation (2.2). Thus, 5% of sequences exceed the threshold 0.83 in a random choice of sequences of four numbers. The generalization of this function to a different number of points and the analysis of its properties can be the object of a more elaborate study.

Model of a plane wave
The brain displays a variety of highly nonlinear, complex dynamics across multiple spatial and temporal scales, including traveling waves. The analysis presented below is focused on detection and characterization of TW, inspired by the techniques presented in [23]. Suppose that the signal measured at a point x and filtered at frequency ω 2π has the form 3) The simplest model of TW is the model of plane wave [4]. The phase offset varies in this model as a linear function of the spatial point as In the two-dimensional case, the vector β = (β 1 , β 2 ) ∈ R 2 is called the wave vector and x = (x 1 , x 2 ) ∈ R 2 is the vector of coordinates of the recording point. In the one-dimensional case, β ∈ R is called the wave number and x ∈ R is the coordinate of the recording point on a line. µ is the phase at the point x = 0. Therefore, the plane wave is described by the expression y(x, t) = A cos(ωt + βx + µ) (1 dim) or (2.5) y(x, t) = A cos(ωt + β 1 x 1 + β 2 x 2 + µ) (2 dims). (2.6)

Relation with the properties of the wave
Consider the model of plane wave defined by (2.5) or (2.6). One may be interested in its properties such as speed, wavelength and direction of propagation in order to interpret the estimated values of parameters.
The wave vector β (or the wave number β) is measured in rad/mm or rad/electrode if measurement is done by electrodes located in nodes of a regular grid. The latter situation is the hypothesis of the present study. The phase µ is given in radians. The wave has the speed c = ω β and the direction of the wave is opposite to the direction of β. The wavelength equals 2π β .

Statistical model
This section follows the presentation of the circular-linear regression in [8]. Equation (2.4) cannot be used as a valid statistical model in this form for two reasons. First, the explained variable φ is angular: angles such as (π + ε) and (−π − ε) should be considered as close to each other for small ε. On the other hand, the independent variables x 1 , x 2 are scalar, and the scalar values (π + ε) and (−π − ε) do not converge to each other as ε tends to 0. Second, the phase φ is measured with some noise in practice. Therefore, the noise should be modeled as well.
The first issue is addressed by introducing a link function g, which takes a scalar input and returns an angular output. Following Section 3.2 of [8], the function is used. The second issue is addressed by introducing an additive angular noise, which follows the Von Mises distribution VM(0, κ) centered at 0, with constant, unknown concentration κ. The Von Mises family of distributions shares many of the key properties for statistical inference that the normal distribution has for linear data (see [8] and Sect. 3.4.10 of [12]).
The resulting model is the following: According to [8], parallels with the theory of generalized linear models are obvious. One can remark that N.I. Fisher provides an equivalent formulation of the model. The form (2.8) intends to help the reader's intuition by underscoring this analogy.

TW parameter identification
The model (2.8) is fit to the training set where x i are the scalar or vector coordinates of the electrodes, and φ i are angular variables, by the maximumlikelihood estimation of the model parameters µ, β, κ. The first two parameters are estimated using the iterative procedure, described in [8].
The estimates of the phase at the point x=0 and of the wave vector are denoted byμ,β. These quantities allow computing the goodness of fit R. The goodness of fit is introduced in the articles [8,23] and measures how close the predicted anglesφ i are to the actual phases φ i by the formula The concentration parameterκ can be estimated using R bŷ where the function A(κ) is the ratio of modified Bessel functions of orders 1 and 0: According to the values in the cited tables, higher goodness of fit R is associated with higher concentration κ, therefore with smaller predicted error . As computingκ from (2.10) is tedious, the goodness of fit is used in this study. When working with scalar or vector data, it is customary to use residual variance (or dispersion) instead of goodness of fit. However, modeling dispersion for an angular deviate presents certain difficulties because there is no natural measure of scale for circular distributions [8].

Local fitting of parameters
The model of plane wave is only valid in a homogeneous physical medium. Taking into account a possible heterogeneity of the ECoG signals, we introduce a space averaging on some subgrids of electrodes, and local values of parameters are estimated independently.
On the other hand, the technique presented here requires data measured at several spatial points. A balance should be achieved between preserving the local nature of the phenomenon and using enough data points for fitting the parameters. One-dimensional strips of size 4 × 1 and the two-dimensional subgrids of size 4 × 4 electrodes are processed in this study, and local values of parameters are estimated for each subgrid.
The data is first processed point by point. The phase information associated to each electrode is recovered by the operations of bandpass-filtering, Hilbert transform, extracting phase at a fixed point of time from the analytic signal. The training set is composed of pairs where x i are the scalar or vector coordinates of the electrodes, and φ i are angular variables.
The points x i represent a subgrid of the ECoG electrodes. Their coordinates are converted to the integer grid ( [1,2,3,4] or [1,2,3,4] 2 ) before passing to the regression. Here φ i are the phases at given time t of the analytic signals associated to the time signals, recorded at these electrodes, filtered at the given frequency band.

Interpretation of the results
This approach gives a matrix of local estimates of wave vectors for every region of size 4 × 1 or 4 × 4. These matrices characterize wave patterns at the given moment of time and at the given frequency range.

Distribution of values of the goodness of fit
Interpreting the goodness of fit R of a circular-linear regression requires information about its distribution for random phases. This distribution is computed for random angles φ i and fixed coordinates x i (values of the independent variable), in an analogous way to the distribution of regularity of phases shown in Figure 2. The number of points and the coordinates x i is chosen in order to match the practical situations of use of the circular-linear regression. A regular strip of electrodes [1,2,3,4] or a square [1,2,3,4] 2 are the two choices used in this work.
The distribution of goodness of fit for 4 random phases has been simulated for a sample of 10 5 random sequences, each of them being composed of four independent uniformly distributed points at the unit circle. The result is shown Figure 3. The 95th percentile of this distribution is close to 0.92, this value is used Section 3.2.2 for the interpretation of the results.

Speed considerations
This section presents two different techniques for analyzing an array of 4 phases: the regularity of phases (Sect. 2.1) and the circular-linear regression (Sect. 2.2).
The circular-linear regression provides additional information about the mean rate of change of phases and can be generalized to arbitrary number and positions of electrodes, as well as to 2D lattices. Due to these features, the circular-linear regression is the main tool used Section 3.3.1.
Nevertheless, the regularity of phases remains useful for checking the parameters of our implementation by comparing the results of two different methods of analysis. Such comparison is performed on real data Section 3.2.2.
As regularity is analogous only to a part of the output of the circular-linear regression, one may expect it to be faster to compute. Indeed, the regularity of phases is defined by a closed formula (2.1) while the circularlinear regression is computed by an optimization procedure. In our implementations, the regularity of phases is computed approximately 10 4 times faster. This implies that the distribution of values of regularity (Fig. 2, sample size 10 6 ) is much faster to compute, than the distribution of values of the goodness of fit (Fig. 3, sample size 10 5 ).

Subject and experimental task
We used data published by [3] and analyzed by [18]. The data was collected after the experimental protocol was approved by the Human Research Protection Program at the University of California, San Francisco. We used the data from one of the three subjects tested (EC2). This subject was a male native English speaking participant who underwent chronic implantation of a high density, subdural, electrocorticographic (ECoG) array over the left perisylvian cortex and specifically on the ventral sensorimotor cortex (vSMC) as part of his clinical treatment of epilepsy. He had a left cerebral language dominance determined by Wada Testing, his MRI was normal, he had normal fluency and no dysarthria. The subject gave his informed written consent before the surgery, had normal hearing, and underwent language testing including the Boston Naming test and had normal verbal fluency tests results. He read aloud consonant-vowel syllables (CVs), composed of one of 19 different consonants followed by one of 3 different vowels at his own pace, and had 2 rest periods during the experiment. The periodicity of reading was about 1.5 s per syllable, and the duration of rest periods about 1 minute. The CVs we studied included Lee and Ree and were repeated 6 times. Finally, microphone recordings were synchronized with the ECoG data. A total of 311 syllables were read, and their data was used for statistical analysis of the results.

Data acquisition and signal processing
ECoG data were recorded at 3052 Hz (for specific details see [3]). Time series from each channel were filtered to eliminate 60 Hz noise, artifacts were removed. The channels retained were used for spectro-temporal analysis. The onset of the consonant-to-vowel transition was used as a common temporal reference point for subsequent analyses because the authors found that this transition was highly reproducible across the syllables. We also considered the transition from silence to speech as well as a pause every 3 to 5 minutes in our analyses.

Examples of filtered signals
An example of filtered signals of four consecutive electrodes (135-138, Fig. 4) is presented in Figure 6 of [17] and reproduced (using second-order Butterworth pass-band filter, function filtfilt, MATLAB, MathWorks) in Figure 5. These ECoG data are registered during syllable pronunciation as described above. The moments of the beginning of pronunciation are shown by black vertical lines and the end of pronunciation by orange lines. Duration of each epoch is approximately 1.5 second. Low frequency signals in the band 0.4-0.8 Hz correspond to the frequency of syllable pronunciation providing activity entrained oscillations. The previous work [17], studying the same data, provides additional arguments based on the measurements. On one hand, the ITC (Inter-trial coherence) plots (see Figs. A2 to A17) show higher values in the band 0.4-0.8 Hz. Higher ITC values indicate that the phases in different trials are close to each other in the frequency band 0.4-0.8 Hz, and this observation is consistent over time and electrodes. On the other hand, the analysis of phase-amplitude coupling shows that high-gamma amplitude oscillations specific for cognition are coupled precisely with low frequency phase oscillations 0.4-0.8 Hz.
One can remark that the maxima of these four signals are reached near the beginning of syllable pronunciation (black lines), while their minima near the end of pronunciation. If we consider electrode 136 providing the maximal amplitude of oscillations (green line) then its maximum is located just before the beginning and the minimum just before the end of pronunciation for all seven epochs shown in the figure  The maxima of these four signals, as well as the whole curves, are shifted with respect to each other due to different spatial positions of the electrodes. This phase shift with ordering 135,136,137,138 is a manifestation of a periodic TW. This question will be discussed in more detail in the following subsections.

Regularity of phases and circular-linear regression
In this subsection and the rest of the article, the phases are extracted from the filtered signals using the Hilbert transform. The oscillation phase θ(t) is extracted from the analytic signal A(t)e iθ(t) .
The regularity of phases can be computed for the four signals according to the description given in Section 2.1.  Figure 6. These observations are compatible with the persistence of a periodic TW during syllable pronunciations and during the pauses between them. The origin of low frequency fluctuations of the regularity of phases is not clear. It can be related either to a difference in articulation or to the interference with other brain activities.
The circular-linear regression (Sect. 2.2) can be applied to the 1-dimensional array of phases of the four signals. Figure 6 shows the goodness of fit R. Similar to the regularity of phases, it has a high value throughout the whole time interval exceeding the value 0.72, which exceeds 74% of the random values (cf. Fig. 3). As before, there is a strongly pronounced minimum between 396 s and 397 s, and a similar tendency of alternation of other maxima and minima.
In order to check this qualitative similarity, the linear correlation between the regularity of phases and the goodness of fit was computed on this data. The array of 201 time points, regularly spaced between 390 s and

Determination of 2D phase waves
The techniques of analysis based on studying a strip of electrodes in one dimension (Sect. 3.2) implies an a priori selection of the electrodes. The grid of 16 × 16 electrodes (Sect. 3.1) contains many possible strips. There are two limitations of the 1D approach. First, a signal may propagate by a path, which changes direction. Such paths become difficult to detect. Second, the available directions, therefore the results, change after rotation of the grid of electrodes. This creates an unwanted dependency of the results on the properties of the measuring device.
In order to overcome these limitations, it is possible to use the circular-linear regression (see Sect. 2.2) on squares or rectangles of electrodes. Such regression takes a matrix of phases of electrodes as input and produces a wave vector β and a scalar goodness of fit R.
Contrary to the analysis Section 3.2, the present approach deals with the phases of the whole ECoG grid at some given moment of time. We call phase portraits the ensemble of wave vectors obtained by the circularlinear regression (Fig. 7). We will analyze the evolution of the phase portraits in time in relation with syllable pronunciation.
The regressions at a single time point are visualized as a diagram, where the electrodes are represented as the squares of a grid of size 16 × 16, and the regression at each square subgrid of electrodes is represented by an arrow. These arrows are drawn according to the following rules.
-Each arrow starts at the center of the corresponding subgrid.
-The direction of each arrow is that of the vector − sgn(sin( β ))β. It is opposite to the vector β if β < π corresponding to the interpretation in Section 2.2.2 (the condition is satisfied in most cases in practice). -All arrows have equal lengths.
-The color of each arrow depends on the speed of the wave estimated by the formulaĉ = ωd | sin( β )| .
The correction factor sgn(sin( β )), used in visualization, is a simplified way of dealing with large wave vectors β. It is compatible with β being measured in rad/electrode (see Sect. 2.2.2), which means that the unit of length is the smallest distance between electrodes. Therefore, in a simplified 1D case, a model withβ = β 1 is equivalent to a model withβ = β 1 + 2π, in the sense that both models predict equal phases at the locations of electrodes. The correction factor leads to visualizing such equivalent circular-linear models by identical arrows.
It should be noted that a large wave vector β may lead to an inverted direction of change of the phases of the electrodes. This rarely happens in practice and leads to an inversion of the inferred direction of the wave.
The size of the subgrids for the 2-dimensional regression can be set as parameter. This choice is discussed in Section 2.2.5 and illustrated in Figure 7. One can remark that the smaller subgrids lead to more detailed results, while the bigger subgrids lead to smoother results. The influence of this parameter is similar to the sharpness of an image.
The size 4 × 4 is the main choice for this study. This size is suitable for the purposes of the paper since it allows detecting the structure of the phase portrait (see next section) and provides more detail than the grid size 5 × 5.

Phase portraits
The part of the cortex covered with the electrodes represents a square with approximately 6 × 6 cm size. We can introduce spatial coordinates (x, y) in this square with the point (0, 0) in the lower left corner. We suppose that electric potential P (x, y, t) is a continuous function of the space variable (x, y) and of time t together with its first partial derivatives. We present a simplified model, which does not take the sulci and gyri into account, and supposes that the cortical waves propagate in a homogeneous medium (see also the analogous hypotheses for the circular-linear regression, Sect. 2.2.5). Under this assumption, the phase Φ(x, y, t) is also a continuous function with continuous first partial derivatives.
Consider the level line Φ(x, y, t 0 ) = Φ 0 for some fixed t 0 and Φ 0 . Then the gradient ∇Φ = (∂Φ/∂x, ∂Φ/∂y) is perpendicular to this level line. It corresponds to the wave vector (Sect. 2.2.1) and it can be determined by the method of circular-linear regression. Hence, we determine the vector field and the corresponding dynamical system at each spatial point (x, y). Intrinsic time τ of the dynamical system is not related to the physical time t which should be considered in the ODE system (3.1) as a parameter. Let us note that the method of optical flow also allows the construction of phase portraits [11,21]. Singular points of the dynamical systems are determined by the equality ∇Φ(x, y) = 0 (parameter t here is omitted for convenience). The minima of the function Φ(x, y) correspond to stable nodes, and the maxima to unstable nodes, and saddles to saddle points. We will interpret below unstable nodes as sources of waves.
If we assume that the function Φ(x, y) has continuous second partial derivatives, then we can consider the corresponding Hessian matrix H. Since it is symmetric, then it has only real eigenvalues. Therefore, its non-degenerate singular points can be only nodes or saddles. Focuses and centers do not exist for such systems.  A discrete set of electrodes provides a discretization of the vector field (3.1) by vectors ∇Φ(x i , y i ), where (x i , y i ) are spatial coordinates of the electrodes. Figure 8 shows two examples of such discrete vector fields. Clearly, the analysis of such vector fields encounters some difficulties due to (a) the discrete vector field can be completed to a continuous vector field in different ways, (b) the fluctuations in the experimental data. We will discuss these two questions.

Discrete and continuous vector fields
We suggest the following way to complete a discrete vector field to a continuous one. Consider two vectors Φ i and Φ j located at two neighboring discretization points (x i , y i ) and (x j , y j ), respectively. We complete the vector field on the straight interval I ij connecting these two points by a linear interpolation of the vectors Φ i and Φ j in the direction of the minimal angle between them. For example, if Φ i = (1, 0) and Φ j = (0, 1), then we consider unit vectors between 0 and π/4 counterclockwise and not clockwise. Certainly, it is an arbitrary choice of interpolation, but it is natural, especially if the difference between the two vectors is small.
At the next stage of transition from a discrete to a continuous vector fields, we need to consider a unit square with four vectors at its vertices. First of all, we complete the vector field at its sides, as described above. Next, we should continue the vector field from the boundary of the domain inside the domain. As before, there is no unique way to do it but the winding number (rotation of the vector field at the boundary) is determined by the vector field on the boundary of the domain if it does not vanish there. In particular, the existence of singular points inside the domain depends on the winding number. If there are no singular points, then it equals 0; if there is a single node, then it equals 1; if there is a single saddle point, then it equals −1. Taking into account that the vector field in the unit square is determined by four vectors at the vertices with the method of completion described above, the winding number can adopt only these three values. We should note that a different number of singular points can exist for the same winding number. For example, two nodes and one saddle point give also 1, while one node and one saddle point give 0. We will assume that only one singular point exists for ±1 and no singular points for 0 since there is a continuous deformation of the vector fields inside the unit square reducing it to these three canonical cases. Moreover, since there are no data inside the vector field inside the unit square, we will consider the simplest one. Finally, since our purpose in what follows consists in the determination of singular points, and they are determined by the four vectors at the vertices, we do not need to construct a continuous vector field explicitly.

Examples, time evolution and fluctuations
Some examples of vector fields obtained for the ECoG data are shown in Figure 8. There are four saddle points and three nodes in the left image. The node in the upper left corner (circle) is clearly unstable. Two other singular points surrounded by closed curves might be focuses. However, due to the analysis presented above, the model defined by the system (3.1) does not contain focuses. Therefore, these singular points should be treated as nodes.
The phase portrait in the right image is obtained from the ECoG data 0.55 s of physical time later. Hence, we can consider it as a continuous deformation of the previous phase portrait. They are quite similar to each other but one of the two nodes in the left image is not present any more in the right image. According to the behavior of dynamical systems, it moved across the boundary of the square domain outside. Indeed, a singular point with a nonzero index (winding number) like a node cannot disappear remaining inside the domain. Either it moves through the boundary or it merges with some other singular point of an opposite index. A node having index +1 can only merge with a saddle point having index −1. However, all saddle points remain the same though the upper left saddle point changed its location.
Thus, the analysis of phase portraits allows us to make the following conclusions. First of all, they are relatively stable during sufficiently long time intervals (physical time). This stability can be partially related to the low frequency band considered here. On the other hand, quite rapid changes of the phase portraits are observed with singular points continuously moving, crossing the boundaries or bifurcating inside the domain (not shown). In particular, this means that singular points are not attached to some particular brain sites but emerge as a result of dynamic brain activity.
ECoG data can have strong fluctuations determined by various reasons. First of all, it is related to the complexity and irregularity of the brain functioning involving simultaneously numerous processes, and many of them cannot be controlled during the denomination experiment. Next, registering the data with ECoG electrodes can involve additional perturbations and imprecision. Finally, the treatment of these data (Hilbert transform, regression) implies some changes in the data representation. We cannot estimate whether this stage can increase or decrease existing fluctuations. We can expect however that spatial averaging in the construction of phase portraits serve to smoothing the data (Sect. 3.3.1). Nevertheless, analyzing the phase portraits we can observe some phase vectors falling out of the surrounding vector field. Therefore, we may need to effectuate some adjustment when we determine singular points and other properties of the phase portraits.

Waves and syllables
Let us recall that the ECoG data presented here are obtained during the experiments on syllable pronunciation. We have discussed above (Sect. 3.2.1) that periodic wave in the array of aligned electrodes can adjust to different phases of pronunciation. In this section, we will study the relation between pronunciation and registered signals in more detail.

Phase portraits
We begin with phase portraits obtained by circular-linear regression in the ensemble of electrodes. Figure 10 shows two phase portraits corresponding to the beginning of two pronunciations of syllable "Lee". These two phase portraits are close to each other except for the parts inside the rectangles. Phase portraits corresponding to different pronunciations of the same syllable are, in general, quite different from each other. We will introduce a method to characterize the level of their similarity below (Sect. 4.3).

Sources of waves
The diagrams resulting from circular-linear regressions can be analyzed automatically in order to identify some properties of phase portraits. The first method of automatic treatment consists in identifying the wave sources as unit squares of the grid such that the arrows at all four vertices are directed outside the square. In particular, this is the case of the unstable node (cf. Fig. 9, left). Let us also note that such points correspond to the maxima of the phase distribution function.
Phase portraits are constructed every 0.05 s and wave sources are automatically detected. An average density of such sources during one syllable pronunciation (speech) is shown in Figure 11 (left). The image on the right contains the densities of wave sources during a time interval, located in the middle of a longer pause in the experiment.
Let us note that the total density of sources is not significantly different between reading and speech, as shown by comparing the data of all 311 syllables. However, their distribution is different. If we count the unit squares with the frequency of sources above 0.85 (yellow in Fig. 11), then there are more of them during pronunciation than during silence (reading) (Fig. 11, left and middle). They are even less frequent during a longer pause (Fig. 11, right). Thus, the spatial position of sources is more stable during a period of pronunciation, than during a period of reading, but the sources can change their location for different pronunciations of the same syllable. Figure 11. Average density of wave sources during pronunciation of syllable "Lee" (left), during a period of silence between two pronunciations (middle), and during a longer pause (right). The squares represent the electrodes, which have at least two other electrodes on each side. This condition is necessary for being able to define arrows at all four vertices of a square, which are used in the definition of a wave source (see the beginning of this section). The positions of the electrodes are the same as at Figure 4. The color code corresponds to the average density. It equals 1 (bright yellow if the corresponding unit square represents a source for all phase portraits during considered time interval. The value 0 (dark blue) corresponds to the case where this square does not represent a wave source for any phase portrait. The sources are more stable during speech, than during reading: there are more yellow squares in the left image than in the middle image, and even more than in the right image.

Similarity
Another method of automatic analysis of phase portraits consists in comparing the results of regressions computed at two different time instants. The comparison is done using the scalar products of unitary vectors, located at same points in the diagrams. Denote by u i,j the normalized wave vectors obtained by the circularlinear regression at some moment of time t 1 and v i,j at time t 2 . Their scalar product (u i,j · v i,j ) adopts a value between −1 (opposite directions) and +1 (same directions). Pairwise comparison of these vectors for all i, j gives a 13 × 13 matrix with the elements between −1 and 1. This matrix characterizes the similarity of two phase portraits. Figure 12 shows a pairwise comparison of phase portraits in seven moments of time. It is a 7 × 7 matrix symmetric with respect to the main diagonal. Time T = 403.63 corresponds to the beginning of pronunciation of syllable "Lee" (one of six pronunciations of the same syllable). The line where this time is indicated shows the comparison with six other times before and after the beginning of pronunciation. The square at the diagonal corresponds to the comparison of the phase portrait for T = 403.63 with itself. Each comparison (square) contains 13 × 13 numbers between −1 and 1 representing a scalar product between two unit vectors. These numbers are shown with a color code where bright yellow corresponds to 1 and dark blue to −1.
The crossed square at the diagonal is composed of numbers 1. It is entirely yellow (color not shown). The squares in the same line show the comparison of the phase portraits at T = 403.63 with six other phase portraits at some moments of times before and after the beginning of pronunciation. Blue spots in these squares show the difference of the phase portraits. Considering three preceding moments of time we observe two persisting blue spots surrounded by red circles. Therefore, there is some stable structure existing at least 0.3 s before the beginning of pronunciation and disappearing when the pronunciation starts. On the other hand, there are some persisting structures appearing at the beginning of pronunciation and remaining during pronunciation. They are shown by red circles in the corresponding column. The analysis of the blue spots suggests that the changes are gradually accumulating and not changing back and forth. Thus, the beginning of pronunciation provokes a transition between persisting phase portraits.  It should be noted that the changes in the phase portraits are not essential, and they remain sufficiently close to each other. Figure 13 shows phase portraits at three consecutive moments of time. There are three circled areas in left phase portrait (T = 403.53) which show the main differences with the middle phase portrait (T = 403.63). They are in agreement with the blue spots in the corresponding square in Figure 12. Similarly, the right phase portrait (T = 403.73) is compared with the middle one. Figure 14. Similarity curves for six pronunciations of syllable "Lee" (left). Each point at the similarity curves is obtained as a normalized sum of scalar products of wave vectors shown in Figure 12 by the color code. An average similarity curve for 311 syllables pronounced by one subject (right). Time is shifted so that time 0 corresponds the beginning of syllable pronunciation (but not normalized). Comparison is done between current time and time 0. The value of the similarity curves equal 1 at time 0 because the phase portrait is compared with itself. Right: each vertical bar shows the confidence interval for the mean (confidence level 95%).
The overall similarity of two phase portraits can be measured by a single number equal the normalized sum of 13 × 13 scalar products. In this representation, two identical phase portraits have similarity 1, randomly generated phase portraits have average similarity 0. Starting form the moment of the beginning of pronunciation (time 0), moving forward and backward, we obtain similarity curves for each syllable pronunciation. Figure 14 (left) shows such curves for six pronunciations of syllable "Lee". Clearly, all these curves have a maximum equal 1 at time 0, and they decay in both directions.
The properties of the similarity curves become clearer for the average similarity curve for all 311 syllables pronounced by one subject (Fig. 14, right). It has a maximum at time 0, fast decay around t = 0, and slower linear decay for larger time difference. We note that this curve is nearly symmetric with respect to zero time. Fast decay near time 0 is observed for the majority of syllables, though this effect can be more or less pronounced. Linear similarity curves outside of switching corresponds to a relatively slow change of phase portraits.
The analysis shown Figure 14 uses the start of syllable (t = 0) as an a priori reference point. One may therefore be interested in having an objective criterion without choosing a reference time point. Such complementary analysis is achieved by comparing the phase portraits at consecutive moments of time (t − ∆t, t + ∆t) for all t (and a fixed half-difference ∆t). The value ∆t = 0.025 s corresponds to time step 0.05 s between phase portraits considered above (see Sect. 4.2). The values at -0.025 s , +0.025 s equal the values of the mean correlation between the time 0 and the adjacent time moments (see Fig. 14, right).
The results are shown Figure 15. We can see that the mean correlations stay relatively close to each other before the time t = 0 (start of speech). This behavior changes after a minimum at time zero, and the mean similarity becomes significantly higher. This indicates that, on average, the changes in the phase portraits happen faster during reading, than during speech. The moment of start of speech corresponds to a slightly higher rate of change (lower correlation), than during both conditions of reading or speech.

TW and their role
Propagation of traveling waves of electric activity in the brain cortex has been known for a long time (see [10] and the references therein) but their functional role is still being actively discussed [13]. In spite of some confirmation about their correlation with memory tasks [23], and some other activities [15], this question is far from being elucidated. Since the discovery of nerve impulses propagating in neuron axons and the mechanism of this propagation based on opening and closing of transmembrane ion channels, it became clear that nerve impulses play a fundamental role in the brain functioning. At the level of tissue, the mechanisms of signal propagation are multiple [6], but it can be expected that they also serve to activate or to inhibit adjacent or distant parts of the brain. However, there is a gap in the understanding of signal propagation in the individual axon and in ensemble of neurons at the mesoscopic or macroscopic levels.
It is now well established that language is a distributed brain function based on coordinated networks and subnetworks activated at the appropriate stage of visual processing [9]. The coordination of different parts of the brain involved in this process requires timely communication between them. This communication can be effectuated along the white matter fibers, or along the cortex, or, more likely, both of them. Involvement of the brain cortex in the process of communication is indirectly confirmed by stroke-induced tissue damage, which may be restricted in some cases to the cortex grey matter outside of the regions processing language. However, this tissue damage can disturb network functioning (connectome), that is, the communication between different brain regions.

Methods of analysis of phase waves
Having different spatial location, EEG or ECoG electrodes register signals with some phase difference which can be measured with respect to some reference electrode, to an average signal, or to the maximal amplitude of the given signal. Linear regression for the phase distribution allows the determination of phase waves [17]. Taking into account phase periodicity, circular-linear regression can be more appropriate in order to avoid the discontinuity of the phase distribution function. We use this method in order to construct instantaneous phase portraits of brain activity every 0.05 s. Visual analysis of these phase portraits is difficult because of their large number. We developed here two methods of their automatic analysis in order to determine waves sources and similarity of phase portraits. Let us note that phase portraits of EEG/ECoG signaling can be also constructed with the optical flow method [11,21] or using the gradients of phase [14].

Wave sources
The method suggested in this work suffers from two limitations. Wave sources are defined as unit squares for which all four wave vectors at their vertices are directed towards the outside of the square. This definition is appropriate for unstable nodes and focuses. However, saddle points or unit squares without singular points can also satisfy this definition depending on the geometry of the vector field. A more precise method requires a completion of the discrete vector field to a continuous one, as discussed in Section 3.3.2. However, this approach is more technically difficult, time consuming, and also allows some uncertainty. Another drawback of this method is that wave sources are likely to change if the grid of electrodes is displaced or if its geometry changes. This dependency of the results on the properties of the measuring device is a common property shared with the regressions in one dimension. One may remark that concepts analogous to wave sources are introduced in the works [14,21].

Similarity
Similarity analysis provides a point-wise comparison of wave vectors. Though it is sensitive to the fluctuations of phase distribution, it allows the determination of persistent structures in phase portraits. The average similarity curve for the ensemble of syllables possesses some interesting properties. First of all, it is practically symmetric with respect to the beginning of pronunciation showing that the changes of phase portraits before and during pronunciation have the same rate. Next, growth-decay rate of the similarity curve is linear outside of small neighborhood of the switching point (time 0). Analysis of similarity diagrams (Fig. 12) shows that this behavior corresponds to a gradual accumulation of changes rather than to a random fluctuations. The complementary analysis of correlations between phase portraits in successive moments (Fig. 15) indicates that the changes happen more often, in average, during reading, than during pronunciation.

Speech characterization
Language brain function is characterized by complex spatiotemporal dynamics in EEG or ECoG data. Its relation with TW is not yet well established. In this work, following [17], we investigate phase waves during syllable pronunciation. We showed that phase portraits of brain activity have some correlation with pronunciation. Namely, wave sources are more stable during pronunciation than during silence. Furthermore, the similarity analysis indicates the beginning of pronunciation makes a switch toward slower-changing phase portraits. On the other hand, there is no one-to-one correspondence between phase portraits and syllable pronunciation. For the same syllable, the corresponding phase portraits can have high or low similarity. The latter can be explained by numerous other simultaneous brain activities imposing large fluctuations on phase distributions. Therefore, other methods of characterization for the phase and for amplitude should also be used.

Trajectories
The phases at the beginning of several consecutive syllables may be close to each other for certain electrodes, as was shown in the example Section 3.2. When this happens, the mean phases over the sequence of syllables can be compared among electrodes in order to detect a path with monotonic evolution of phases (i.e., a TW in the "typical" syllable). This perspective has been studied in [18]. The same sequence of 37 consecutive syllables as in [18] was analyzed. The phases were extracted from the signal using the filter-Hilbert method.
A non-rectilinear path of 13 electrodes can be found such that the mean phases are decreasing, and the phases at each electrode are concentrated around the mean (according to the criterion of [18]). Surprisingly, a time-shift does not preserve this path: the mean phases of the same electrodes at the moments (start + 0.05 s) or (start − 0.05 s) are no longer monotonic. This result can be interpreted as a long TW, which appears in the "typical" syllable at the moment of start of speech, and is specific to this moment.
The analysis of mean phases, presented here, is an extension of the one-dimensional TW analysis (see Sect. 3.2), which considers a mean of several trials, and takes the inter-trial variability into account. The mean phases at the beginning of several consecutive syllables replace the instantaneous waves of filtered signals. A longer non-rectilinear path replaces the path of four electrodes. A weaker criterion of monotonically decreasing phases replaces the regularity of phases. The inter-trial variability is used here as a thresholding criterion for accepting or rejecting an electrode.
Both approaches suffer from a common limitation: too many paths exist in the 2D grid and their analysis is time consuming.
Taking the inter-trial variability into account in the analysis of 2D waves is more challenging. It can be done on the level of final results (directions estimated by the circular-linear regression), in an analogous way to Figure  1-figure supplement 6 of [14].

Limitations and perspectives
Registered ECoG signals represent an enormous amount of data obtained from 256 electrodes with frequency 3052 Hz for 3 subjects pronouncing more than 300 syllables. Analysis of these data should be carried out in different frequency bands, which adds one more dimension to the problem. However, the main difficulty of the analysis of these data is not even their quantity but their variability. One of the main conditions of successful modeling and explanation of some phenomenon consists in its reproducibility. However, analyzing the data obtained for the same subject pronouncing the same syllable, we often obtain different results. This is related to the intrinsic variability of brain functioning and, possibly, to the missing part of information, taking into account that only a relatively small area 6 × 6 cm 2 is covered by the electrodes, not even mentioning more profound brain structures. Therefore, our goal is to find some hidden regularities in the ECoG data characterizing brain functioning during syllable pronunciation. It should be noted that the area of size 6 × 6 cm 2 was determined by the location of the grid during the ECoG experiment.
In this work we analyze traveling waves of electric potential which can complement existing studies of sources and networks [9]. Moreover, we consider here only phase or periodic waves, contrary to the amplitude waves describing propagation of activation along the brain tissue. Let us note that, following [17], we consider in this work low frequency range 0.4-0.8 Hz imposed by the periodicity of syllable pronunciation. Higher frequency α and γ oscillation specific for cognitive tasks require further investigations.
Overall the results obtained corroborate those of [17][18][19] but the analyses presented here have an additional focus which is to characterize the dynamics of the cortical electrical potential during syllable pronunciation Analysis of phase portraits constructed by the method of circular-linear regression allows us to characterize syllable pronunciation by wave sources, which are more stable during pronunciation, and by the similarity curves manifesting switching behavior at the beginning of pronunciation. These properties can be observed for the individual syllable pronunciation and for the averaged data. However, they cannot be used to identify the syllables. In the other words, we cannot make a correspondence between each syllable and the properties of phase portraits. If we take, for example, average density distribution of wave sources for six pronunciation of syllable "Lee", it will be different from such distribution for syllable "Ree". From this point of view we can make a difference between them. However, taking an individual density distribution, we cannot classify it to one of these classes, at least for such a limited number of training set. The question of classification of individual syllables with ECoG data will require much more investigations and new approaches to data analysis.
The presented methods can be used for EEG data as well. Indeed, [16] successfully applies a simpler, but analogous method (linear regression in 1D) to lines of 6 electrodes in an EEG grid of 26 electrodes in total. A higher number of electrodes is required for running a 2D circular-linear regression: 16 electrodes in a single region. Therefore, a dense EEG grid is preferred. The method of circular-linear regression can be applied either to irregularly placed electrodes in 2D, or to reconstructed EEG sources in 3D.