Analysis of Biomedical Signals by Flicker-Noise Spectroscopy: Identification of Photosensitive Epilepsy using Magnetoencephalograms

The flicker-noise spectroscopy (FNS) approach is used to determine the dynamic characteristics of neuromagnetic responses by analyzing the magnetoencephalographic (MEG) signals recorded as the response of a group of control human subjects and a patient with photosensitive epilepsy (PSE) to equiluminant flickering stimuli of different color combinations. Parameters characterizing the analyzed stochastic biomedical signals for different frequency bands are identified. It is shown that the classification of the parameters of analyzed MEG responses with respect to different frequency bands makes it possible to separate the contribution of the chaotic component from the overall complex dynamics of the signals. It is demonstrated that the chaotic component can be adequately described by the anomalous diffusion approximation in the case of control subjects. On the other hand, the chaotic component for the patient is characterized by a large number of high-frequency resonances. This implies that healthy organisms can suppress the perturbations brought about by the flickering stimuli and reorganize themselves. The organisms affected by photosensitive epilepsy no longer have this ability. This result also gives a way to simulate the separate stages of the brain cortex activity in vivo. The examples illustrating the use of the"FNS device"for identifying even the slightest individual differences in the activity of human brains using their responses to external standard stimuli show a unique possibility to develop the"individual medicine"of the future.

, demonstrate that the real time series V(t), where t is the time, for measured dynamic variables, which characterize the current condition of a living system, on some specific time intervals T usually contain chaotic components with a "long memory". In other words, there are correlation links on large intervals Т tot of the time series. The links are commonly determined by analyzing the power spectrum S P (f), where f is the frequency, when S P (f) is described by the flicker-noise function f -n with exponent n ~ 1 [6,7]. In some cases, special "memory functions" characterizing the long-range correlations are introduced [8,9]. The long-range links virtually indicate that the living organism, which is an open nonstationary system operating under variable external conditions, can quickly and efficiently reorganize itself, thus manifesting its property of biological adaptability [8]. Specifically, the analysis of heartbeat interval series shows that the mechanisms of neurohormonal cardiac regulation are brought into effect as dynamical rearrangements in the visual chaos of the studied time series. The decrease in chaoticity and the loss of long-range correlations in the measured biomedical signals may sometimes be associated with unfavorable changes in the organism, deviation of its functioning state from the normal state, and pathological changes in the organs [2][3][4][5].
As the real signals generated by living systems contain both chaotic and regular components [1,8,9], the above conceptual conclusions may not be directly used for evaluating the state of a specific organism and the effect of various factors such as medications and stimulation on its operation. In this case, it is necessary to separate the contribution of each component from the informative medical characteristics of the system before making any conclusions about the degree of loss in the correlation links. For example, stochastic quantifiers of a statistical memory were used to describe the phenomenological regularity accounting for the fundamental role of chaoticity and robustness in the functioning of living systems [8]. Based on the fundamental laws and concepts of memory functions formalism and FNS phenomenology, an original algorithm in which the effects of dynamic intermittency, nonstationarity, and characteristic frequencies in the original time series of a pathological tremor are separated out was proposed to evaluate the effectiveness and quality of different methods of treating Parkinson's disease [9]. This separation procedure can be used to advance in solving the general problems of medical diagnosis by providing the basis for the standards of medical signals through relating the values of specific signal characteristics to the particular state of the organism [10]. It is obvious that the general analysis of such problems is complicated by the presence of individual features in each organism and the specificity of its response to various treatment techniques.
The individual features of a living organism manifest themselves primarily in the lowfrequency components of its biomedical signals, which account for the collection of characteristic and stimulus-initiated frequencies particular to each organism, as well as in the interferential contributions of these resonances. In this case, the low-frequency "envelopes" are always accompanied by high-frequency chaotic ("noise") components the series of which is used to identify the informative correlation links individual to each organism. The state of the living organism exposed to external, including therapeutic, stimuli and the dynamics of its subsystems can be adequately evaluated only by sequential separation of the contributions of these components in the physiological time series for different frequency bands and introduction of the corresponding parameterization.
The phenomenological scheme for representing the information stored in various complex signals developed by Nicholis [11] provided a conceptual foundation for making some progress in the practical extraction of both chaotic and resonant components from medical time series. This scheme assumes that there is an infinite number of levels in the evolution hierarchy of the system under study and that there are some recurrent rules that generate information on a specific hierarchy level and compress it on a higher level. The general ideas of Nicholis' scheme were used to develop flicker-noise spectroscopy (FNS) [12,13], the phenomenological framework in which the concept of the information contained in the signals generated by open dissipative systems is generalized. According to the basic idea of FNS, the correlation links existing in the sequences of different irregularities such as spikes and "jumps" in original signals and discontinuities in their derivatives of different orders on all levels of the spatiotemporal hierarchy of the system under study are treated as main information carriers. In this theory, the tools for the extraction and analysis of information are the cosine-transform S(f) of an autocorrelation function, where f is the frequency, and the difference moment, or structural function, Φ (2) (τ) of the second order, where τ is the time lag, which complement each other in the extraction and analysis. The chaotic component Φ c (2) (τ) of the structural function is formed exclusively by the "jumps" of the dynamic variable while the chaotic component S c (f) of the cosine-transform is formed by both spikes and "jumps" on every level of the hierarchy. The information "passport" characteristics determined by fitting the derived expressions to the experimental curves of S c (f) and Φ c (2) (τ) are interpreted as the correlation times and parameters describing the rate of "memory loss" on these correlation time intervals for different irregularities. The number of extracted parameters depends on the problem under study [12][13][14][15]. simultaneously measured at different points in space. The FNS approach has already been applied to the analysis of the structure and dynamics for various physicochemical, electrochemical, biological, geophysical, and astrophysical processes [12][13][14][15].
The major problem emerging in the analysis of biomedical signals is its ability to provide an adequate evaluation of the dynamic states of the organism. As the spatiotemporal structure of the signals of a living system generated as the response to given stimuli is unique, the adequacy of the analysis can be provided only when the method in use makes it possible to identify even the slightest individual differences in the measured signals. In the present paper, the study of neuromagnetic responses by analyzing the magnetoencephalographic (MEG) signals recorded as the response of a group of control subjects and a PSE patient to equiluminant flickering stimuli of different color combinations [16][17][18][19] will show that the FNS method can be used to identify even the slightest differences in the recorded individual responses.

