Quantum Theory of EEG with Application to the Single-Trial ERP Analysis

Dmitriy Melkonian ^1,2,*

1. Kaoskey Pty. Ltd. Sydney, Australia

2. Westmead Clinical School, Faculty of Medicine and Health, University of Sydney, Sydney, Australia

* Correspondence: Dmitriy Melkonian

Academic Editor: Raul Valverde

Special Issue: Quantum Brain Dynamics

Received: July 17, 2020 | Accepted: January 17, 2021 | Published: January 27, 2021

OBM Neurobiology 2021, Volume 5, Issue 1, doi:10.21926/obm.neurobiol.2101084

Recommended citation: Melkonian D. Quantum Theory of EEG with Application to the Single-Trial ERP Analysis. OBM Neurobiology 2021; 5(1): 084; doi:10.21926/obm.neurobiol.2101084.

© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

The electroencephalogram (EEG) is a recording of brain’s electrical activity as measured on the scalp which serves as a common probe of brain’s activity, both in clinical and cognitive research settings. It is generally accepted that the EEG originates within the cerebral cortex and incorporates events acting on multiple levels of organization, from molecular to cellular. Single trial event-related potential (ERP) is a specific EEG segment time-locked to some cognitive event.

The processes involved are complex and their interpretation rests mainly on an empirical understanding based on bringing together pieces of multidisciplinary research into the phenomenological hypotheses [1,2].

The attempts to establish quantitative relationships between the global scale EEG and underlying cellular sources are widely researched problems which have modelled the cortical neurons as elementary electric generators of extracellular potentials. The assumption common to almost all approaches is that the membrane potentials produced in some way by the cortical neurons are the “building blocks” the linear summation of which creates the mass potential. The linear model means that a global scale potential contains the parameters of all participating microscopic scale sources of electricity. This leads to intractably huge number of degrees of freedom and prevents a unique determination of the mass effect.

A radically different approach to the treatments of mass potentials, the category of signals to which the EEG belongs, has been introduced by the quantum theory of mass potentials that relocated the elementary sources of the global scale potentials from the cellular to molecular level [3]. The theory preserves a widely accepted concept that measurable changes in potential differences between various global scale locations are produced by synchronous activation of large populations of closely located excitable cells. However, instead of membrane potentials, the role of elementary sources of electricity is prescribed to charged particles, the ions which cross the cell membranes in both directions.

A crucial step is creating a link between the global scale mass potential and the underlying microscale events using the probabilistic formalism of quantum mechanics. A general solution is supported by a particle model using non-homogenous birth and death processes (BDP) for description of trans-membrane transport of ions.

The purpose of this work is to use the quantum theory and corresponding probabilistic tools to bridge the global scale EEG with the underlying molecular sources and create, on this basis, the EEG and ERP models that combine deterministic and probabilistic notions.

The paper is organized as follows.

We first consider the basic models from which the dynamics of the EEG can be synthesized. On the microscopic scale such model is deduced as a stochastic ensemble of extracellular ions called the STION. The corresponding global scale model is deduced as the half wave function (HWF), the model of EEG half waves and ERP cognitive components. A principal step is the bridging of the STION and HWF, i.e. the establishment of analytical relationships between the micro- and macro -scales.

We next consider the results of the creation of an overall EEG model using the time-frequency analysis and fragmentary decomposition. On this basis the universal elements of the EEG dynamics are deduced.

We next consider the results of the developed methods applications to the analysis of event related cognitive potentials, particularly the P3a and P3b.

Finally, we discuss the general properties of the HWF which allow us to consider this function as a transient deterministic chaos.

2. Materials and Methods

2.1 Basic Models

2.1.1 STION

EEG is a global scale potential which resembles electrical phenomena occurring on a microscopic scale consisting of ensembles of multiple cortical cells immersed in an interstitial fluid. There is general agreement that the principal generators of the EEG are cortical neurons, more particularly pyramidal neurons. Numerous EEG models considered the membrane potentials as the origins of the scalp-recorded potentials. Accordingly, the EEGs have been treated in previous studies as summands of transient membrane potentials described by continuous functions of time. Supported by the methods of classical physics, such functions usually represent deterministic solutions of differential equations designed as the models of different membrane configurations. An emerging research field is devoted to a quantum mechanics-based framework for EEG signal feature extraction and classification [4]. In this context a crucial step of our radical departure from deterministic solutions of classical physics to the probabilistic reasoning of the quantum theories is the relocation of elementary bioelectric sources from the cellular to molecular level. Thus, the role of elementary sources of the EEG is prescribed to the ions, both positively and negatively charged particles, which cross the cell membranes in both directions. The probabilistic nature of these microscopic scale events has been discovered by the patch-clamp technique, which provided a means of measuring ion currents through individual ion channels in a cellular membrane [5]. A fundamental finding was that individual ion channels are essentially stochastic entities that open and close in a random way. Thus, the behavior of ions is governed by the probabilistic rules.

Contribution of a single ion to the changes in potential differences between various global scale locations is vanishingly small. Our treatments propose that only cumulative effects of massive trans-membrane transport of ions caused by a coherent activation of a large number of neurons can produce observable changes of the EEG dynamics. The two requirements must be met for this integration to occur: a) the neurons must be closely located and comprise an ensemble, i.e. to have a functional connectivity; (b) the neurons must be active synchronously.

The interior and the exterior of a neuron are both various saline solutions (water with ions dissolved in it) separated by the membrane. Due to a high resistance, compared with the resistances of saline solutions, the cell membranes are rather good electrical insulators [6]. What we have then is two conductors (the inside and the outside of the cell), separated by an insulator (the membrane).

Serving as a border that separates interior of neurons from the extracellular space, the membranes prevent the ions located inside the neurons from producing measurable changes of electric fields in the extracellular space. In contrast, the cumulative effects of the charges of the ions released from the cells during synchronized activation of neuronal ensembles may underlie the dynamics of the global scale EEGs. Thus, we presume that potential fields in the extracellular space are linked to the extracellular charges, i.e. the extracellular populations of ions. In other words, the global scale potential is associated with the net charge of positive and negative ions released from ensembles of pertinent cells and located in the space surrounding their membranes. Such populations of ions comprising the STION is measured at the time t by the integer-valued time-dependent random variable X(t). Here and throughout the paper random entities are indicated by italic or sloping bold type.

Physically the situation represented is this: The population of particles capable to produce a STION develops as a thin cloud of cations and anions spread over the outer surfaces of the cell membranes of the closely located and functionally linked neurons. During the resting conditions the transmembrane ion transport is balanced. This means that, given a local volume, the numbers of positive and negative particles randomly fluctuate over the mean values. Under these conditions the STION does not produce measurable changes of the macroscale voltages. The situation is changed as the result of induced synchronized activity in an ensemble of closely located and functionally linked neurons. Such activity induces transient changes in the membrane machinery that controls the ion exchange between the intracellular and extracellular compartments. The interference of deterministic and random factors produces transient deterministic chaos the limiting distributions of which can be associated with the EEG component waveforms.

There are several probabilistic methods of the quantum theories an essential goal of which is to describe the behavior of multi-particle systems with many degrees of freedom in terms of global systems with few “macroscopic” degrees of freedom [7].

We refer to the central limit theorem as a rule that defines the limiting behavior of ensembles of random variables [8]. This theorem states that the sum of a large number of random samples from independent sources always tends to a limit in the form of normal distribution.

Let n charged particles composing the STION are described by the random variables $f_1,f_2,...,f_n$ with means $\eta_1,\eta_2,...,\eta_n$ and variances,$\sigma_1^2,\sigma_2^2,...,\sigma_n^2$. With increasing number of particles, i.e. n⇒∞, the limiting behavior of their exponential Fourier transforms $F_1(i\omega),F_2(i\omega),...,F_n(i\omega),...$ converges under quite general conditions to the limit probability distribution in the form of the following function of complex variable [9]

\[ G(i\omega)=\exp\{-[(\sigma\omega)^2 / 2]-i\eta\omega\}\tag{1} \]

where $\sigma^2=\sigma_1^2+\sigma_2^2+...+\sigma_n^2,\eta=\eta_1+\eta_2+⋯+\eta_n$ and ω is an angular frequency.

The real and imaginary parts of G($i\omega$) are given by

\[ G_C(\omega)=Re[G(i\omega)]=\exp[-(\sigma\omega)^2 / 2]\cos(\eta\omega), \]

\[ G_S(\omega)=Im[G(i\omega)]=\exp[-(\sigma\omega)^2 / 2]\sin(\eta\omega).\]

As G(iω) belongs to the frequency domain, it can be regarded in terms of the probabilistic notions of the quantum theories as a characteristic function [10]. Accordingly, the time domain counterpart of the characteristic function can be considered as distribution.

Mathematically, the relationship between the frequency domain characteristic function F(iω) and the corresponding time domain distribution function f(t) is established by the reciprocal Fourier integrals:

\[ f(t)=\frac{1}{2\pi}\int_{-\infty}^\infty F(i\omega)\exp(i\omega t)d\omega, \]

\[ F(i\omega)=\int_{-\infty}^\infty f(t)\exp(-i\omega t)dt. \]

These integrals convert one of these representations into the other and vice versa. If F(iω) in the first of these integrals is expressed over an infinite frequency scale by (1), i.e. F(iω) =G(iω), then f(t) is a normal distribution defined on an infinite time scale. However, we deal with a causal process which starts at the time when the triggering event is applied. Such condition of causality establishes specific relationships between the real and imaginary parts of G(iω). On assumption that f(t)=0 at t<0, both $G_C(\omega)$ and $G_S(\omega)$ appear as the counterparts of one and the same time domain function, and either one alone is sufficient to find the time domain counterpart, i.e. the f(t) at t≥0 [9].

Using the imaginary part, we obtain the time domain solution at t≥0 in the form of the following sine Fourier transform of G_S(ω)

\[ \psi(t)=\frac{2}{\pi}\int_0^\infty \exp\left[-[(\sigma\omega)^2 / 2]\right]\sin(\beta\omega)\sin(\omega t)d\omega. \]

This integral has an analytical solution [11]. For t≥0

\[ \psi(t)=\big(\sigma\sqrt{2\pi}\big)^{-1} / [\psi_P(t)-\psi_S(t)], \tag{2} \]

where

\[ \psi_P(t)=\exp[-(t-\eta)^2 / 2\sigma^2],\tag{3} \]

\[ \psi_S(t)=\exp[-(t+\eta)^2 / 2\sigma^2].\tag{4} \]

Taking $\sigma=1\ and\ \eta=1$,these functions are illustrated in the upper panel of Figure 1 by the solid lines: ψ(t) - black, ψ_P(t) - blue and ψ_S (t) – red. The dotted lines are the Gaussian functions which indicate that ψ_P(t) and ψ_S (t) are the fragments of the shifted normal distributions. With σ=1, the lines in the bottom panel show ψ(t) with the various values of η.