BASIC PRINCIPLES OF FLICKER-NOISE SPECTROSCOPY
(1) A hierarchy of spatiotemporal levels in complex open dissipative systems whose chaotic evolution is described by measured dynamic variable V(t) on time interval Т is introduced.
(2) The main information hidden in stochastic signals V(t) is provided by specific "resonant" components as well as by sequences of different types of irregularities such as spikes and jumps in the original signals and discontinuities in their derivatives of different orders on all levels of the spatiotemporal hierarchy of the system under study. In FNS, all the introduced information is related to one of the fundamental concepts of statistical physics, the autocorrelation function where τ is the time lag parameter (0 ≤ τ ≤ T/2). This function characterizes the correlation of the values of dynamic variable ( ) V t at higher and lower values of the argument. The angular brackets in relation (1) stand for the averaging over time interval T : The averaging over interval T implies that all the characteristics that can be extracted by analyzing ( ) ψ τ functions should be regarded as the average values on this interval. If the interval T is a section of the larger interval tot T (T < tot T ), the value of function ( ) ψ τ can depend on the position of the interval T within the larger interval tot T . If there is no such dependence and ( ) ψ τ is a function only of the difference of the arguments of the dynamic variables involved in (2), the evolution process being analyzed is defined as stationary. In this case, ( ) ψ τ = ( ) ψ −τ .
(3) To extract the information contained in ( ) ψ τ ( ( ) 0 V t = is assumed), one should analyze the transforms, or "projections", of this function, specifically the cosine transform (power spectrum function) ( ) S f , where f is the frequency: and its difference moments (transient structural functions) of the second order Φ (2) (τ): It is obvious that for the stationary process we have [ ] implying that Φ (2) (τ) depends linearly on ψ(τ). The introduction of Φ (2) (τ) as the "projection" of autocorrelator ψ(τ) is convenient because the functions Φ (2) (τ) are formed solely by the jumps of the dynamic variable at different spatiotemporal hierarchy levels of the system and S(f) is formed by both spikes and jumps.
To illustrate this statement, consider the process of one-dimensional "random walk" with small "kinematic viscosity" ν ( Fig. 1). The small value of ν implies that when the signal passes from position V i to V i+1 , which are |V i+1 -V i | apart (in value) from each other, the system first overleaps ("overreacts") due to inertia and then "relaxes". We assume that the relaxation time is small compared to the residence time in the "fluctuation" position. It is obvious that when the number of walks is large, the functions Φ (2) (τ) will be independent of the values of "inertial overleaps" of the system and depend only on the algebraic sum of walk "jumps". At the same time, the functions S(f), which characterize the "energy side" of the process, will depend on both spikes and jumps. It should be emphasized that it is the intermittent character of the evolution dynamics that accounts for the differences in the information stored in various irregularities. In other words, if there is no intermittency, the information contents of S c (f) and Φ c (τ) will coincide with each other.
(4) The "passport" parameters extracted by analyzing the functions S(f) and Φ (2) (τ) built using the time series V(t) have the meaning of relaxation times and the parameters characterizing the loss of "memory" (correlation links) at these correlation times for irregularities such as spikes and jumps. The corresponding parameters for irregularities such as "discontinuities in derivatives" are extracted from the power spectra and difference moments built using time series for the values of the dynamic variable recorded at discrete times t k .  which are obvious for continuous signals (Fig. 1) and caused by the difference in the contributions made by the high-frequency (spike) and low-frequency (jump) chaotic components into these functions, can also manifest themselves in the functions S(f) and Φ (2) (τ) build using discretely recorded time series [15]. This supports the validity of the fundamental FNS hypothesis that the difference in the information character of functions S(f) and Φ (2) (τ) can exist even when there are no "singularities" in the measured discrete signals. In the general case, when a complex stochastic signal is measured at some sampling frequency f d , we can speak only about finding the collection of parameters characterizing the correlation links in the sequences of jumps and spikes inherent to the given signal determined at sampling frequency f d .
Consequently, the sampling interval should be considered as another parameter in general parameterization algorithms for real signals.