Figure 1 Half-wave function and its components. The black, blue, and red solid lines in the top panel show ψ(t), ψ_P(t), and ψ_S(t) with σ=η=1 parameters. The dotted lines are Gaussian functions which indicate that ψ_P(t) and ψ_S(t) are the fragments of the two shifted curves of normal distributions. The bottom panel illustrates HWFs with different values of η.

The equation (2) is consistent with the wave function in a general form of the d'Alembert's solution to the wave equation [12]. In this context ψ(t) is the sum of a right traveling wave ψ_P (t) and a left traveling wave ψ_S(t). A crucial difference is that the d'Alembert's wave function is defined on an infinite time scale while ψ(t) is zero at t<0. In the context of this specific feature we call ψ(t) a half-wave function (HWF).

The frequency domain counterpart of HWF is defined by (1). We repeat this equation here in terms of the amplitude and phase frequency characteristics:

\[ G(i\omega)=A(\omega)\cdot exp[-i\varphi(\omega)],\tag{5} \]

where

\[ A(\omega)=exp[-(\sigma\omega)^2 ⁄ 2]\tag{6} \]

Is the amplitude spectrum and

\[ \varphi(\omega)=\eta\omega\tag{7} \]

is the phase function.

2.1.2 HWF

To identify the microscale origins of the HWF, we now turn from the role of the membrane as an electrical insulator and consider it as a system that controls the transport of the ions from intracellular to extracellular space and vice versa. We wish to describe these transport processes in terms of the birth and death process (BDP). Physically, a particle moving out of a cell would constitute a ‘death’ for the inside of the cell and a ‘birth’ for the extracellular space, i.e. the X(t) ensemble. To describe the temporal evolution of X(t) on the probabilistic basis we refer to the Markov processes.

The main assumption of the Markov process is that during a sufficiently small element of time, Δ, the probability of the change of X(t) by more than one particle is negligibly small [13]:

\[ \mathrm{P}[X(t+\Delta)=X(t)+k]=o(\Delta) \text { if }|k|>1,\tag{8} \]

where P denotes probability and k is an integer.

Remarkable property of the Markov process seen from this equation is that the future state of the population, i.e. X(t+Δ), depends only on its present state X(t), and not on how that state was reached. Thus we deal with a Markov chain, a stochastic model describing a sequence of possible events in which the probability of each event depends only on the state attained in the previous event.

The particle system can change its state only through transition to the nearest neighbours. An increase of the population size by a unit represents birth, X(t+Δ)=X(t)+1, whereas a decrease by a unit represents death, X(t+Δ)=X(t)-1.

Wide classes of BDPs deal with the constant transition probabilities. However, the EEG is a nonstationary signal which indicates the changing structure of the underlying particle ensembles. As the most adequate probabilistic tool we use the BDP of a “transient” type, i.e. a nonhomogeneous BDP in which the birth and death rates may be any specified functions of the time t [14]. The probabilities of population changes for nonhomogeneous BDP are:

\[ P[X(t+\Delta)=X(t)+1]=X(t)\cdot \Delta\cdot \lambda(t)+o(\Delta),\text{(birth)}\tag{9} \]

\[ P[X(t+\Delta)=X(t)]=1-X(t)\cdot \Delta\cdot [\lambda(t)+\mu(t)]+o(\Delta),\text{(no change)}\tag{10} \]

\[ P[X(t+\Delta)=X(t)-1]=X(t)\cdot \Delta\cdot \mu(t)+o(\Delta),\text{(death)}\tag{11} \]

where μ(t) and λ(t) are the birth and death rates, respectively.

Let us choose Δ in compliance with (8) and denote by $x_i=X(t_i)$ the state of the STION at the time t_i. Using these terms we can represent the time evolution of X(t) as the succession of discrete states x_i The permitted states of the particle population at the time $t_{i+1}$ are: x_i+1=x_i+1 (birth), x_i+1=x_i (unchanged size), or x_i+1=x_i-1 (death). The probabilities of the corresponding inter-state transitions are:

\[ P[x_{i+1}=x_i+1]=\hat{P}(i)+o(\Delta),\ \text{(birth)} \]

\[ P[x_{i+1}=x_i-1]=\check{P}(i)+o(\Delta),\ \text{(death)} \]

where $\hat{p}\left(i\right)$ and $\check{p}\left(i\right)$ denote the time dependent transition probabilities for the birth and death, respectively. From (10) and (11) we get:

\[ \hat{p}\left(i\right)=x_{i\ }\cdot\Delta\cdot\lambda\left(t_i\right),\ \ \ \ \ \ \text{(birth)}\tag{12} \]

\[ \check{p}\left(i\right)=x_{i\ }\cdot\Delta\cdot\mu\left(t_i\right),\ \ \ \ \ \text{(death)}\tag{13} \]

We presume that the transient behavior of the STION started from the time instant when some triggering event moves the particle population from the resting to the transient conditions.

Taking t=0 as the beginning of the transient, the expected trajectory of the mean population size e(t) from 0 to t is given by [14]

\[ e\left(t\right)=E\left\{{X}\left(t\right)\right\}=\exp\left[-\mathrm{\rho}\left(t\right)\right], \]

where E {} denotes expected value and

\[ \mathrm{\rho}\left(t\right)=\int_{0}^{t}{\left[\mathrm{\mu}\left(\xi\right)-\mathrm{\lambda} \left(\xi\right)\right]d\xi.} \]

Thus,

\[ e\left(t\right)=\mathrm{exp}\left\{-\int_{0}^{t}[\mu(\xi)-\mathrm{\lambda}(\xi)]{d}\xi\right\}\tag{14} \]

According to the theory of the field potentials within the frequency range of physiological interest the extracellular space may be regarded with reasonable accuracy as a resistive medium [6]. This means that the tissue that lies between the source and the scalp acts as a volume-conductor. Thus, we may assume that temporal evolution of the global scale transient potential v(t) is proportional to the expected population size, i.e. u(t)=k∙e (t), where k is the weighting coefficient.

2.1.3 Bridging Microscopic and Macroscopic Scale Activities

We wish to use the integral from (14) to bridge the global scale EEG expressed by e(t) with the λ(t) and μ(t) variables that govern the microscopic scale processes. We consider the appearance of the two terms in (2) as the indication of the two identifiable particle sub-populations. We define the primary particle population associated with ψ_P(t) and the secondary particle population associated with ψ_S (t).

We first consider the primary particle population, using as a model the BDP with the birth and death rates $\lambda_P\left(t\right)$ and $\mu_P\left(t\right)$ , respectively. Inserting $e\left(t\right)=\psi_P(t)$ in (14) we obtain

\[ \mathrm{exp}\left\{-\int_{0}^{t}\left[\mathrm{\mu}_P\left(\xi\right)-\mathrm{\lambda}_P\left(\xi\right)\right]d\xi\right\}=\mathrm{exp}[-(t-\eta)^2/2\sigma^2]. \]

Consequently,

\[ \int_{0}^{t}{\left[\mathrm{\mu}_P\left(\xi\right)-\mathrm{\lambda}_P\left(\xi\right)\right]d\xi={\left(t-\eta\right)^2}/{2\sigma^2.}} \]

The solution of this equation is a simple matter that provides

\[ \mathrm{\lambda}_P\left(t\right)\mathrm{=}\mathrm{\lambda}_P=\eta/\sigma^2\ and\ \mu_P(t)=t/\sigma^2. \]

Turning to the secondary particle population, we express ψ_S(t) in terms of the birth and death rates denoted by $\lambda_S\left(t\right)$ and $\mu_S\left(t\right)$.

The replacement of e(t) in (14) by ψ_S(t) gives

\[ \mathrm{exp}\left\{-\int_{0}^{t}\left[\mathrm{\mu}_\mathrm{S}\left(\xi\right)-\mathrm{\lambda}_\mathrm{S}\left(\xi\right)\right]d\xi\right\}=\mathrm{exp}[-(t+\eta)^2/{2\sigma^2}] \]

Consequently,

\[ \int_{0}^{t}{\left[\mathrm{\mu}_\mathrm{S}\left(\xi\right)-\mathrm{\lambda}_\mathrm{S}\left(\xi\right)\right]d\xi={\left(t+\eta\right)^2}/{2\sigma^2}}. \]

It follows from solution of this equation that

\[ \mathrm{\lambda}_\mathrm{S}\left(t\right)=0\ \ \mathrm{and} \ \mathrm{\mu}_\mathrm{S}\left(t\right)={\left(t+\eta\right)}/{\sigma^2=\lambda_P+\mathrm{\mu}_\mathrm{P}\left(t\right).} \]

As in [3], we summarize these relationships in the form of the following rules governing the birth and death rates for the primary and secondary particle populations.

Rule 1. After activation at t=t₀ by the triggering event the transient behavior of the primary particle population develops as a non-homogenous BDP with the constant rate of birth

\[ \lambda_P={\eta}/{\sigma^2}.\tag{15} \]

and the time dependent rate of death

\[ \mathrm{\mu}_\mathrm{P}\left(t\right)={\left(t-t_0\right)}/{\sigma^2}.\tag{16} \]

Rule 2. After activation at t=t₀ by the triggering event the transient behavior of the secondary particle population develops as a non-homogenous death process with the time-dependent rate of death

\[ \mathrm{\mu}_\mathrm{S}\left(t\right)={\left(t-t_0+\eta\right)}/{\sigma^2}=\mathrm{\mu}_\mathrm{P}\left(t\right)+\lambda_P.\tag{17} \]

Turning to the resting conditions, it is essential that the ion transport is balanced for both the cations and anions. This means that in the resting states the sizes of the primary and secondary particle populations fluctuate over the mean values. We consider development of these events as simple BDPs with the constant rates of the birth and death. To define parameters of these processes we note that (15)-(17) contain the time dependent term μ_P(t) and the time independent term λ_P. On the physical basis we associate the time dependent μ_P(t) with the gated ion channels and the constant λ_P with the resting ion channels. In this context λ_P participates in both the resting and transient conditions. Thus, we consider λ_P as a universal estimate for the birth and death rates during the resting conditions. Accordingly, the ion transport, balanced for both particle populations, is supported by the condition:

\[ \lambda_P^r=\mu_P^r=\lambda_S^r=\mu_S^r={\beta}/{\sigma^2},\tag{18} \]

where $\lambda_P^r,\ \mu_P^r,\ \lambda_S^r$ and $\mu_S^r$ denote the resting state rates of the birth and death for the primary and secondary particle populations.

It is important to note that resting state conditions are not recognizable from the global scale. The parameters included to (18) are deductions from the rules governing the transient regimes. Thus we may expect existence of additional resting state parameters which do not affect the transient components.

2.1.4 Algorithm of Numerical Simulations