BASIC RELATIONS
Let us write the basic interpolation expressions for chaotic components used in the analysis of experimental time series. The parameters characterizing the dynamic correlations on every level of the evolution hierarchy are assumed to be the same. Consider the simplest case, in which there is only one characteristic scale in the sequences of spikes and jumps: where Γ(s) and Γ(s, x) are the complete and incomplete gamma functions (x ≥ 0 and s > 0), respectively; σ is the standard deviation of the measured dynamic variable with dimension [V]; H 1 is the Hurst constant, which describes the rate at which the dynamic variable "forgets" its values on the time intervals that are less than the correlation time T 1 [12][13][14][15]. In this case, T 1 may be interpreted as the correlation time for the jumps in the stochastically varying time series V(t).
For asymptotic cases, we obtain the formulas: The interpolating function for power spectrum component S cS (f) formed by spikes can be written as: Here, S cS (0) is the parameter characterizing the low-frequency limit of S cS (f) and n 0 describes the degree of correlation loss in the sequence of spikes on the time interval T 0 .
The interpolating function for the power spectrum component S cR (f) formed by jumps is written as : where S cJ (0) is the parameter characterizing the low-frequency limit of S cJ (f).
Although the contributions to the overall power spectrum S c (f) given by Eqs. (9) and (10) are similar, the parameters in these equations are different: S cS (0) ≠ S cR (0) , T 1 ≠ T 0 , and 2H 1 + 1 ≠ n 0 .
This implies that the parameters in the expressions for the power spectrum and structural function of the second order generally have different information contents when the experimental time series V(t) is analyzed. For example, the characteristic times T 0 and T 1 are usually much different because they correspond to different frequency bands of the power spectrum. The fact that the jumps are more regular than the spikes implies that the contribution of jumps to the power spectrum will be concentrated in its lower frequency band. At the same time, its higher frequency band, which is often characterized by the flicker-noise function S c (f) ~ 1/f n , is generated mostly by spikes.
As noted above, the slowly varying components specific to each complex system under study, which are characterized by their own collection of frequencies, can affect the chaotic dynamics of the sequences of informative irregularities in such systems. This gives rise to interferential frequencies, which are produced by the addition of different resonant components, in the power spectra S(f) to be analyzed. During the evolution of open systems, the entire set of these resonant and interferential frequencies can be rearranged. For convenience, in our further consideration all specific frequencies manifesting themselves in the oscillations of the dynamic variable V(t) under study, no matter the genesis of these frequencies used in functions S(f), will be considered as resonant.

SEPARATION OF FREQUENCY COMPONENTS FROM A SIGNAL
In the general case, the signal V(t) under study can be formally written as where Here, V cS (t) and V cJ (t) are the chaotic components formed by spikes (mostly, the highestfrequency band) and jumps (mostly, the intermediate-frequency band), respectively; V r (t) is the low-frequency signal formed by resonant components, which are characterized by a "gradual" variation against the background of mostly high-frequency chaotic components.
Assume that there is no correlation of the high-frequency components V cS (t) with V r (t) and V cJ (t): and the later two components V cJ (t) and V r (t) can correlate with each other: In addition, it is assumed that Relations (12) Therefore, the determined value of T 0 should be corrected below. In this case, this value depends on the contribution made by the jumps in the band of intermediate frequencies.
The (3) and (12)- (14) can be written where S r (f), S cS (f), and S cS (f) are the cosine transforms corresponding to the correlators ψ r (τ), ψ rcJ (τ), and ψ cJ (τ). The function Q rcJ (f) can be used to write the "incomplete" inverse cosine transform: where f max = ½ f s = ½ ∆t -1 = N/2T; f s is the sampling frequency; ∆t is the time interval between the adjacent digitized values of signal V(t), the total number of which is equal to N. It is obvious that even if conditions (12)- (14) are satisfied, there is no way to write expression (16) as a linear combination of introduced correlators ψ r (τ), ψ rcJ (τ) and ψ cJ (τ) because of the bounded value of the upper integral limit in (16). Moreover, one can expect that the value of autocorrelator ψ cJ (τ), which corresponds to the high-frequency (against the background of resonant frequency) component of the signal, will mostly contribute to φ rcJ (τ) at small values of τ in the interval [0 ≤ τ max ≤ T/2] under consideration. At the same time, the autocorrelator ψ r (τ) and interferential cross-correlator ψ rcJ (τ), which depend on the low-frequency resonant component, can be virtually reproduced [mainly, ψ r (τ)] using transform (16) and are able to determine the function φ rcJ (τ) in the most part of interval τ. The above is supported by the further analysis of signals V(t) of different nature given in [12,14] and below.
In the parameterization procedure, φ rcJ (τ) in view of (5) and (14) is used to find the function: [ ] which characterizes the total contribution made by the resonant component and interferential "resonant-jump" terms into the difference moment of the second order Φ (2) (τ). In this case, the complete function Φ (2) (τ) can be found by adding interpolation expression (6) to (17): The analysis made in [12,14] and below shows that the use of expression (18) (17) and (18) using the method of least squares is used to determine the values of parameters σ, T 1 , and H 1 , which characterize the contribution made by jumps into the structural function. In this case, the values of parameters H 1 , σ, and T 1 are chosen by providing the best agreement between the experimental and calculated curves Φ (2) (τ) in the entire interval of τ under study. As a result, the above parameter H 1 is somewhat different in meaning from the Hurst constant, which is usually introduced for describing the functions Φ (2) (τ) at small values of τ. After that, the parameters H 1 and T 1 are used to find the interpolation function S cJ (f), formula (10), which characterizes the contribution made by jumps into the spectrum S(f) of the autocorrelator. In this case, the value of S cJ (0), like the value of S cS (0) introduced above, is regarded as a free parameter.
When the interpolation spectrum found above, is compared with the overall spectrum S(f), which was calculated by (3)