In order to construct numerical simulations using the above formulated rules we consider the numerical counterpart of (2) in the form: X_N(t)= X_P(t)-X_S(t), where X_P(t) and X_S(t) are the numbers of the particles in the primary and secondary populations, respectively. Accordingly, X_N(t) is the number of particles producing the net charge. Often it is not necessary to link the number of particles with concrete population. In such general cases the subscript is omitted for brevity. Thus, x_i denotes the number of particles at time instant t_i, i.e. x_i=X (t_i).

We consider the time instants t_i=iΔ (i=0, 1, …), where t₀ is the time at which the EEG transient starts. With the reference to (12) and (13) we deduce the following formulas of transition probabilities under the transient conditions.

Primary particle population:

\[ {\hat{p}}_p\left(i\right)={x_i\cdot\Delta\cdot\eta}/{\sigma^2},\ \ \ \ \ \ \ \ \text{(birth)}\tag{19} \]

\[ {\check{p}}_P\left(i\right)={x_i\cdot\Delta^2\cdot i}/{\sigma^2}.\ \ \ \ \ \ \ \text{(death)}\tag{20} \]

Secondary particle population:

\[ {\check{p}}_S\left(i\right)={x_i\cdot\Delta\cdot\left(i\Delta+\eta\right)}/{\sigma^2}.\ \ \ \ \ \text{(death)}\tag{21} \]

According to (15) the birth and death processes in the primary and secondary particle populations are governed under the resting conditions by the birth and death rates equal to η/σ². The corresponding transition probability is

\[ p_r\left(i\right)={x_i\cdot\Delta\cdot\eta}/{\sigma^2.\ }\tag{22} \]

Now we have the explicit equations (19)-(21) for the transition probabilities that provide means for numerical reconstructions of the time evolution of the particle populations during the transient and resting conditions. A general procedure is organized as the succession of standard steps dealing in sequential order with the time intervals [t_i, t_i+1] for i taking values from -M to N-1. The time interval from t_-M to t₀ corresponds to the resting conditions. At t=t₀ the resting state is switched to the transient conditions.

For the primary particle population the procedure for each step is as follows:

(1) Set x_i to be the initial state of the particle population. For the first interval beginning at t_-M define the initial state by an arbitrary integer N₀, i.e. $x_{-M}=N_{\mathrm{0.}}$.

(2) Estimate the transient probabilities from (19) and (20) for the transient conditions or use (22) for the resting conditions.

(3) Pick out random real numbers $\hat{r}$ and $\breve{r}$. using a random number generator to produce the real numbers in the range from 0 to 1.

(4) Estimate the size of the particle population at the end of the interval $x_{i+1}=x_i+b-d$,

where b and d are the binary numbers defined as follows.

Resting conditions:

b=1 if r_b<p_r(i) and is zero otherwise,

d=1 if r_d<p_r(i) and is zero otherwise.

Transient conditions:

b=1 if $\mathbf{r}_b<{\hat{p}}_P\left(i\right)$ and is zero otherwise,

d=1 if $\mathbf{r}_d<{\check{p}}_p\left(i\right)$ and is zero otherwise.

For the secondary particle population the procedures are similar, with the exception that components of the birth process under the transient conditions are excluded from considerations.

2.2 Fragmentary Decomposition

There is a general agreement that EEG is a complex signal composed from multiple components produced by activation of various cortical processors. As objective markers of the components we consider the positive and negative peaking waveforms in the EEG time course. Such conceptual approach is in harmony with a widely accepted notion that various positive and negative deflections of electrophysiological signals are functional entities that reflect the summed activity of the underlying cellular populations.

To identify a component we need to estimate segments over which the EEG deflections develop. The methodology adopted here is that of the method of the fragmentary decomposition suggested in earlier works [15,16,17].

Considering the EEG as a time function v(t), we deal with a series of samples $v_m=v\left(t_m\right)$ at regular, discrete time points $t_m=m\Delta$ where Δ is the sampling interval.

The segmentation points are defined as zero-crossings and points of global and local minima of |v(t)|. More particularly, if,

\[ \left(\mathrm{v}_{m-1\ }\le0\ \ \mathrm{AND} \ \mathrm{v}_{m+1}>0\ \right)\ \ \mathrm{OR}\ \mathrm{v}_{m-1\ }\geq0\ \ \mathrm{AND\ \ }\mathrm{v}_{m+1}<0\ , \]

then τ_k =t_m is qualified as zero-crossing (segmentation point), where k is the number of the segmentation point.

The point τ_k =t_m (k is the number of the segmentation point) specifies a global or local minimum of |v(t)| if |v_m-1|≥|v_m|≤|v_m+1|. By ordering the segmentation points, both zero-crossings and minimums, as consecutive time points the sequence of the segmentation points τ₀, …, τ_k, …, τ_K is formed.

The EEG segment between sequential segmentation points is called an empirical half wave (EHW). Given a segment of the length T_k = τ_k-τ_k-1 between the points τ_k-1 and τ_k (i=1, …, K), the EHF is defined as

\[ \mathrm{w}_k\left(t\right)=\left\{\begin{matrix}\mathrm{v}\left(t+\tau_{k-1}\right)\ if\ \ 0\le t\le T_{i\ },\\0\ \ \mathrm{otherwise.\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\\\end{matrix}\right. \]

In the interval from 0 to T_i this function reproduces the EEG fragment between the segmentation points τ_k-1 and τ_k. Since the following manipulations with w_k(t) are universal, we omit the “k” subscript and consider w(t) as EHW defined on the interval [0, T].

The frequency domain counterpart of the w(t) is defined by the exponential finite Fourier transform

\[ W\left(i\omega\right)=\mathrm{W}_C\left(\omega\right)-{\mathrm{jW}}_S\left(\omega\right)=\int_{0}^{T}w\left(t\right)\mathrm{exp} \left(-i\omega t\right)dt,\tag{23} \]

where ω=2πf, f is frequency and j=√-1.

Since w(t) is an empirical entity, the calculations of such integrals are performed numerically. The most readily available computational technique of standard Fourier analysis is an efficient algorithm of the fast Fourier transform (FFT) [18]. The critical aspect of this technique is that the FFT is supported by the Fourier series model of the data. This distinction with the Fourier integral poses a number of limitations, the most important of which is that the FFT is applied to 2ⁿ equidistant points of the function to be transformed. The integer “n” and the sampling interval define the set of equidistant points at which the transformation results are computed. A strict grid of computational points creates a number of methodological difficulties. For example, it is incompatible with the widely accepted logarithmic frequency scales for computations and display of the frequency characteristics. The distinction of the FFT with the Fourier integral is also the source of specific computational errors.

The alternative numerical method we use in this study is the SBF algorithm [19], a novel version of the Filon-type methods that provide maximum precision in the estimation of trigonometric integrals. The SBF algorithm deals with the continuous Fourier spectrum instead of a discrete spectrum defined by the FFT. This allows us to use a logarithmic frequency scale for calculations and display of frequency characteristics. Another point is that dealing with the EEG segments of different length the SBF algorithm does not demand windowing and zero-padding of the time series in question. This is especially important for our algorithms dealing with EHWs the lengths of which are variable. Due to these advantages the algorithm is applicable for the EEG short-time spectral analysis using segments of arbitrary length. The explanations of the SBF algorithm may be found in the appendix.

The SBF algorithm is addressed to the estimation of the real and imaginary parts of the complex spectrum (23):

the finite cosine Fourier transformation

\[ W_C\left(\omega\right)=\int_{0}^{T}{w\left(t\right)cos\left(\omega t\right)dt,} \]

and the finite sine Fourier transformation

\[ W_S\left(\omega\right)=\int_{0}^{T}{w\left(t\right)sin\left(\omega t\right)dt,} \]

where ω=2πf is the angular velocity, and f is the frequency.

The corresponding amplitude spectrum and phase function are defined as:

\[ \begin{gather} W\left(\omega\right)=\left|W\left(i\omega\right)\right|=\sqrt{W_C^2\left(\omega\right)+W_S^2\left(\omega\right)},\\ \xi\left(\omega\right)=\mathrm{arctg}\left[{W_S\left(\omega\right)}/{W_C\left(\omega\right)}\right]. \end{gather} \]

These are continuous functions of the angular frequency, the values of which can be calculated for arbitrary sets of points. We employ a logarithmic frequency scale for the calculations of the amplitude spectra and a natural frequency scale for calculations of the phase functions.

In the case of the logarithmic scale, the frequency characteristics are calculated for the angular frequencies ω_i=ω₀C^i-1 (i=1, …, N), where ω₀ is initial angular frequency and C>1 is the parameter that defines the sampling rate. For selection of this parameter it is convenient to use the formula C=exp(ln10/N_D), where N_D is the number of samples per decade. In the case of a natural frequency scale, ω_i=ω₀+iΔω (i=0, …, N), where Δω is the discretization step.

We normalize the amplitude spectrum and the frequency in order to express these entities in dimensionless units.

The normalized amplitude spectrum is defined as

\[ W^\ast\left(\omega\right)={W\left(\omega\right)}/{W\left(\omega_0\right)}, \]

where W(ω) is computed amplitude spectrum and ω₀ is sufficiently small value of the angular frequency selected to satisfy the condition: $W\left(\omega_0\ \right)\approx W\left(0\right)$.

Typical normalized amplitude spectrum is illustrated in the middle panel of Figure 2. It corresponds to the half-wave indicated by the black dotted lines in the upper panel.

Figure 2 The black line in the top panel illustrates typical EEG. The middle and bottom panels show the amplitude apectrum and phase function, respectively.

We may regard W^*(ω) as the frequency response of a low pass filter, conventional parameter of which is the cut-off frequency F_C. At this frequency the attenuation of the amplitude spectrum drops by 3dB, i.e.

\[ W^\ast\left(\omega_C\right)={1}/{\sqrt{2\ }} \]

where ω_C=2πF_C.

For normalization of the frequency scale we use F_C as a basis unit and define the normalized dimensionless relative frequency γ=ω/2πF_C=ω/ω_C.

Appropriately scaled in the magnitude and frequency, the standard empirical amplitude spectrum is defined as $Z\left(\gamma\right)=W^\ast\left(\omega_C\gamma\right)$

Associating empirical W(ω) with the theoretically deduced analytical (6) we consider an expected form of Z(γ) as the Gaussian spectrum

\[ G(\gamma)=\exp(-\gamma^2)\tag{24} \]

Note that $Z\left(\gamma\right)=G\left(\gamma\right)\ at\ \gamma=1$ the relative frequency which corresponds to f=F_C.

The combining Z(γ) and G(γ) at γ=1 provides an excellent fit of analytical G(γ) to empirical Z(γ). We use this quite universal feature to create a procedure for estimation the σ using G(γ) as a frequency domain template for comparison with the amplitude spectra of EHWs.

Template matching is performed by moving the Gaussian amplitude template throughout the abscissa scale. This is a step by step procedure during which the template samples remain unchanged while their locations on the abscissa scale are jointly shifted by the multiplication of the frequencies by the C constant.

The first shift starts from the sequence of points $\nu_m=\nu_AC^{m+1}\ \left(m=0,\ 1,\ \ldots,\ M\right)$. The sequence of the frequency points after the kth shift is $\nu_m=\nu_AC^{k+1}\ \left(m=0,\ 1,\ \ldots,\ M\right)$.

For evaluation of the accuracy of the fits the reference is made to the following discrete values of the relative frequency

\[ \gamma\left[i\right]={\omega_0C^{i-1}}/{2\pi\ \ \left(i=1,2,\ldots\right)} \]

As a measure of accuracy we use the MSE[m,n] which represents the mean square error calculated at the range of frequencies γ[i] from i=m to i=n.

Numerous trials with various EHFs revealed that typically the fits virtually coincide with the amplitude spectra in the range of standard frequencies from 0 to 1. At γ>1 the errors increase with increase of the frequency. Taking into account these peculiarities, the tests of accuracy are organized as a two-step procedure supported by the empirically established error thresholds T1=0.0001 and T2=0.002.

Step 1. The range of the relative frequencies from γ [0] to γ[J]=1 is selected. It corresponds to the frequency band from F₀=ω₀/2π to F_C. If the MSE [0, J] < T1 the fit is accepted and followed by the step 2.

Step 2. Starting from i=1 the MSE [J+i-2, J+i+3] is calculated for increasing numbers of i until MSE [J+i-2, J+i+3]>T2. This is a six point window which proceeds from F_C to the higher frequencies. The point m=i at which the procedure stops defines the boundary frequency F_B=γ[m]F_C. Thus, an acceptable fit extends over the frequency range from F₀ to F_B.

For comparison of these spectra the reference is made to the discrete values of the relative frequency γ denoted as γ[i]=γ₀C^i-1 (i=1, 2, …) where γ₀=ω₀/ω_C. These are evenly distributed points on the logarithmic abscissa scale.

The larger F_B is in comparison with F₀, the more accurate is the Gaussian model of the amplitude spectrum. To employ this attribute for assessing the goodness of fit, we use the dimensionless extension ratio ε=F_B/F_C. The fit in the middle panel of Figure 2 gives a visual idea of how the F_C and F_B are estimated. Taking F_C from the best fit, the $\mathrm{\sigma=}{\sqrt{\mathrm{ln2}}}/{\mathrm{2\pi}\mathrm{F}_\mathrm{c}}$.

The calculations of σ are followed by the estimation of the phase functions. We found that initial part of the phase function φ(ω) shows consistency with a simple linear model

\[ \phi\left(\omega\right)=\xi\omega\tag{25} \]

where ξ is a parameter.

As numerous calculated linear fits indicate, the deviation from linearity can be neglected over the frequency range from 0 to 1.4∙F_C. A typical result is illustrated in the bottom panel of Figure 2. Thus, the estimation of η is reduced to the calculation of the linear regression lines using Φ(f) samples from f=f₀ to f=1.4∙F_C. The slope of the regression line ξ serves as the estimate of the parameter η.

2.3 Kolmogorov Smirnov Test

The above procedures of the parameter estimation define the frequency range of the best fit of the theoretical amplitude spectrum to the empirical amplitude spectrum of EHW from F₀ to F_B. An opportunity to compare different fits is provided by the dimensionless extension ratio $\varepsilon={F_B}/{F_C}$ Given the samples of ε in the form of two different ensembles, $\varepsilon_1=\left\{\varepsilon_{1,\ }^1,\ ..,\varepsilon_n^1,..,\varepsilon_N^1\right\}$ and $\varepsilon_2=\left\{\varepsilon_{1,\ }^2,\ ..,\varepsilon_k^2,..,\varepsilon_K^2\right\}$, we use the Kolmogorov-Smirnov one and two sample tests in order to decide whether $\varepsilon_{1}$ and $\varepsilon_2$ are produced by the same or different distributions. Each of the data sets $\varepsilon_{1}$ and $\varepsilon_2$ is converted to the cumulative frequency distribution. The test is based on the evaluation of the maximum vertical deviation D between the cumulative frequency distributions. The null hypothesis that the two distributions are the same is rejected if the value of D exceeds the critical value defined by the tables of D statistics.

2.4 Materials

The EEG and ERP data have been collected and described in the study of single-trial auditory event related potentials [16]. A standard auditory “oddball” paradigm was employed. The EEGs were recorded from 40 healthy subjects (20 males and 20 females; age 20-50 years). Stereo headphones conveyed in a random order 1500-Hz task-relevant tones and 1000-Hz task-irrelevant tones. EEGs were recorded from 19 electrode sites according to the 10-20 international system and off-line artefact corrected. The EEG segments for analysis have been extracted from a digitally stored data sets from both the pre- and post- stimuli intervals.

An initial signal is electrical potential recorded from certain scalp location. Such continuous time signal represents a mixture of the EEG, i.e. the potentials produced by electrical activity of the brain, with the noise and various artefacts. Numerous algorithms have been developed to exclude the components unrelated to the brain activity and to deal with the EEG in the form of a time series. The technique we employ has been described earlier elsewhere [16,20].

3. Results

3.1 Universal Elements of the EEG Dynamics

Typical feature of EEG is the tendency toward repetition of some characteristic activity patterns over time. Accordingly, most of the methods of EEG quantitative analyses are supported by the concept that EEG is a composite of several oscillatory processes the most important of which are the delta, theta, alpha and beta band rhythms. The strongest EEG waves are occipital alpha rhythms recorded in the frequency range of 8–12 Hz during the periods of eyes closed. The blue line in the upper panel of Figure 3 exemplifies typical alpha rhythm.

Figure 3 Typical EEG with α rhythm in the top panel and its amplitude spectrum in the bottom panel.

Traditional toolkits of the EEG rhythm exploration are the frequency domain methods. A common methodology called a standard frequency analysis (SFA) [21] consists in the estimation of the EEG power or amplitude spectra using numerical calculations of the Fourier transforms. Usually the FFT is employed. The advantage of the SBF algorithm is that the frequency characteristics are treated as continuous functions of the frequency and can be estimated for any selected EEG segments.

Since the SBF algorithm is applicable to the signals of an arbitrary length, it is universal and can be used as a tool of the SFA. In this role the SBF algorithm has been utilized in the study of patients with the first episode and chronic Schizophrenia [20]. Methodological innovations of the SBF algorithm supported the overarching finding that increased delta band activity remains the essential difference between the age matched groups of subjects with schizophrenia and the normal population, but that this is most pronounced for subjects with chronic schizophrenia.

Particular advantage of the SBF algorithm is the employment of logarithmic frequency scales which are suitable because frequency domain parameters tend to vary linearly with the logarithm of the frequency. This is seen from typical EEG amplitude spectrum depicted in the bottom panel of Figure 3. The major feature is well pronounced large resonant peak at the frequency 9.55 Hz indicated by the arrow 1. This is an unambiguous indication of the alpha rhythm. The three peaks of smaller amplitudes indicated by the arrows 2-4 correspond to the frequencies 2.75, 3.89 and 18.6 Hz, respectively. It is tempting to speculate that such secondary peaks indicate some additional periodic EEG components.

Unfortunately, such expectations are not supported by an adequate theoretical basis. A critical point is that the mathematical backgrounds of the SFA are addressed to the stationary processes, both random and periodic. In the case of a periodic signal the advantage of the SFA is that it decomposes the signal into the sum of harmonically related trigonometric components and establishes the relative intensity of each component.

This elegant mathematics of Fourier methods does not, however, adequately work in the case of the EEG the activity patterns of which are highly irregular and non-stationary. A heuristic, rather than a mathematically rigorous approach, considers EEG segments with strong rhythmicity as quasi-periodic processes. The weakness of this approach is that besides changes in the morphology of peaks and troughs the periodic tendencies underlying the EEG temporal development are combined with different forms of irregularities.

In such cases, similar to the situation illustrated in Figure 3, the major outcome of the SFA is the identification of a dominant resonant peak. As a rule, additional peaks of smaller amplitudes are not uniquely defined and testable.

The power of the short time-frequency analysis applied by the SBF algorithm to each EHW is that it works like a “mathematical microscope” and allows us to retain a full extent of the information contained in the time series of non-stationary EEG. Figure 4 illustrates typical results of the time frequency analysis applied to the EEG signal with 10 peaking waveforms shown in the top panel. Each waveform from 1 to 10 separated by vertical dotted black lines has been transformed to the frequency domain using the SBF algorithm. Normalized by both the amplitude and frequency these empirical amplitude spectra are depicted in the bottom panel. The black line denoted by G is the Gaussian spectrum defined by (24).

Figure 4 EEG recording in the top panel is divided into ten segments indicated by the dotted black lines. The corresponding amplitude spectra normalized by both the magnitude and frequency are shown in the bottom panel.

The key observation is a mutual coincidence of the empirical amplitude spectra in the wide range of the relative frequency γ from 0.1 to 100. The Gaussian function appears as a limiting form of these empirical dependencies. This result can be qualified as a phenomenon known in the quantum theories as universality [22].

Conceptually, the universality means that despite a profound diversity of complex systems observed in nature and technology, their topology may possess universal objects.

The formalism of our theory links EHWs to the class of universal elements the appearance of which is due to the statistical factors governing the collective behavior of molecular events underlying the EEG genesis. Such universal elements are identifiable in the EEG dynamic structure without requiring the knowledge of the details of the microscale system. The problem we now consider is how to create adequate models of EEG waveforms using universal objects as the building blocks.

The major point seen from Figure 4 is a universal profile of the amplitude spectrum of EHW and its agreement with the Gaussian spectrum (24). Thus, we can consider theoretical HWF as a model counterpart of EHW. Consequently, the model of overall EEG with N EHWs can be made up as the following sum of HWFs

\[ h\left(t\right)=\sum_{i=1}^{N}{\kappa_i\cdot\mathrm{\psi}_\mathrm{i}\left(t-\tau_{i-1}\right)},\tag{26} \]

where index “i” labels different HWFs with corresponding σ_i and η_i parameters, κ_i is the weighting coefficient and τ_i is the time instant from which the corresponding HWF starts.

The weighting coefficient is derived from the interpolation conditions separately formulated for each EHW

\[ \kappa_i\cdot\psi_i\left(L_i\right)=w_i\left(L_i\right), \]

where L_i is the peak latency of the ith EHW. Consequently,

\[ \kappa_i={w_i\left(L_i\right)}/{\psi_i\left(L_i\right).} \]

Thus, the peak latencies and amplitudes of the model are equalized to the peak latencies and amplitudes of original EEG.

The results of numerous tests revealed remarkable accuracy of the model (26) at those EEG segments where EHWs durations didn’t significantly exceed the intervals between the corresponding segmentation points. If these conditions are not satisfied, the model is likely to be contaminated by the temporal overlap of EEG components.

To exclude the effects of the component overlap we implement the procedure of the step-by-step resolution of the component overlap. We assume that analytical ψ_m(t) started from t=τ_m-1 captures the dynamics of empirical w_m(t) at the interval from τ_m-1 to τ_m and contributes to the EEG development at t>τ_m. By contrast, at $t>\tau_m$ the impact of the ψ_m-1 (t) is vanishingly small. Thus, we may consider the EEG model in the segment from $\tau_{m\ }$ to $\tau_{m+1}$ as the following extract from the general model (26)

\[ \theta\left(t\right)=\kappa_m\cdot\mathrm{\psi}_m\left(t-\tau_{m-1}\right)+\kappa_{m+1}\cdot\mathrm{\psi}_{m+1}\left(t-\tau_m\right). \]

On assumption that $\theta\left(t\right)\approx v(t)$ in the $\left[\tau_m\ ,\ \tau_{m+1}\right]$ interval, we define on this interval the overlap corrected EHW in the form

\[ {\widetilde{w}}_{m+1}\left(t\right)=v\left(t\right)-\kappa_m\cdot\psi_m\left(t-\tau_{m-1}\right). \]

This formula supports processing of the mth segment taking into account the trajectory of the previous segment.

Schematic illustration of the significance of the overlap correction procedure is provided in Figure 5. The black line shows EEG fragment on the time interval of 800 ms. The green line is a model the two elements of which, HWF1 and HWF2, are shown by the red and blue lines, repsectively. Consider the EEG maximum at the time instant indicated by the black vertical line. At this point it is larger by ΔA the amplitude of the HWF2. This difference is compensated by the overlap of the right flank of the HWF1 with the left flank of the HEF2. It is important to note that the overlap account and correction are only possible if the components are expressed, as in the case of our theory, by sufficiently accurate models.

Figure 5 Illustrates how the temporal overlap of subsequent positive components is eliminated and provides means for creation of adequate models of ERP components.

To utilize the described method of the overlap correction, the identification of GMP parameters is organized as the recursion procedure performed in successive steps from the first segment to the last segments of the EEG in question. At the first step the effect of the previous HWF is neglected. Each following step of the recursion includs renewal of the segmentation point and identification of all parameters of the overlap corrected EHWs.

The high accuracy is a remarkable outcome of the fragmentary decompositions supported by the overlap corrections of the EHWs building the model. Typical result is illustrated in the top panel of Figure 4 where the red line shows the model.

An important aspect of the theory is that the models are derived entirely from EHWs without any adjustments to the various origins of the underlying molecular and cellular systems.

An obvious explanation is that in practically all artefact and noise free cases an empirical $W(\omega)\ and\ \xi(\omega)$ show the features of the normal distribution. Given that the pair $W(\omega) and \varphi(\omega)$ is in every way equivalent to their time domain counterpart w(t), it is clear that the time domain universal objects should also follow to the predictions of the central limit theorem.

Because we deal with empirical entities, an important problem was to determine whether the calculated empirical distribution functions are different or identical in statistical terms. We applied the Kolmogorov-Smirnov tests to examine the null hypothesis that the samples follow the normal distribution. The test has an advantage of making no assumption about the distribution of the samples, i.e. it is non-parametric and distribution free. The parameter ε has been selected as an adequate parameter for these tests because this single parameter is sufficient to evaluate the fitting results.

We demonstrate the typical results using collections of ε which were identified using EEG recordings from Fz, Cz and Pz cortical sites. For each location we have collected values of ε from 20 different identified EHWs. The EHW was accepted as eligible for testing the D-statistics if the number of the EHW’s samples was ≥8.

The means, standard deviations (SD) and D (K-S statistic) were as follows.

Fz: mean =1.756 (SD=0.205), D=0.106.

Cz: mean=1.712 (SD=0.152), D=0.076.

Pz: mean=1.721 (SD=0.182), D=0.079.

Normalized by the sample size, empirical cumulative distribution functions evaluated from each of these sample sets are illustrated in Figure 6. The blue lines in upper panels show the pertinent EEGs from Fz, Cz and Pz cortical sites. The red lines are the models calculated using the fragmentary decomposition with overlap correction. The lower panel shows the corresponding cumulative distribution functions normalized by the sample size. The line denoted by CND is the cumulative normal distribution.

Figure 6 Upper panels show EEGs from Fz, Cz and Pz recording sites and their remarkably accurate models calculated by the method of fragmentary decomposition with facilities of the elimination of a component temporal overlap. In the bottom panel the cumulative distribution functions of these processes, i.e. empirical functions, are compared with the curve of cumulative normal distribution to accomplish one-sample Kolmogorov-Smirnoff test.

The greatest discrepancy between the CND and empirical cumulative distribution, called the D-statistic, serves as a criterion to reject the null hypothesis. Given that all empirical cumulative distributions have been supported by equal numbers of ε parameter (40 values of ε employed in the tests), the null hypothesis is rejected if D≥0.189 (5%).

The presented estimates are well below this value and do not provide any reasons to reject the null hypothesis. It is important to note a highly stereotypical character of the test results. The outcomes of multiple tests didn’t reveal any dependency between the empirical distributions and the cortical sites from which the EEGs were recorded.

3.2 Event Related Cognitive Potentials

Exquisite sensitivity of EEG to the changes in mental activity has been documented in numerous studies. The fundamental criterion used to link the EEG dynamics with cognitive events is the concomitant variation of neurophysiological and psychological processes. Consequently, the key to the understanding of the information processing context of EEG signals is provided by detection of the changes in the ongoing EEG activities that are time locked (or at least closely related in timing) to particular cognitive event. Specific changes of the EEG dynamics elicited by cognitive variables are known as event related potentials (ERPs). Typical activation task used to elicit ERP in human subjects is the ‘oddball’ paradigm in which the task-related ‘oddball’ (or target) stimuli occur randomly among other ‘background’ (or non-target) stimuli. See [16] for details.

Because ERPs are embedded in ongoing EEG activity, signal extraction methods are absolutely essential in ERP analyses. The most common method is the ensemble averaging based on the model of a single trial ERP as a deterministic transient potential strictly time-locked to the pertinent stimulus. This is well known oversimplification which discards the probabilistic nature of ERPs, and creates ambiguity with regard to the analysis and interpretation of the endogenous potentials. Fundamental aspect of this problem is that non-stationary single trial ERPs are expected to have changing patterns from trial to trial. The generality of such understanding is particularly seen from the following psychophysiological definition of endogenous ERP components: “The components must be nonobligatory responses to stimuli. The same physical stimulus, presented to the same subjects, sometimes will and sometimes will not elicit the component” [23].

As a result, conventional ensemble average would not necessarily correspond to any of the individual single-trial responses, which means that the averaged waveforms may obscure an important information contained in the specific changes of the EEG signal.

3.2.1 Single Trial ERPs

To improve the situation, a lot of effort was invested to invent the methods of estimating the ERPs directly from single trials. Mathematical difficulties arise from the fact that ERPs belong to the category of non-stationary processes. This class of signals is not supported by a general theory, and relies on heuristic, rather than mathematically rigorous methods. The most common approach is based on the assumption that EEG may be decomposed into a number of characteristic waveforms linked to the mechanisms of the underlying systems. In this context the analysis of a non-stationary process is synonymous with creation of the corresponding models.

To the best of our knowledge, the presented quantum theory is the first mathematical and computational tool which deduces empirically testable adequate models of ERP components directly from the EEG recordings. This development is supported by the finding of universal objects in the EEG dynamics and creation on this basis the EEG model composed from analytical HWFs. Being linked to the pertinent EHWs, these entities are considered as potential models of cognitive components.

On the grounds of significant experimental and clinical evidence, the EEG peaking waveforms that appear at specific time instants after the application of a cognitive stimulus are regarded as ERPs, i.e. the functional entities associated with various cortical neuronal networks. In the auditory oddball paradigm, such peak potentials are known as the late latency ERPs denoted as P50, N100, P200, N200, p300a (P3a) and P300b (P3b) components. The number indicates an approximate time (in milliseconds) of a negative (N) or positive (P) peak appearance after the stimulus application.

Using the amplitude and latency as major parameters, we support the choice of “true” cognitive components by the following conditions:

Condition 1: L₁≤L<L₂ (latency window),

Condition 2: A₁≤A<A₂ (amplitude window),

where A and L are the amplitude and latency parameters of the pertinent HWF.

The following L₁ and L₂ (ms) parameters of the latency windows have been selected: P50 -20 versus 60, N1(00) - 80 versus 120, P2(00) - 160 versus 220, N2(00) - 180 versus 235, P3a - 240 versus 299, P3b - 300 versus 360. The A₁ and A₂ (μV) parameters of the amplitude windows are -45 versus -2 for negative and 2 versus 45 for positive components, respectively.

The top panel of Figure 7 shows the EEG which can be considered as a typical example of a single trial ERPs elicited in an auditory oddball paradigm. The dynamics of the EEG signal shown in the upper panel of Figure 7 is affected at t=0 by the target auditory stimulus. This event is followed by specific succession of positive and negative waveform deflections. The EEG segment beginning from this time instant is redrawn with the higher temporal resolution in the bottom panel. It represents a single trial ERP.

Figure 7 Top panel illustrates how ongoing EEG transforms in an auditory oddball paradigm to the succession of cognitive components after target stimulus application at t=0. The bottom panel shows models of these cognitive components and their ensemble which appears as a model of the whole single trial ERP.

On the basis of above conditions, six HWFs shown by the colored lines in the lower panel of Figure 7 are qualified as P50, N100, N200, P300a (p3a) and P300b (P3b) late cognitive components.

3.2.2 Selective Component Averaging

An explicit analysis of the single trial ERP variability necessitates consideration of various ensembles of the identified single trial parameters. The corresponding procedure is called selective component averaging (SCA). We address the SCA to particular component defined by an ensemble of corresponding HWFs from single trials. From these entities the HWFs are selected which satisfy conditions 1 and 2.

The SCA of selected components is defined by the following sum

\[ u^D\left(t\right)=\frac{1}{M_i}\sum_{i=1}^{M_i}{e_i^D\left(t-\tau_{Di}\right)}, \]

where the symbol “D” stands for the name of the component, M_i is the number of selected HWFs and τ_Di is the time instant from which the HWF starts. To create a composite of several ERP components such sums are calculated for each selected component. For example, selective average of the late positive complex (LPC) consisting of the P3a and P3b components is computed as follows,

\[ u^{LPC}\left(t\right)=\frac{1}{M_a}\sum_{i=1}^{M_a}{e_i^{P3a}\left(t-\tau_{P3ai}\right)}+\frac{1}{M_b}\sum_{i=1}^{M_b}{e_i^{P3b}\left(t-\tau_{P3bi}\right)}, \]

where $M_a$ and $M_b$ are the numbers of $p3_a$ and $p3_b$ components delivered to the SCA.

The SCA improves accuracy of the average measures because it selects the trials with identified, i.e. meaningful components, and ignores the trials with missing component.

For account of missing responses we introduce a novel parameter called an elicitation rate (ER). This parameter takes into account an actual number A of the trials in which the component was defined: A=T-M where T is the total number of single trials and M is the number of the trials with missing component. The ER is defined as $P_E={A}/{T}$, i.e. this parameter is the probability of the elicitation of defined component in a single trial.

Typical differences in probabilities with which the endogenous potentials appear in single trials of healthy subjects are illustrated in Figure 8 which refers to the P50, N1(00), P2(00), N2(00), P3a(00) and P3b(00) late component ERPs from Fz, Cz and Pz recording sites. The left side panels show conventional averages, i.e. the sums of all EEG segments time locked to the auditory stimuli application in a target detection task. Taken the total number of single trials as 100%, the bars show percentages of sweeps in which the components were identified, i.e. the elicitation rates. The results of the SCA are shown in the right side panels.

Figure 8 Left side panels show conventional average ERPs calculated from 40 single trials. Right side panels show average ERPs calculated from the same set of single trials by a new method of selective component averaging. The bars of the histograms show the elicitation rates.

We see drastic differences between the results of the two methods of averaging. The two important points are as follows.

First, compared with conventional averages, the outcomes of the SCA show significantly increased voltages. This indicates a crucial oversimplification of conventional ensemble averaging which assumes that ERP components in question are present in all single trials. Consequently, the average is divided by the number of single trials. By contrast, the SCA deals with actual numbers of identified components. As the histograms in Figure 8 show, these numbers are different for different components and smaller that numbers of trials. A crucial drawback of conventional averages is not an amplitude reduction itself, but its diverse and unpredictable effect on the morphologies of different components.

Second, significant differences in the waveform structures can be characterized as a better disclosure of the components by the SCA. The plots from the left side panels show that the peaking waveforms of conventional averages can be qualified as components just in the cases of N1, P2 and P3. By contrast, the SCA recognizes the whole complex of the late ERP components. Of the crucial importance is that model based algorithms of the SCA eliminate the temporal overlap of the P3a and P3b components and provide opportunities to treat these ERPs as separate entities.

To emphasize the differences in the reconstructions of the late positive components the 300 ms time instants are indicated in Figure 8 by the vertical blue lines. In all the SCA based component reconstructions the blue lines appear as borders between the P3a and P3b. In the reconstructions provided by conventional averaging the p3a component is hidden. We may propose that the temporal overlap of this component with the P3b causes increase of the P3b latency.

3.2.3 P3a and P3b

It is widely accepted that among the late component ERPs the most important marker of cognitive functions is the P3 (or P300) component, a positive ERP wave of approximately 300 ms latency which is evoked by rare, task-related ‘oddball’ stimuli that occur randomly among other meaningless ones [24]. The earliest contribution found that the P3 is made of two separate components, was the study of Squires and colleges who found a late-positive wave which peaks 60-80 ms earlier than the P3 [25]. They termed this component the “P3a” to distinguish it from the longer-latency classical P3 (latency 300-500 ms) that they la­beled “P3b”.

P3a generators are localized in cingulate, frontal and right parietal areas, and temporo-occipital regions [26]. Functionally, the P3a is associated with the orienting reflex and the automatic allocation of attention [27].

The temporal overlap of P3a with P3b creates numerous methodological complications in a reasonably precise quantification of the complex component composition of the late positive waveforms. If two positive peaks occur in the P3 latency range of average ERP, many investigators evaluate the complex average waveforms only for the maximum amplitude and latency which presumes a monolithic average P3. This approach a priori ignores the P3a as potentially meaningful endogenous component.

A simplified consideration of the late positive waveform as a monolithic P3 is also common for a single trial analysis using procedures with the correlation template [28]. This technique models P3 by a 2 Hz half-sine wave template or a similar simple waveform [29]. A possible double-peak composition of P3 is concealed by a narrow band filtering of single trial records with a typical bandpass from 0.5 to 4.43 Hz [28]. These methodological difficulties and related simplifications may explain why the P3a was not observed in ERP studies across different sub­ject and patient populations as consistently as was the larger and more prom­inent P3b.

Reliable recognition of both the P3a and P3b unfulfilled by the previous studies is achieved in our theory by creation of the model based approach to the component identification. General features of this technique are illustrated in the bottom panel of Figure 7 which shows separate models of both overlapping the P3a and P3b components.

Identification of the P3a-P3b complex instead of a single P3 has shown significant potential to enrich the diagnostic power of ERPs.

The study of the patients with borderline personality disorder (BPD) found distinctive disturbances in the P3a: abnormally enhanced amplitude, failure to habituate and a loss of temporal locking with the P3b. The reference to the normative age dependencies suggests that natural age-related decline in the P3a amplitude is reduced in BPD patients indicating the failure of frontal maturation [30]. Further investigation found that P3a amplitudes, but not latencies, were significantly larger in the right hemisphere compared with the left [31]. This symptom together with the failure of the habituation of the P3a is consistent with the deficient inhibitory activity. It is argued that these abnormalities may reflect impeded maturation of the fronto-medial processing systems of the brain.

From investigation of the ERPs in children and adolescents with functional neurological symptoms it is clear that conventional averaging does not disclose the late positive component as a composite potential [32]. In contrast single trial analysis using fragmentary decomposition provides clear reconstructions of the P3a and P3b. Both components show increased amplitudes in the patient groups.

3.3 Numerical Simulations of the EEG Microscale Basis

Conceptual basis for simulations of the microscale origins of the EEG is the proposition that EEG waveforms are produced by the ions released from cortical neurons to the extracellular space. The intracellular compartments in both the top and bottom panels in Figure 9 contain positive and negative ions which, in principle, can be released to the extracellular space. During the resting conditions the net charges of the particles released from the cells and captured by the cells are balanced. Given a local volume this means that net charges of the particles located outside of the cells in close proximity to the cell membranes fluctuate over some mean values.

Figure 9 The blue regions show particles in the extracellular space. Transition from the resting to the transient conditions creates STION, the ensemble of extracellular charges powerful enough to produce specific field potentials and EEG.

Transient conditions arise from the activation of a pertinent ensemble of the cells by a certain triggering events arising in the brain networks. Such processes induce massive exchange of the ions between the intracellular compartments and extracellular space. Due to a high electrical resistance of the membrane the movements of the ions inside the cells are not detectable in the extracellular space.

By contrast, the particles outside the cells are located in the conductive media. In the light of this fact, the generation of the EEG waves is associated in our theory with the transient changes of the net charges of the ions synchronously released from ensembles of closely located cells. Such population of the ions, i.e. a STION, is shown schematically in the top panel of Figure 9.

It is presumed that electric charge of the STION is enough powerful to induce the transient changes of electric fields in the extracellular medium including the scalp. This effect is schematically indicated by the red arrows in the top panel of Figure 9.

Important properties of these processes identified by our theory are the differences between the behavior of the positive and negative charges considered as the primary (P) and secondary (S) populations. A crucial effect following from the rule 2 is the blockage of the “birth” processes in the secondary particle population. Thus, during the transient conditions the particles in the primary population behave as a non-homogenous birth and death process while the particles in the secondary population behave as a non-homogenous death process.

Numerical simulations are unique tools to reconstruct the time courses of particle populations. The simulations illustrated in Figure 10 extend over the time interval from -10 to 70 ms with t=0 corresponding to the instant at which the transient starts.

Figure 10 Numerical reconstructions of the temporal evolution of particle populations in a typical trial. Resting conditions computed from -10 ms are switched at t=0 to the transient conditions. Blue, red, and black lines in the top panel show functions X_P(t), X_S(t) and X_N(t). The time courses of underlying transition probabilities are shown in the middle and bottom panels.

As an initial condition, an equal size N₀=50 was prescribed for both populations. The HWF was defined by the parameters σ=13.3 ms and η=26.2 ms.

Using defined parameters, it was important to satisfy (8) by the choice of a sufficiently small time interval Δ during which the expected change of the population size by more than one particle is negligibly small. Using Monte Carlo simulations for estimation of the numbers of particles crossing the membrane, different values of Δ were tested. In this way the value Δ=0.0001 ms was selected for the following simulations. The corresponding segmentation points over 80 ms time interval were t_i=i∙Δ with i taking values from 0 to 800,000.

The transition from the resting to transient condition was simulated as the change of the resting state transition probabilities (18) to the transient state transition probabilities in (19)-(21). The change occurs in a “smooth” fashion. This means that the sizes of the primary and secondary particle populations developed under the resting conditions serve as initial conditions for the transient regimes.

During the resting conditions (interval from -10 ms to t=0) the transport of particles between the primary and secondary populations is balanced. The particle numbers fluctuate over the mean value equal to N₀. The manner in which the transition from the resting to

transient conditions contributes to a rapid change of the net charge is due almost entirely to the change of the birth and death rates in the primary particle population. In the general case the size of the primary population is governed by the complex interplay of the birth and death transition probabilities. Onset of the transient conditions gives rise to both probabilities. Initially, from t=0 to the time instant indicated in the middle panel of Figure 10 by the arrow, the birth probability prevails over the death probability. At this stage nearly a tenfold increase of the size of the primary population occurs. After the peak, the death probabilities take a progressively larger share. As a result, the size of the primary population declines and returns to the initial conditions.

The effect of the transients in the secondary particle population on the net charge is minor and brief compared with the primary population. As shown in bottom panel of Figure 10, at t=0 the probability of the birth in the secondary population drops to zero. This blocks the particle transfers from the inside to outside of the cells.

In order to decide whether the mass effects of particle movements are sufficient to a full account for dynamics of macroscale processes, it is necessary to compare the results of computer simulations with the theory. We consider the transient conditions starting at t=0. Since ψ(0)=0, we use ψ₀=ψ_P(0)=ψ_S(0) as the initial condition for theoretical solutions. The corresponding initial condition for numerical simulations is $\mathbf{X}_{P\ }\left(0\right)=\mathbf{X}_S\left(0\right)=N_0$.

Defined η₀ and N₀ parameters allow us to create dimensionless functions ψ^*(t)=ψ(t)/ψ₀ and $\mathbf{X}_N^\ast\left(t\right)={\mathrm{X}_N\left(t\right)}/{N_0}$. After this normalization we can directly compare numerically calculated $X_n^*t$ with theoretical ψ^*(t).

Taking values of σ and η parameters from the previous simulations, the samples of $X_n^*t$ computed with N₀ equal to 10, 20, and 100 particles are illustrated by the colored lines in the top panel of Figure 11.

It is evident that an increase in the number of particles brings single trial samples of ψ^*(t) to better proximity with the theoretical model. Thus, the role of deterministic factors in statistical samples of $X_n^*(t)$ increases with increase in the number of particles involved.

To emphasize this tendency we have estimated absolute values of residuals between ψ^*(t) and $X_n^*(t)$. Given N₀, 20 single trials were calculated. In each trial the residuals were estimated in 100 equidistant points. From these measurements the average residuals were calculated. The points in the plot in the bottom panel of Figure 11 show dependency of average residuals on the number of particles N₀. The increase of the particle numbers to several hundredths makes single trials perfectly identical to the theoretical solution. This provides convincing evidence that, under the transient conditions the particle behavior in both primary and secondary particle populations develops as an amalgamation of deterministic and stochastic processes. These are two broadly defined classes of processes with specific properties. We now consider such situation from the point of view of deterministic chaos, i.e. processes that share attributes with both deterministic signals and stochastic processes.

Figure 11 Coloured lines in the top panel show typical computer reconstructions of the temporal evolution of particle populations with different initial sizes of 10, 20, and 100 particles. The black line is the theoretical solution η^*(t). Data points in the bottom panel show average residuals between the theoretical and numerical solutions calculated for the particle populations with different initial sizes.

3.4 Transient Deterministic Chaos

Research into chaos phenomena is progressing steadily, extending to a diverse range of applications. No universal definition of the system producing chaos has been accepted in such situations of evolving hypotheses and concepts. Necessary conditions are the non-linearity of the system generating the chaos and its sensitivity to the initial conditions.

In terms of the specifics of biological systems, in most cases the equations of the sources of deterministic chaos are unknown [33]. Accordingly, the identification of chaotic behavior is usually based on the analysis of empirical time series. Most popular methods for characterizing chaotic behaviors are the Lyapunov exponents and Kolmogorov-Sinai entropy. These are essentially indirect methods which do not define the equations but diagnose whether the process in question is chaotic or non-chaotic. This limitation applies equally to other diagnostic methods, among which are the Poincare map, correlation dimensions, fractal dimensions, and attractor reconstructions from the time series.

Our approach actually does not need such indirect methods because the theory provides complete models of the HWF generation on both the macro- and micro-scales. The two major criteria which allow us to qualify the system producing these processes as chaotic are as follows.

1) The non-linearity of the system producing HWF.

2) Sensitive dependence of the system behavior on initial conditions.