RELATION BETWEEN FNS AND DIFFUSION PARAMETERS
To conclude this section, we will discuss the relation between the introduced difference moments of the second order Φ (2) (τ), which in the general case characterize the random walks of system states (Fig. 1), and the processes of anomalous diffusion. For these processes, the average squared deviation of system states V varying over the whole set of possible states, (-∞ < V < ∞), from the average value found for the time τ can be written as [21][22][23][24][25] ( ) 1 Here, D is the diffusion coefficient; t 0 is the characteristic time; H 1 is the Hurst constant; the averaging, denoted by the symbol <(…)> pdf , is effected by introducing the probability density function W(V, t), which accounts for the probability of the system state being within the given interval of states. It is assumed that the system was found To find the relation between the phenomenological FNS parameters σ , T 1 and H 1 and the parameters characterizing the diffusion dynamics for a stationary process, we will first consider the simplest case of Fick's diffusion, for which H 1 = ½. Assume that the behavior of the probability density W(V, τ) for the random variable V on the segment [-L, + L] over time τ can be described by the diffusion equation: with the reflection conditions at the end points of the segment and the initial condition The expressions for the solution W(V, τ) to this system, the average value of random values at x <<1 and x >> 1 coincide: It is difficult to compare the functions <V 2 > pdf and Φ  (7) and (20), which correspond to the case of anomalous diffusion at small values of τ, by choosing T 1 as the characteristic time t 0 and assume that in this case the FNS parameter σ corresponds to some model parameter L a determining the region in the domain of values of the dynamic variable to be analyzed where the states of the system can be localized: Here, b is the dimensionless parameter. Below, biomedical signals will be analyzed to show how the concept of anomalous diffusion can be used to identify the individual functioning features of a human brain. ). It was found that among the sensors recording MEG response signals, sensor 10, which is localized at the frontal lobe on the head, shows the highest sensitivity to these color stimuli [18,19].
It was shown that the zones of the human cerebral cortex located in the vicinity of sensor 10, as well as some others, are the zones where the pathological PSE mechanisms are localized and can be the propagation centers for the abnormal group excitement of neurons in the cortex and subcortical structures, which ultimately causes an epileptic seizure [18,19]. Therefore, the present paper was focused mostly on the analysis of the MEG response signal recorded by this sensor. Also, to compare the results, we analyzed the signal recorded by nonspecific sensor 43, which was located on the scalp over the occipital lobe of the cerebral cortex. First, we will compare the FNS parameters characterizing the responses of the healthy subjects and patient to the RB stimulus. The differences between the recorded responses of the subjects and patient primarily manifest themselves in the spectra S(f). The low-frequency bands in the spectra of subjects 7 and 9 show a collection of specific peaks at frequencies from 2 to 30 Hz (Figs. 4b, 5b) whereas the patient spectrum shows clearly seen 50 and 100 Hz peaks in addition to the intense 2.5 Hz peak (Fig. 6b).
We can also see the difference in the characteristic curves for the components of (2) ( ) r Φ τ , which determine the total contribution of specific frequencies and the interferential "resonantjump" component into the structural functions Φ (2) (τ). There are clearly cut differences not only in the values of the parameters, mainly in Т 1 , which are two orders of magnitude lower for the patient, but also in the presence of large intervals on the abscissa axis τ with (2) ( ) r Φ τ < 0 for either subject (Figs. 4c, 5c). At the same time, the corresponding curve for the patient response signal shows (2) ( ) r Φ τ > 0 throughout the entire range of τ (Fig. 6c). In this case, it is necessary to note the values of relaxation times Т 0 found by analyzing the response signal produced by the patient, which are one order of magnitude shorter than those for the subjects.
A rather unexpected result of the analysis is the fact that the function typical for anomalous Levy diffusion can adequately approximate the chaotic (less the total contribution of resonances and interferential resonant-jump component) component in the signals produced by subjects 7 and 9 (Figs. 4d, 5d). In this case, the value of parameter Н 1 for subject 9 is not much higher than the value of Н 1 = 0.5, which is characteristic of Fick's diffusion. Consequently, this implies that even in a very complex process such as the cerebral neuroactivity we can separate out the components that can be described by physical and mathematical models. The construction of this chaotic component became possible only after the resonant components that continuously transform in vivo during the evolution and therefore defy mathematical simulation were "removed" from the analyzed functions. At the same time, the analysis of the patient signal showed that the function typical for anomalous diffusion cannot be used to describe the behavior of this chaotic component (Fig. 6d), which can be attributed mainly to the presence of a large number of high-frequency components in the original signal. Actually, anomalous diffusion has nothing to do with the dynamics of the patient signal, and the value of coefficient D given in the caption to Fig. 4d, which is three orders of magnitude higher than the corresponding values for the signals of subjects 7 and 9, was used only for the formal comparison of the values of D. It is also necessary to note the anomalously high value of parameter Н 1 = 7.1, which characterizes the nearly complete absence of correlations in the sequence of jumps in the signal (Poissonian process) recorded by sensor 10 in the patient as the response to the RB stimulus.
The data given in Figs The same behavioral pattern is observed for the parameters of the MEG response to the RG stimulus when we analyze the patient responses. The action of the RG stimulus leads to some decrease in the height of the 100 Hz high-frequency peak relative to the 50 Hz peak in the autocorrelator spectrum (Fig. 9b), as compared to the action of the RB stimulus, for which the 100 Hz peak is somewhat higher than the 50 Hz peak (Fig. 5b). As in the RB case, the RGresponse chaotic component Φ (2) (τ) cannot be described by the Levy (anomalous diffusion) approximation (Fig. 9d). It should be noted that the value of component (2) ( ) r Φ τ in the patient MEG signals recorded by sensor 10 is positive throughout the interval τ for both RG and RB stimuli (Fig. 9c).  (Figs. 10b, 11b, 13b, 14b) whereas the collection of these frequencies for the patient is complemented by a rather weak peak at a frequency of 50 Hz (Figs. 12b, 15b). The curves (2) ( ) r Φ τ plotted for sensor 43, which determine the total contribution of specific frequencies and interferential resonant-jump components into the structural functions Φ (2) (τ), have large intervals on the abscissa axis τ in which the values of (2) ( ) r Φ τ are negative not only for the signals recorded in subjects 7 and 9 in response to the RB and RG stimuli (Figs. 10c, 11c, 13c, 14c) but also for the patient signal recorded as the response to the RG stimulus (Fig. 15c). In all these cases, the chaotic component of the structural function Φ (2) (τ) is quite adequately described by the approximation typical for anomalous diffusion (Figs.   10c, 11c, 13c, 14c, 15c). It should be noted that even for the RB stimulus, the curve (2) ( ) r Φ τ plotted for patient's sensor 43 is negative in a small segment of the interval of τ (Fig. 12c).
Although the corresponding chaotic component for the patient cannot be represented in the form typical for anomalous diffusion (Fig. 12c) MEG responses of subject 9 and the patient show a higher internal correlation in the "chaos" of signals recorded as the response to the RB stimulus. In this case, however, the entire set of reported data implies a relatively low diagnostic effectiveness of analyzing the signals produced by sensor 43 for PSE, which is in agreement with the previous general findings [18,19].
Consequently, the FNS analysis of MEG signals recorded as the response to a color stimulus can be used to reliably identify the differences between the parameters of these signals for a patient and healthy subjects and quantify the symptoms of PSE, such as the anomalous group activity of cerebral cortex neurons, best identifiable by the analysis of the signals of sensor 10, which is accompanied by various clinical and paraclinical symptoms. These results are in agreement with the general findings of previous papers [18,19]. For instance, the neuromagnetic responses of the cerebral cortex of a PSE patient show the obvious manifestation of additional quasi-periodic processes, which are absent in the signal of the magnetoelectrical activity of the brain of healthy people. Also, the above results are of interest to the problem of the regulating role of chaos in the dynamics of native systems, which has been discussed for a long time [2,8].
Here, of importance is the behavior of found (2) we can find the solution to Eq. (21) subject to the above initial and boundary conditions: Now we can obtain the expressions for the average value of random variable V and the average squared deviation of this variable from the average value in time τ: From (24), we can derive the asymptotic expressions: Asymptotic expression (A.5) was found using the formula [28]: where ξ ≡ π 2 Dτ/L 2 << 1, in view of the relation: where Ф(x) is the error integral [28]. In (A.7), the change in the discrete values of n in the summation was formally taken to be ∆n = 1, followed by the transition from summation to integration using the integration variable x = ξ n and the differential dξ = ξ when ∆n << 1.