The nonlinearity of the system is evident on the macroscopic scale where the system producing HWF is expressed by the following system of nonlinear differential equations:

\[ \begin{gather} {d\left[\mathrm{\psi}_P\left(t\right)\right]}/{dt=\left[-{\left(\eta-t\right)}/{\sigma^2}\right]\mathrm{\psi}_P\left(t\right),}\\ d{\left[\mathrm{\psi}_S\left(t\right)\right]}/{dt}=\left[-{\left(\eta+t\right)}/{\sigma^2}\right]\mathrm{\psi}_S\left(t\right),\\ \mathrm{\psi}\left(t\right)=\mathrm{\psi}_P\left(t\right)-\mathrm{\psi}_S\left(t\right). \end{gather}\tag{28} \]

We may qualify the system as a non-autonomous nonlinear system, that is, a system driven by the time-varying processes.

The chaotic behavior is hidden in the EEG dynamics and is only recognizable by specific Gaussian profiles of the HWF shapes. It is important to note that specific statistical background of such patterns allowed us to uncover the rules 1 and 2 that govern chaotic phenomena taking place on the microscale. An important additional factor for the confirmation of chaos is that the chaotic system defined by our theory arises solely from the equations of the system without influence of additional factors.

The strong dependency of the chaotic processes in question on initial conditions is seen from the computer simulations illustrated in Figure 12. The parameters are the same as in the simulations shown in Figure 11. The simulations started from 10 ms before the arrival of the triggering event. At that instant both P and S populations contain 50 particles. Because on the initial 10 ms interval the sizes of P and S populations were changing randomly as simple BDPs, the initial sizes of the whole population were slightly different in different trials. As the results of simulations indicate, these minor differences in the initial conditions lead to significantly different future behavior of the realizations of X_N(t).

Figure 12 Using the same conditions as the simulations in Figure 10, illustrates the sensitive dependence of X_N(t) realizations in different trials (coloured lines) on initial conditions. The dotted line is a theoretical solution.

Another important feature of the chaotic processes is the short term predictability. It is clear from the simulation results in Figure 12 that despite substantial differences in the trajectories of X_N(t) realizations, it is possible to give rough estimates of some parameters, for example the times at which the trajectories reach the maximums. We may also claim that the features of predictability are also incorporated in the time courses of X_N(t) realizations because the averages of these trajectories from different trials converge to the analytical solutions.

An advantage of the analytical methods of our theory is the possibility to make a clear distinction between chaos and pure stochasticity. During the resting periods we deal with purely stochastic processes. The start of the HWF generation is specified by the time instant when the triggering event arrives. Thus the deterministic chaos develops as a transient process. We label this amalgamation of deterministic and stochastic processes with the term “transient deterministic chaos”.

4. Discussion

The major assumption of EEG models supported by differential and integral equations of the classical physics presumes that the EEG dynamics can be deduced using the superposition of membrane potentials produced by the underlying neuronal ensembles. Severe difficulties are created by intractably huge multiplicity of the cellular elements and an insufficient knowledge of their parameters. Under such uncertainty little confidence can be placed on particular solutions from a large number of different models.

The idea of an alternative phenomenological approach using statistical considerations has been put forward by Elul [34]. Based on the analysis of the synchronization of EEG sources, Elul proposed that evolution of the brain waves may be governed by statistical regularities following from the central limit theorem. Thus, the EEG may simply be accounted in terms of a normally distributed output resulting from a combination of the activity of many independent neuronal generators.

The mathematics behind these propositions has never been worked out in detail. The major difficulty is necessity to support probabilistic propositions by adequate models of pertinent microscale events. To the best of our knowledge, this study is the first to provide an empirically grounded probabilistic model of the microscale events and bridge them with the EEG dynamics. The crucial step is the change of basic microscale units from the continuous time membrane potentials to elementary particles, positive and negative ions.

This permits a great deal to be inferred about the mass effects of very complicated cellular structures without even mentioning the cells and channels, or even being very explicit about internal makeup. A fundamental point is the consideration of deterministic components of EEG potentials as the limits of the underlying microscale processes. Vanishingly small impact of individual particles to the generation of the macroscopic scale processes reduces the problem to the study of the limiting behavior of large numbers of random variables. This is a classical problem in the probability theory [8]. Some purely mathematical aspects of how deterministic potentials, such as membrane potentials, emerge as a limit from the underlying stochastic sources have been recently considered by Austin [35].

Referring to various probabilistic models of ion transport [36], the new feature of the probabilistic framework of our theory is the implementation of non-homogenous BDPs with time dependent rates of birth and death. There are at least two important consequences of this extension of the probabilistic tools employed.

First, it became possible to extend purely stochastic models of microscale machinery by emergence of the field potentials with deterministic trends which we classify as the transient deterministic chaos.

Second, a link was established by (6) between the microscale transients and the macroscale components of the mass potentials.

This link provides means to estimate the microscale parameters on empirical grounds by consideration of the EEGs as non-stationary processes. On the macroscopic scale we use fragmentary decomposition, a model based method of non-stationary signal analysis which decomposes EEG into the functionally meaningful components on an adaptive segmentation basis. At this step we follow conventional empirical understanding of the functional components as positive and negative peaking waveforms that are observed in the EEG time course. The possibility of transferring these empirical pieces of ongoing EEG to the universal analytical models has been originally supported by the results of the single trial analysis of the late component ERPs recorded in a conventional auditory oddball paradigm [16]. It appeared that though different components are associated with different neural sources their normalized amplitude spectra are practically identical. Such results are somewhat difficult to reconcile with our physiological intuitions and concepts about different structure and functions of underlying cellular systems.

As a clue to the understanding of this seemingly paradoxical situation we need to refer to the role of the Gaussian function as a normal distribution, which perfectly describes statistical regularities independently of the physical nature of the objects from which the samples are drawn. In accordance with the central limit theorem, any process of random sampling tends to produce a normal distribution of sample values, even if the whole population from which the samples are drawn does not have a normal distribution. It was therefore reasonable to propose that appearance of a Gaussian amplitude spectrum indicates that mass effect produced by multiple microscale events is governed by the rules of the central limit theorem. In order to accept this proposition it was necessary to reveal how statistical factors identified in the frequency domain work in the time domain.

The time domain counterpart of the frequency domain Gaussian function is another Gaussian function defined over an infinite time scale. Such function is not compatible with the semi-infinite extent of the transient potentials. The ability of our theory to deal with transients is due to the specific composition of the HWF from the two Gaussian functions with identical parameters. We regard these Gaussian functions as products of separate but operationally linked particle populations with identical birth and death rates. The primary population contains particles of the same polarity as the macroscopic component. The particles from the secondary population have an opposite polarity. Thus, an important finding of our study is that during development of the transient deterministic chaos on the microscopic scale the two types of participating particles with opposite polarities are recognizable from the macroscopic scale.

We may state that half-wave function ψ(t) emerging as the macroscopic scale effect from the synchronized chaotic ion movements on the microscopic scale can be regarded as a universal building block from which various mass potentials are composed. In terms of quantum theories this effect can be qualified as universality, the concept that there are properties for a large class of systems that are independent of the exact structural and functional details of the system. This probabilistic basis implemented by our theory provides the tools which are able to reduce intractably huge numbers of microscale variables to a universal macroscale object in the form of a HWF with the fewest parameter set (κ, σ, and η). Being limit distributions for the sums of random variables, these parameters discriminate those aspects of the molecular machinery that are significant on the global scale from those that are not.

The possibility to decompose a complex signal into a set of universal “building blocks” plays an important role in different fields of physics and biology. Thus, synaptic physiology has greatly benefited from the quantum analysis which deduced a miniature end potential as a building block from which the postsynaptic potentials are composed. An important outcome is that this concept is applicable to various types of synaptic junctions with different molecular organization. This is a remarkable feature of the quantal universality which applies equally to our treatments.

To our knowledge, our study is the first to put the models of various types of EEGs and ERPs on a common theoretical and computational basis. Using HWF as a universal building block, we obtain a general model of the mass potential in the form of (26). This equation suggests that a mass potential evolves as the succession of transient components induced at consecutive time instants. Each component is a HWF described by the system of non-linear equations (27) with specific sets of estimated parameters.

Conventional procedures of parameter estimation consider a mass potential as the set of components in the form of peaking waves. On consideration of the component as being not just a peak in the waveform, but a whole deflection (positive or negative) with specific shape, the HWF is an adequate analytical model of such an entity. The timing (τ) and magnitude (κ) of HWF correspond to conventional component parameters. By contrast, no parameters were accepted so far as adequate measures of the component shapes. In this regard the σ and η parameters may serve as universal shape parameters with the potential to provide some specific information about underlying microscale processes.

However the range of the microscale phenomena to which our models of single trial ERP components are relevant is limited in two respects: in the first place, they do not contain information about spatial distributions of ERP sources, and in the second, they apply in their present form only to the summary effects produced by large ensembles of closely located ions.

To approach these problems we are planning to consider ensembles of quantum models supported by EEG data from multiple cortical sites.

5. Appendix. Time-Frequency Analysis Using the SBF Algorithm

The time-frequency analysis is applied to the EEG which is presented as the following sum of EHWs:

\[ v\left(t\right)=\sum_{i=1}^{N}{w_i\left(t-\tau_i\right).} \]

EHW w_i (t) in this formula reproduces EEG v(t) in the interval between the segmentation points τ_i-1 and τ_i.

The goal of the time-frequency analysis is to transform each w_i(t) to the frequency domain using numerical calculation of the exponential finite Fourier transform. Procedures applied to various HWs are universal. Thus, we omit subscripts in EHW descriptions and consider w(t) as a target for the following transformations. The corresponding complex spectrum is defined by the following Fourier transform

\[ W\left(i\omega\right)=\int_{0}^{T}w\left(t\right)\mathrm{exp} \left(-i\omega t\right)dt. \]

Numerical estimation of Fourier transforms is usually based on procedures employing various algorithms of the fast Fourier transform. The latter is supported by a Fourier series model of the data, i.e. the addressees are periodic signals. This distinction with the Fourier integrals is a troublesome problem when aperiodic signals of short duration, like EHW, are analyzed. When a waveform is not periodic in time, the spectral leakage may cause significant distortions of the frequency characteristics. The remedies of windowing and zero-padding usually introduce problems of their own.

As an adequate tool of the spectral analysis of HWs we use the SBF algorithm [19]. The algorithm is an original version of Filon-type methods that provide maximum precision in the estimation of trigonometric integrals using interpolation polynomials of different degrees.

Using the SBF algorithm we deal with numerical calculations of the real and imaginary parts of the complex spectrum W(iω):

\[ \mathrm{W}_C\left(\omega\right)=Re\left[\mathrm{W}\left(i\omega\right)\right]=\int_{0}^{T}\mathrm{w}\left(t\right)\mathrm{cos}(ωt)dt, \]

\[ \mathrm{W}_S\left(\omega\right)=Im\left[\mathrm{W}\left(i\omega\right)\right]=\int_{0}^{T}\mathrm{w}\left(t\right)\mathrm{sin}(ωt)dt, \]

where W_C(ω) and W_S(ω) are the finite cosine and sine Fourier transforms of w(t) and [0, T] is the interval of integration.

For numerical calculations w(t) is specified by its sampled values w_i=w(t_i), where t_i=iΔ, i takes values from 0 to N and t_N=T. Over these samples a piece-wise linear polynomial h(t) is created.

Figure 13a illustrates the principle of the piecewise-linear interpolation. Given w(t) (blue line), the interpolant, h(t), is the broken red line created by joining the nodal points 0, 1, 2, 3 and 4 by the straight lines. In the same fashion the interpolation can be continued for any number of the following nodal points.

Figure 13 Graphical illustration of the SBF algorithm.

A specific aspect of the SBF algorithm is that h(t) is decomposed into the weighted sum of similar basis functions of the following form

\[ \mathrm{h}\left(t\right)=\sum_{i=0}^{N-1}{a_ir\left({t}/{t_{i+1}}\right),}\tag{A1} \]

where a_i is the interpolation coefficient and r(t) is a similar basis function defined in the form of the triangular basis function (TBF):

\[ \mathrm{r}\left(t\right)=\left\{\begin{matrix}1-t\ \ \ \mathrm{if}\ \ \ 0\le t\le1\\0\ \ \mathrm{otherwise\ \ \ \ \ \ \ \ \ \ \ \ \ }\\\end{matrix}\right. \]

The simple geometric form of TBF (unit right-angled triangle) is illustrated in Figure 13.

The interpolant, h(t), h(t) is created as a piece-wise linear polynomial which satisfies the interpolation condition:

\[ h_i=w_i for\ \ \ i=0,1,\ldots,N-1,\tag{A2} \]

where $h_i=h(t_i)$.

A general procedure of the evaluation of the interpolation coefficients consists of writing Eq (A1) for each interpolation point and finding the solution of the resulting system of N linear equations by conventional matrix methods.

While this set of linear equations is easily solved by standard methods, there is a further simplification possible due to the fact that $r\left({t_m}/{t_n}\right)=0$ for m≥n. To make this mode of interpolation as intuitive as possible, we use as an intermediate element the hat function, one of the geometrical objects employed by the method of finite-elements. The hat function is defined as

\[ \vartheta_{i}(t)=\left\{\begin{array}{l}\frac{t-t_{i-1}}{t_{i}-t_{i-1}}, \text { if } t_{i-1} \leq t<t_{i} \\\frac{t_{i+1}-t}{t_{i+1}-t_{i}}, \text { if } t_{i} \leq t<t_{i+1} \\0, \text { otherwise }\end{array}\right. \]

on the mesh t₀<t₁<...<t_i<...<t_N.

Triangles a1c, b2d, c3e in Figure 13a and abc in Figure 13b are the examples of the weighted hat functions.

The interpolation capability of the hat function is supported by the two properties. First, ϑ_i(t) vanishes everywhere except on the two subintervals to which t_i belongs. Second, ϑ_i(t) is unity at the node i and zero at all other nodes.

Using the hat function, the interpolation condition (A2) may be presented in the form

\[ \mathrm{h}\left(t\right)=w_0r\left(\frac{t}{t_1}\right)+\sum_{i=1}^{N-1}{w_i\vartheta_i\left(t\right).}\tag{A3} \]

Now the interpolation coefficients are equal to the samples of w(t). This simplification is the result of the change of the basis function from the TBF in Eq (A1) to the hat function in Eq (A3). The geometric principle found useful in this connection may be seen by reference to Figure 13. Each nodal point from 0 to 3 in Figure 13a is linked with a particular triangle. The nodal point 0 is the top of the right-angled triangle which corresponds to the first terms in (A1) and (A3). The nodal points from 1 to 3 are the tops of the hat functions multiplied by the interpolation coefficients (triangles a1c, b2d and c3e). The graph in Figure 13b shows that the hat function (triangle abc with bd=1) may be decomposed into the sum of the three TBFs triangles ohc, ofd and oga). Thus, in the interval from t_i-1 to t_i+1 the corresponding hat function is defined as

\[ \vartheta_i\left(t\right)=\alpha_ir\left({t}/{t_{i+1}}\right)-\beta_ir\left({t}/{t_i}\right)+\gamma_ir\left({t}/{t_{i-1}}\right) \]

where

\[ \begin{align}\alpha_i=\frac{t_{i+1}}{\Delta t_{i+1}}\ \ \ \ \ \left(0\le i\le N-1\right),\ \ \ \ \ \ \ \ \ \ \ \ \ \\\beta_i=t_i\frac{\Delta t_{i+1}+\Delta t_i}{\Delta t_{i+1}\Delta t_i}\ \ \ \ \ \left(1\le i\le N-1\right),\\\gamma_i=\frac{t_{i-1}}{\Delta t_i}\ \ \ \ \ \left(2\le i\le N-1\right),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\\end{align} \]

and Δt_i=t_i-t_i-1.

In these terms

\[ h\left(t\right)=w_0r\left({t}/{t_1}\right)+w_1\left[\alpha_1r\left({t}/{t_2}\right)-\beta_1r\left({t}/{t_1}\right)\right]+\sum_{i=2}^{N-1}{w_i\vartheta_i\left(t\right).} \]

Comparison with Eq (A1) shows that interpolation coefficients are readily found to be

\[ a_i=\alpha_iw_i-\beta_{i+1}w_{i+1}+\gamma_{i+2}w_{i+2}\ \ \ \ \ \ \ \ \ \ \left(i=0,\ldots,N-1\right), \]

where α₀=1 and it is assumed that w_o=0 and $w_{N+1}=0$.

The system of simultaneous linear equation has thus been reduced to a set of independent linear equations. Now h(t) is described in terms of a linear combination of TBFs, the finite Fourier cosine- and sine- integrals of which are given by

\[ \mathrm{R}_C\left(\omega\right)=\int_{0}^{1}{\mathrm{r}\left(t\right)\mathrm{cos} \left(\omega t\right)dt=\frac{1-\mathrm{cos} \left(\omega\right)}{\omega^2},} \]

\[ \mathrm{R}_S\left(\omega\right)=\int_{0}^{1}\mathrm{r}\left(t\right)\mathrm{sin} \left(\omega t\right)dt=\frac{\omega-\mathrm{sin} \left(\omega\right)}{\omega^2}. \]

These functions are illustrated in Figure 13d.

According to the similarity theorem of the theory of Fourier transforms, the compression of the abscissa in the time domain corresponds to the expansion of the abscissa plus contraction of the ordinate in the frequency domain. These relationships reduce the entire issue of the transform calculations to some standard manipulations with relatively simple functions R_C(ω) and R_S(ω). The cosine- and sine- Fourier integrals from h(t) obtain the following representations

\[ \mathrm{W}_C\left(\omega\right)=\sum_{i=0}^{N-1}{a_i t_{i+1}\frac{1-\mathrm{cos} \left(t_{i+1}\omega\right)}{\left(t_{i+1}\omega\right)^2},} \]

\[ \mathrm{W}_\mathrm{S}\left(\omega\right)=\sum_{i=0}^{N-1}{a_it_{i+1}\frac{\omega t_{i+1}-\mathrm{sin} \left(t_{i+1}\omega\right)}{\left(t_{i+1}\omega\right)^2}.} \]

The corresponding amplitude spectrum W(ω) and the phase function δ(ω) are defined as:

\[ \mathrm{W}\left(\omega\right)=\left|\mathrm{W}\left(i\omega\right)\right|=\sqrt{\mathrm{W}_C^2\left(\omega\right)+\mathrm{W}_S^2\left(\omega\right),} \]

\[ \mathrm{\delta}\left(\omega\right)=\mathrm{arctg}\left[W_S\left(\omega\right)/W_C\left(\omega\right)\right]. \]

These are continuous functions of angular frequency. Accordingly, the values of the frequency domain characteristics can be calculated for arbitrary sets of points. We employ a logarithmic frequency scale for calculations of amplitude spectra and a natural frequency scale for calculations of phase functions.

In the case of a logarithmic scale, the frequency characteristics were calculated for angular frequencies ω_i=ω₀C^i-1 (i=1, …, N), where ω₀ is initial angular frequency and C>1 is the parameter that defines the sampling rate. For selection of this parameter it is convenient to use the formula C=exp(ln10/N_D), where N_D is the number of samples per decade. In the case of a natural frequency scale, ω_i=ω₀+iΔω (i=0, …, N), where Δω is the discretization step.

Author Contributions

Dmitriy Melkonian did all the research work of this study.

Competing Interests

The authors have declared that no competing interests exist.