Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

A harmonic function method for EEG source reconstruction

  • In this paper we study a harmonic function method for dipolar source reconstruction, and implemented the numerical simulations. We propose a new error estimate and provide a rigorous proof of the estimate. Then, we validate our method in computer-simulated data and study its numerical stability in different noise levels. It is shown that the harmonic function method can be used to quickly and accurately locate the active regions in EEG source reconstruction.

    Citation: Hongguang Xi, Jianzhong Su. A harmonic function method for EEG source reconstruction[J]. Electronic Research Archive, 2022, 30(2): 492-514. doi: 10.3934/era.2022026

    Related Papers:

    [1] Jiyu Sun, Jitao Zhang . Muti-frequency extended sampling method for the inverse acoustic source problem. Electronic Research Archive, 2023, 31(7): 4216-4231. doi: 10.3934/era.2023214
    [2] Xinlin Cao, Huaian Diao, Jinhong Li . Some recent progress on inverse scattering problems within general polyhedral geometry. Electronic Research Archive, 2021, 29(1): 1753-1782. doi: 10.3934/era.2020090
    [3] Youzi He, Bin Li, Tingting Sheng, Xianchao Wang . Generating geometric body shapes with electromagnetic source scattering techniques. Electronic Research Archive, 2020, 28(2): 1107-1121. doi: 10.3934/era.2020061
    [4] Jie Wang . A Schwarz lemma of harmonic maps into metric spaces. Electronic Research Archive, 2024, 32(11): 5966-5974. doi: 10.3934/era.2024276
    [5] Weishi Yin, Jiawei Ge, Pinchao Meng, Fuheng Qu . A neural network method for the inverse scattering problem of impenetrable cavities. Electronic Research Archive, 2020, 28(2): 1123-1142. doi: 10.3934/era.2020062
    [6] Xiaoping Fang, Youjun Deng, Zaiyun Zhang . Reconstruction of initial heat distribution via Green function method. Electronic Research Archive, 2022, 30(8): 3071-3086. doi: 10.3934/era.2022156
    [7] Yao Sun, Lijuan He, Bo Chen . Application of neural networks to inverse elastic scattering problems with near-field measurements. Electronic Research Archive, 2023, 31(11): 7000-7020. doi: 10.3934/era.2023355
    [8] Yujie Wang, Enxi Zheng, Wenyan Wang . A hybrid method for the interior inverse scattering problem. Electronic Research Archive, 2023, 31(6): 3322-3342. doi: 10.3934/era.2023168
    [9] Hao Wang, Wei Yang, Yunqing Huang . An adaptive edge finite element method for the Maxwell's equations in metamaterials. Electronic Research Archive, 2020, 28(2): 961-976. doi: 10.3934/era.2020051
    [10] Youjun Deng, Hongyu Liu, Xianchao Wang, Dong Wei, Liyan Zhu . Simultaneous recovery of surface heat flux and thickness of a solid structure by ultrasonic measurements. Electronic Research Archive, 2021, 29(5): 3081-3096. doi: 10.3934/era.2021027
  • In this paper we study a harmonic function method for dipolar source reconstruction, and implemented the numerical simulations. We propose a new error estimate and provide a rigorous proof of the estimate. Then, we validate our method in computer-simulated data and study its numerical stability in different noise levels. It is shown that the harmonic function method can be used to quickly and accurately locate the active regions in EEG source reconstruction.



    Neuronal activities generate the electrical current in the brain, and further result in the potential changes over the scalp. Electroencephalography (EEG) is a technique used to record the potential changes on the scalp. Even though fMRI, PET, MEG and other brain-imaging tools are widely used in brain research, they are limited by low spatial/temporal resolution, cost, mobility and suitability for long-term monitoring. For example, fMRI has the advantage of providing spatially-resolved data, but suffers from an ill-posed temporal inverse problem, i.e., a map with regional activations does not contain information about when and in which order these activations have occurred [1]. In contrast, EEG signals have been successfully used to obtain useful diagnostic information (neural oscillations and response times) in clinical contexts. Further, they present the advantage to be highly portable, inexpensive, and can be acquired at the bedside or in real-life environments with a high temporal resolution. Because of the lack of significant patient risks, EEG is additionally suited for long-term monitoring.

    EEG offers the possibility of measuring the electrical activity of neuronal cell assemblies on the sub-millisecond time scale [2,3,4]. EEG source imaging further identifies the positions or distributions of electric fields based on EEG signals collected on the scalp [5]. This new tool is widely used in cognitive neuroscience research, and has also found important applications in clinical neuroscience such as neurology, psychiatry and psychopharmacology [6,7]. In cognitive neuroscience, the majority of the studies investigate the temporal aspects of information processing by analyzing event related potentials (ERP). In neurology, the study of sensory or motor evoked potentials is of increasing interest, but the main clinical application concerns with the localization of epileptic foci. In psychiatry and psychopharmacology, a major focus of interest is the localization of sources of certain EEG frequency bands. Localizing the activity sources of a given scalp EEG measurement is achieved by solving the so-called inverse problem [8]. These kinds of inverse problems are usually ill-posed and their solutions are non-unique [9,10].

    Leahy et al. [11] investigated the accuracy of forward and inverse techniques for EEG and MEG dipole localization using a human skull phantom. El Badia and Ha-Duong [12] established an algebraic method to identify the number, locations and moments of electrostatic dipoles in 2D or 3D domain from the Cauchy data on the boundary. Chafik et al. [13] further provided an error estimate without proof. Nara and Ando [14] provided a new projective method for 3D source reconstruction by projecting the sources onto a Riemann sphere. Kang and Lee [15] proposed an algorithm for solving the inverse source problem of a meromorphic function and apply their method to an electrical impedance tomography (EIT) problem. El Badia [16] established a uniqueness result and a local Lipschitz stability estimate for an anisotropic elliptic equation, assuming that the sources are a linear combination of a finite number of monopoles and dipoles. The author also proposed a global Lipschitz stability estimate for dipolar sources. Baratchart et al. [17] solved the inverse source problem by locating the singularities of a meromorphic function from the 2D boundary measurements using best rational or meromorphic approximations.

    Chung and Chung [18] proposed an algorithm for detecting the combination of monopolar and multipolar point sources for elliptic equations in the 2D domain from the Neumann and Dirichlet boundary data. Kandasmamy et al. [19] proposed a novel technique, called "analytic sensing", to estimate the positions and intensities of point sources in 2D for a Poisson's equation. Analytic sensing also used the reciprocity gap principle, but with a novel design of an analytic function which behaved like a sensor. The authors evaluated their estimation accuracy by Cramér-Rao lower bound. Nara and Ando [20] proposed an algebraic method to localize the positions of multiple poles in meromorphic function field from an incomplete boundary. They investigated the accuracy of the algorithm for the open arc or the closed arc, and for the arc enclosing the poles or not enclosing the poles. El Badia and Nara [21] established the uniqueness and local stability result for the inverse source problem of the Helmholtz equation in an interior domain, assuming the source is composed of multiple point sources.

    Clerc et al. [22] applied best rational approximation techniques in the complex plane to EEG source localization and offered stability estimates. Mdimagh and Ben Saad [23] identified the point sources in a scalar problem modeled by Helmholtz equation, using reciprocity gap principle and assuming the sources are harmonic in time. They proved local Lipschitz stability by two methods: one was derived from the Gâteaux differentiability, and the other used particular test functions in the reciprocity gap functional. Vorwerk et al. [24] studies the important role of head tissue conductivity in EEG dipole reconstruction. Rubega et al. [25] estimated EEG source dipole orientation based on singular-value decomposition. Michel and Brunet provided a thorough review on EEG source imaging. There exist several reconstruction methods, such as minimum norm estimates (MNE) [26], low resolution electrical tomography (LORETA) [27,28] or multiple-signal classification algorithm (MUSIC) [29,30], etc.

    Recently, Muñoz-Gutiérrez et al. [31] managed to improve the accuracy of EEG source reconstruction by decomposing the EEG signals into frequency bands with different methods, such as empirical mode decomposition (EMD) and wavelet transform (WT). Kaur et al. [32] presented a new method of EEG source localization using variational mode decomposition (VMD) and standardized the low resolution brain electromagnetic tomography (sLORETA) inverse model. Their VMD-sLORETA model could locate EEG sources in the brain in a very accurate way. Oikonomou and Kompatsiaris [33] developed a novel Bayesian approach for EEG source localization. They incorporated a new sparse prior for the localization of EEG sources with the variational Bayesian (VB) framework and obtained more accurate localization of EEG sources than state-of-the-art approaches.

    In our study we need new methods to detect small changes in EEG source for which dipole methods have advantage. We followed the analytic dipole method by El Badia and Ha-Duong [12] and derived a new error estimate for this source localization method. We provided a mathematical proof of this estimate. We then use simulated data to validate the method. The simulation results support our error estimation, which has a different distance power than a similar error estimate in [22].

    We organize the rest of the paper as follows. In section 2, we introduce the method and its formulation. In section 3, we provide the error estimate of an inverse EEG source localization problem in a bounded domain and its mathematical proof. In section 4, we use simulated data to valid the method and error estimate. A brief conclusion and discussion is in Section 5.

    The electric field E is the negative gradient of the potential u.

    E=u. (2.1)

    The quasi-static approximation means all time derivatives in the equation are set to zero. By quasi-static approximation of Maxwell equation ×HDt=J, we have

    ×H=J

    where H is the magnetizing field, J is the total current density, and D is the displacement field.

    Since the divergence of a curl is always zero, we have

    (×H)=J=0.

    EEG problem can be modeled by a Poisson equation.

    (σu)=(σE)=(JJp)=J=0Jp=Jp=F,

    where σ is the conductivity, Jp is the primary current density, and F is the source term.

    If we assume the source is composed of a finite number of point charges, then by linear combination, we have

    F=mk=1qkδ(rrk), (2.2)

    where m is the number of point charges, qk are values of charges, and rk are the locations of the point charges.

    If we assume the source is composed of a finite number of dipoles, we have

    F=mk=1pkδ(rrk),

    where m is the number of dipoles, pk are the moments (or strengths) of the dipoles, and rk are the centers of dipoles.

    The dipolar source reconstruction problem can be viewed as a Poisson problem.

    Δu=mk=1pkδ(rrk) in Ω, (2.3)
    u=f on Γ, (2.4)
    uν=φ on Γ, (2.5)

    where f and φ are known, and ν is the outer unit normal vector.

    We will use the concept of reciprocity gap functional [34]:

    R(v)=uν,vH1/2(Γ),H1/2(Γ)u,vνH1/2(Γ),H1/2(Γ)=φ,vH1/2(Γ),H1/2(Γ)f,vνH1/2(Γ),H1/2(Γ), (2.6)

    where v is a harmonic function in Ω:

    vH(Ω)={wH1(Ω)Δw=0}. (2.7)

    By Green's formula, we have

    R(v)=mk=1pkv(rrk),vH(Ω). (2.8)

    Let m be the number of dipoles in the brain. Assume mM in our problem, i.e., there is an upper bound for the number of dipoles.

    Let us consider the harmonic polynomials

    vj(x,y)=(x+iy)j,jN.

    Then, in 2D case

    R(vj)=mk=1pkvj(rk)=mk=1[pk1pk2](xk+iyk)j=mk=1[pk1pk2][x(x+iy)jy(x+iy)j]x=xk,y=yk=mk=1[pk1pk2][j(xk+iyk)j11j(xk+iyk)j1i]=mk=1[pk1pk2][1i]j(xk+iyk)j1=jmk=1(pk1+ipk2)(xk+iyk)j1.

    We define

    βj:=R(vj)j=Mk=1(pk1+ipk2)(xk+iyk)j1,j=1,2,...,2M1. (2.9)

    Let

    ηj=[βjβj+1βj+M1]CM,1jM, (2.10)

    and

    Zi=[ηi,ηi+1,...,ηi+M1]=[βiβi+1βi+M1βi+1βi+2βi+Mβi+M1βi+Mβi+2M2],iN.

    Then,

    Z1=[η1,η2,...,ηM]=[β1β2βMβ2β3βM+1βMβM+1β2M1].

    The number m of dipoles is estimated as the rank of Z1.

    Now we can reduce the size of the matrix by recalculating βj and ηj with M replaced by m. Then, the m vectors η1,...,ηm are independent.

    To get the estimates of the positions we need to construct an m×m matrix T such that ηj+1=Tηj,j=1,...,m. Then,

    [η2,...,ηm+1]=T[η1,...,ηm].

    So,

    T=[η2,...,ηm+1][η1,...,ηm]1=[β2β3βm+1β3β4βm+2βm+1βm+2β2m][β1β2βmβ2β3βm+1βmβm+1β2m1]1=Z2Z11.

    The positions of dipoles are estimated as the eigenvalues of T.

    We now show that the eigenvalues of T are the positions of dipoles. Let us first look at an example η2=Tη1.

    Tη1=T[β1β2βm]=T[p1+p2++pmp1S1+p2S2++pmSmp1Sm11+p2Sm12++pmSm1m]=p1T[1S1Sm11]+p2T[1S2Sm12]++pmT[1SmSm1m],

    where pk=pk1+ipk2,k=1,2,...,m is the moment and Sk=xk+iyk,k=1,2,...,m is the position.

    η2=[β2β3βm+1]=[p1S1+p2S2++pmSmp1S21+p2S22++pmS2mp1Sm1+p2Sm2++pmSmm]=p1S1[1S1Sm11]+p2S2[1S2Sm12]++pmSm[1SmSm1m],

    where pk=pk1+ipk2,k=1,2,...,m is the moment and Sk=xk+iyk,k=1,2,...,m is the position.

    Since [1S1Sm11],[1S2Sm12],...,[1SmSm1m] are independent and the results are similar for ηj+1=Tηj,j=1,2,...,m, we know S1,S2,...,Sm are just the eigenvalues of T.

    Now the question is how to get T. Only η1 and η2 are not enough to determine T because vectors have no inverse. So, we use the redundant information to construct the matrices Z1 and Z2 such that T=Z2Z11, where Z1 is invertible because η1,...,ηm are independent.

    To estimate the moments of dipoles we will write Eq (2.9) in matrix form. Notice that now we use m instead of M.

    [β1β2βm]=[S01S02S0mS11S12S1mSm11Sm12Sm1m][p1p2pm], (2.11)

    where pk=pk1+ipk2,k=1,2,...,m is the moment and Sk=xk+iyk,k=1,2,...,m is the position.

    We can write Eq (2.11) in matrix form

    b=Sp, (2.12)

    where b=[β1β2βm],S=[S01S02S0mS11S12S1mSm11Sm12Sm1m], and p=[p1p2pm]. Then, the moments of dipoles in 2D are estimated as

    p=S1b. (2.13)

    Equation (2.13) works in the ideal case of no noise. In reality, due to the noise in the measurements and in the sources, we need find a linear operator L to estimate the moments, i.e.,

    ˜p=Lb (2.14)

    where ˜p represents the estimates of the moments, and b represents the quantities obtained from the measurements.

    Considering the noise accompanied in the measurements, we rewrite Eq (2.12) as

    b=Sp+n,

    where n is a random vector of mean 0. Let N be the covariance matrix of n. Also, assume that ˜p is normally distributed with mean p and its covariance matrix is P.

    Using multiple measurements and the statistical estimation theory we can find the linear operator L which minimizes the expected difference ErrL between the estimated moments ˜p and the exact moments p.

    ErrL=˜pp2=Lbp2=L(Sp+n)p2=(LSI)p+Ln2=Mp+Ln2(where M=LSI)=Mp2+Ln2(by independence of p and n)=Tr(MPMT)+Tr(LNLT).

    Setting the gradient of ErrL to 0 and solving for L, we get the optimal linear operator

    L=PST(SPST+N)1. (2.15)

    Then, by Eq (2.14) we get the best estimates of the moments.

    Theorem 2.1 (Uniqueness of solutions). Let ui,i=1,2 be the solutions of the problems

    (σui)=mik=1pk(i)δS(i)k in Ω,
    uiν=φ on Γ,

    such that

    u1=u2 on Γ,

    then

    m1=m2=m,
    pk(1)=pk(2),k=1,2,...,m,
    S(1)k=S(2)k,k=1,2,...,m.

    The solution of Poisson equation is the convolution of the fundamental solution of Laplace equation and the source function.

    w(x)=12π[m2k=1pk(xSk)|xS(2)k|2m1k=1pk(xSk)|xS(1)k|2],n=2.
    w(x)=14π[m2k=1pk(xSk)|xS(2)k|3m1k=1pk(xSk)|xS(1)k|3],n=3.

    As EEG imaging data are typically noisy, especially determining the rank of a near singular matrix is very unstable, the error of the numerical reconstruction method needs to be studied. Chafik et al. [12,13] proposed that when the norms of the perturbations (g=˜ff,h=˜φφ) are small in H1/2×H1/2, there exist a>0 and b>0 such that k=1,2,...,m,

    ˜SkSk2m(1Rm)dm1(1R)max{(m1j)Rj,0jm1}(agH1/2(Γ)+bhH1/2(Γ)), (3.1)

    where Sk=xk+iyk is the exact position of the kth dipole, ˜Sk=˜xk+i˜yk is the estimated position of the kth dipole, d is the minimal distance between Sk and ˜Sk, and R1 is a real number bigger than the norm of any point on Γ. However, the analysis is not given by Chafik et al.

    Here we present a new error estimate and provide a proof.

    Theorem 3.1. Suppose m dipoles are enclosed in a circular boundary of radius R. The potential f on the boundary and the gradient of the potential φ perpendicular to the boundary are known. If T is the measurements without noise, and ˜T is the measurements with noise, then the error estimate is given by

    T˜T2m(φ2R2m2πR+f2R2m2πR)(m!mm1pm1maxRm(m1)pmmindm(m1))+2m2(φ2R2m2πR+f2R2m2πR)2(m!mm1pm1maxRm(m1)pmmindm(m1))2, (3.2)

    where p is the moment of dipoles and d is the smallest distance between any two dipoles.

    Proof. We define

    Zi=[βiβi+1βi+m1βi+1βi+2βi+mβi+m1βi+mβi+2m2],iN.

    Then,

    Z1=[β1β2βmβ2β3βm+1βmβm+1β2m1].

    where

    βj=mk=1pkSj1k=mk=1(pk1+ipk2)(xk+iyk)j1,j=1,2,...,2m1.
    det(Z1)=|β1β2βmβ2β3βm+1βmβm+1β2m1|=|pkpkSkpkSm1kpkSkpkS2kpkSmkpkSm1kpkSmkpkS2m2k|=m1m2mmτ(m1,m2,...,mm)pm1pm2pmm|1Sm2Sm1mmSm1S2m2SmmmSm1m1Smm2S2m2mm|
    =m1m2mmτ(m1,m2,...,mm)pm1pm2pmm|111Sm1Sm2SmmSm1m1Sm1m2Sm1mm|S0m1S1m2Sm1mm=p1p2pm|111Sm1Sm2SmmSm1m1Sm1m2Sm1mm|(m1m2mmτ(m1,m2,...,mm)S0m1S1m2Sm1mm)=p1p2pm|111Sm1Sm2SmmSm1m1Sm1m2Sm1mm||111Sm1Sm2SmmSm1m1Sm1m2Sm1mm|=p1p2pm1i<jm(SiSj)2.

    Here, (m1,m2,...,mm) is any permutation of (1,2,...,m) and τ(m1,m2,...,mm) is the sign determined by the permutation.

    The maximum absolute row sum norm is defined by

    A=maxij|aij|,

    where A is a matrix. When A is a vector, A=maxi|ai|.

    In the following proof we will use an important inequality:

    a(x)b(x)a(y)b(y)a(x)a(y)b(x)+b(x)b(y)a(x)

    where a(x) and b(x) can be scalar, vector, or matrix.

    By Cauchy-Schwarz inequality, we have

    R(vj)=φ,vjf,vjν=ΓφvjdsΓfvjνds=Γφ(x+iy)jdsΓf(x+iy)jνds(Γφ2ds)1/2(Γ(x+iy)2jds)1/2+(Γf2ds)1/2(Γ((x+iy)jν)2ds)1/2(Γφ2ds)1/2Rj2πR+(Γf2ds)1/2jRj12πRjφ2Rj2πR+jf2Rj12πR.
    |βj|=|R(vj)j|φ2Rj2πR+f2Rj12πRφ2R2m2πR+f2R2m2πR

    where R>1.

    Let

    T=Z2Z11=Z2adj(Z1)det(Z1)

    where Z1=[β1β2βmβ2β3βm+1βmβm+1β2m1] and Z2=[β2β3βm+1β3β4βm+2βm+1βm+2β2m].

    We can view R(vj) as the measurement obtained by the "detector" vj, while βj is just a constant multiple of R(vj). So, βj is still a measurement of another form, which contains the information about the moment and the position of the dipole source. Since Z1 and Z2 are constructed by different measurements βj, T is also a matrix of measurements.

    Assume T is the measurements without noise, and ˜T is the measurements with noise. Then,

    T˜T=Z2Z11˜Z2˜Z11Z2˜Z2Z11+Z11˜Z11Z2.

    We will analyse the four norms in the above inequality one by one.

    Z2˜Z2mφ˜φ2R2m2πR+mf˜f2R2m2πR.

    To find Z11 we need to estimate adj(Z1). We first observe the results for m=3, then prove the results to the arbitrary m using mathematical induction.

    If Z1=[β1β2β3β2β3β4β3β4β5], then the absolute value of the first element of adj(Z1) would be

    abs(|β3β4β4β5|)=|β3β5β24||β3||β5|+|β24|=(p1S21+p2S22+p3S23)(p1S41+p2S42+p3S43)+(p1S31+p2S32+p3S33)2(3pmaxR2)(3pmaxR4)+(3pmaxR3)2=2(3pmaxR3)2=(31)!331p31maxR3(31)=:max(abs(|β3β4β4β5|)).

    Then,

    adj(Z1)max(abs(|β3β4β4β5|))+max(abs(|β2β4β3β5|))+max(abs(|β2β3β3β4|))3max(abs(|β3β4β4β5|))=3(31)!331p31maxR3(31)=3!331p31maxR3(31).

    Assume when m=n1, we have

    abs(|β3β4βn+1β4β5βn+2βn+1βn+2β2n1|)(n1)!nn1pn1maxRn(n1).

    In fact, this inequality is also true for other minors with matrix size (n1)×(n1).

    Then, when m=n we have

    abs(|β3β4βn+1βn+2β4β5βn+2βn+3βn+1βn+2β2n1β2nβn+2βn+3β2nβ2n+1|)max|β2n+1|max(abs(|β3β4βn+1β4β5βn+2βn+1βn+2β2n1|))++max|βn+2|max(abs(|β4β5βn+2β5β6βn+3βn+2βn+3β2n|))nmax|β2n+1|max(abs(|β3β4βn+1β4β5βn+2βn+1βn+2β2n1|))nmax|p1S2n1+p2S2n2++pnS2nn|(n1)!nn1pn1maxRn(n1)nnpmaxR2n(n1)!nn1pn1maxRn2n)=n!nnpnmaxR(n+1)nn!(n+1)npnmaxR(n+1)n.

    Then, for any m we have

    adj(Z1)=adj([β1β2βm1βmβ2β3βmβm+1βm1βmβ2m3β2m2βmβm+1β2m2β2m1])max(abs(|β3β4βm+1β4β5βm+2βm+1βm+2β2m1|))++max(abs(|β2β3βmβ3β4βm+1βmβm+1β2m2|))mmax(abs(|β3β4βm+1β4β5βm+2βm+1βm+2β2m1|))=m(m1)!mm1pm1maxRm(m1)=m!mm1pm1maxRm(m1).

    Thus,

    Z11=adj(Z1)det(Z1)m!mm1pm1maxRm(m1)p1p2pm1i<jm(SiSj)2m!mm1pm1maxRm(m1)pmmindm(m1)

    where d is the smallest distance between any two dipoles.

    Notice that

    Z1(Z11˜Z11)+(Z1˜Z1)˜Z11=0.
    Z11˜Z11=Z11(Z1˜Z1)˜Z11.
    Z11˜Z11Z11Z1˜Z1˜Z11(m!mm1pm1maxRm(m1)pmmindm(m1))2(mφ˜φ2R2m2πR+mf˜f2R2m2πR).

    Based on the above results, we have

    T˜TZ2˜Z2Z11+Z11˜Z11Z2(mφ˜φ2R2m2πR+mf˜f2R2m2πR)(m!mm1pm1maxRm(m1)pmmindm(m1))+(m!mm1pm1maxRm(m1)pmmindm(m1))2(mφ2R2m2πR+mf2R2m2πR)(mφ˜φ2R2m2πR+mf˜f2R2m2πR)2m(φ2R2m2πR+f2R2m2πR)(m!mm1pm1maxRm(m1)pmmindm(m1))+2m2(φ2R2m2πR+f2R2m2πR)2(m!mm1pm1maxRm(m1)pmmindm(m1))2. (3.3)

    We can further simplify it as

    T˜TE+E2 (3.4)

    where

    E=2mR2m2πR(f2+φ2)(m!mm1pm1maxRm(m1)pmmindm(m1)). (3.5)

    When 0<E<1, the error in the position estimate is mainly controlled by E; when E>1, the error in the position estimate is mainly controlled by E2.

    Let Ω be a circular disk centered at the origin and of radius r=1. Then, the numerical implementation can be simplified as follow.

    vjν=(x+iy)jr=(reiθ)jr=jrj1eiθj=jrjeiθjr=jvjr.
    R(vj)=f,vjν=ΓfvjνdΓ=2π0fjvjrrdθ=j2π0fvjdθ=j2π0f(reiθ)jdθ,

    where f is a function of θ on the boundary. We do not need to know the explicit form of f, but we can measure as many points as possible on the boundary to get enough discretized function values of f. Then, the above integral can be approximated by a Riemann sum.

    The measurable values we want to use in the following are

    βj=R(vj)j=2π0f(reiθ)jdθ.

    The Romberg algorithm is used to calculate the integral numerically.

    We compare the efficacy of the harmonic function method in dipolar source reconstruction when the perturbation level is 0,0.001,0.01,0.1 and the number of dipoles is 1,2,3,4,5. It is shown that as the perturbation level increases, the reconstruction error increases (see Figures 15).

    Figure 1.  The effect of the perturbation level on the reconstruction error of 1 dipole. As the perturbation level increases, the reconstruction error increases. Here, the perturbation means adding noise to the exact measurement. If the perturbation level is σ, then the perturbed measurement is the exact measurement times (1±σ), where plus or minus signs are randomly assigned to each channel. Also, the error is defined as the sum of position errors.
    Figure 2.  The effect of the perturbation level on the reconstruction error of 2 dipoles. As the perturbation level increases, the reconstruction error increases. Here, the perturbation means adding noise to the exact measurement. If the perturbation level is σ, then the perturbed measurement is the exact measurement times (1±σ), where plus or minus signs are randomly assigned to each channel. Also, the error is defined as the sum of position errors.
    Figure 3.  The effect of the perturbation level on the reconstruction error of 3 dipoles. As the perturbation level increases, the reconstruction error increases. Here, the perturbation means adding noise to the exact measurement. If the perturbation level is σ, then the perturbed measurement is the exact measurement times (1±σ), where plus or minus signs are randomly assigned to each channel. Also, the error is defined as the sum of position errors.
    Figure 4.  The effect of the perturbation level on the reconstruction error of 4 dipoles. As the perturbation level increases, the reconstruction error increases. Here, the perturbation means adding noise to the exact measurement. If the perturbation level is σ, then the perturbed measurement is the exact measurement times (1±σ), where plus or minus signs are randomly assigned to each channel. Also, the error is defined as the sum of position errors.
    Figure 5.  The effect of the perturbation level on the reconstruction error of 5 dipoles. As the perturbation level increases, the reconstruction error increases. Here, the perturbation means adding noise to the exact measurement. If the perturbation level is σ, then the perturbed measurement is the exact measurement times (1±σ), where plus or minus signs are randomly assigned to each channel. Also, the error is defined as the sum of position errors.

    In the following we show the results of source estimation, assuming there are 3 dipolar sources (m=3).

    ● Dipole 1: position (0.3,0.3) and moment (0,1).

    ● Dipole 2: position (0.6,0.2) and moment (1,1).

    ● Dipole 3: position (0.5,0.4) and moment (2,2).

    In the graphs (see Figure 6) we use a small circle and a red line segment to indicate the true value, and use a cross sign and a green line segment to indicate the reconstructed values.

    Figure 6.  The effect of the perturbation level on the reconstruction error of 3 dipoles. As the perturbation level increases, the reconstruction error increases.

    From error estimates we know that as the distance between two dipoles gets closer, the reconstruction error for the positions of dipoles gets larger (see Table 1 and Figure 7). This is verified by the numerical simulations.

    Table 1.  The effect of dipole distance on the reconstruction error. As two dipoles get closer, the mean reconstruction error in the positions of the dipoles gets larger, which is consistent with the result in the error estimate.
    Exact Dipole Distance Reconstructed Dipole Distance
    0.03 0.7429
    0.05 0.3084
    0.10 0.1200

     | Show Table
    DownLoad: CSV
    Figure 7.  The effect of dipole distance on the reconstruction error. As two dipoles get closer, the reconstruction error in the positions of the dipoles gets larger, which is consistent with the theoretical analysis in the error estimate. When dexact=0.10, ¯dest=0.1200; when dexact=0.05, ¯dest=0.3084; when dexact=0.03, ¯dest=0.7429.

    We randomly assign two dipoles with fixed distance, say 0.1, in the unit disk, then reconstruct their positions. We fix the noise level for all experiments at σ=0.001.

    Let di (i=1,2) be the distance between the ith exact dipole and the ith estimated dipole, and dmax be the largest d.

    We repeat the experiment 10 times and show their performance on average over different dipole distances.

    The above experiment also provides a numerical example to show that the estimate given by us in Theorem 3.1 provides a better error bound when the two poles are very close.

    When the number of dipoles is m=2, Chafik's estimate is bounded by C1d (see Inequality (3.1)), while our estimate is bounded by C2d2 (see Inequality (3.4) and Eq (3.5)) where d is the smallest distance between two dipoles and Ci (i=1,2) are constants independent of d. That is, when the distance is halved, the error bound will be amplified by 2 in Chafik's estimate and by 4 in our estimate.

    From the data simulation, we see that

    0.050.03=1.67<0.74290.3084=2.41<1.672=2.79.
    0.100.05=2<0.30840.1200=2.57<22=4.
    0.100.03=3.33<0.74290.1200=6.19<3.332=11.09.

    For example, when the distance between the two dipoles is reduced from 0.10 to 0.05, by Chafik's estimate the error should be amplified by 2, but in fact, the error is amplified by 2.57, which is bounded by 4 in our estimate.

    In this paper we studied a harmonic function method for the dipolar source reconstruction, derived error estimate for the harmonic function method and compared our result with Chafik's estimate. By numerical simulations it is shown that the harmonic function method can quickly and accurately locate active regions in EEG source reconstruction. In the future, we plan to extend the harmonic function method to 3D case and applied this method to some real EEG data. The brain's conductance variation in different brain regions also leads to additional challenges in source localization [35]. Although these tissue properties can be quantified through MRI methods, numerical methods such as finite element method will be needed to solve the inverse problems. Since the estimation of the number of dipoles relies on the calculation of the rank of the measurement matrix, which is significantly affected by the noise, we hope to find some way to solve or circumvent this problem. In addition, the situation that the number of exact dipoles is not equal to the estimated value could also be considered. Furthermore, when two dipoles get close enough, it may be better to regard them as an equivalent dipole to avoid increased error.

    We thank anonymous reviewers for providing us with valuable suggestions.

    All authors declare no conflicts of interest in this paper.



    [1] N. K. Logothetis, What we can do and what we cannot do with fMRI, Nature, 453 (2008), 869–878. https://doi.org/10.1038/nature06976 doi: 10.1038/nature06976
    [2] B. He, J. Lian, High-resolution spatio-temporal functional neuroimaging of brain activity, Crit. Rev. Biomed. Eng., 30 (2002), 283–306. https://doi.org/10.1615/CritRevBiomedEng.v30.i456.30 doi: 10.1615/CritRevBiomedEng.v30.i456.30
    [3] P. L. Nunez, Electric Fields of the Brain: The Neurophysics of EEG, Oxford University Press, 1981.
    [4] R. Grave de Peralta Menendez, S. L. Gonzalez Andino, S. Morand, C. M. Michel, T. Landis, Imaging the electrical activity of the brain: ELECTRA, Hum. Brain Mapp., 9 (2000), 1–12. https://doi.org/10.1002/(SICI)10970193(2000)9:1<1::AIDHBM1>3.0.CO;2# doi: 10.1002/(SICI)10970193(2000)9:1<1::AIDHBM1>3.0.CO;2#
    [5] C. M. Michel, M. M. Murray, G. Lantz, S. Gonzalez, L. Spinelli, R. Grave de Peralta, EEG source imaging, Clin. Neurophysiol., 115 (2004), 2195–2222. https://doi.org/10.1016/j.clinph.2004.06.001
    [6] A. Hunter, B. Crouch, N. Webster, B. Platt, Delirium screening in the intensive care unit using emerging qeeg techniques: A pilot study, AIMS Neurosci., 7 (2020), 1–16. https://doi.org/10.3934/Neuroscience.2020001 doi: 10.3934/Neuroscience.2020001
    [7] G. V. Portnova, O. Ivanova, E. V. Proskurnina, Effects of eeg examination and aba-therapy on resting-state eeg in children with low-functioning autism, AIMS Neurosci., 7 (2020), 153–167. https://doi.org/10.3934/Neuroscience.2020011 doi: 10.3934/Neuroscience.2020011
    [8] R. Grech, T. Cassar, J. Muscat, K. P. Camilleri, S. G. Fabri, M. Zervakis, et al., Review on solving the inverse problem in EEG source analysis, J. NeuroEng. Rehabilitation, 5 (2008), 25. https://doi.org/10.1186/1743-0003-5-25 doi: 10.1186/1743-0003-5-25
    [9] V. Isakov, Inverse Problems for Partial Differential Equations, Springer, 2nd edition, 2005.
    [10] S. Vessella, Locations and strengths of point sources: stability estimates, Inverse Probl., 8 (1992), 911. https://doi.org/10.1088/0266-5611/8/6/008 doi: 10.1088/0266-5611/8/6/008
    [11] R. M. Leahy, J. C. Mosher, M. E. Spencer, M. X. Huang, J. D. Lewine, A study of dipole localization accuracy for meg and eeg using a human skull phantom, Electroencephalogr. clin. neurophysiol., 107 (1998), 159–173. https://doi.org/10.1016/S0013-4694(98)00057-1 doi: 10.1016/S0013-4694(98)00057-1
    [12] A. El Badia, T. Ha-Duong, An inverse source problem in potential analysis, Inverse Probl., 16 (2000), 651–663. https://doi.org/10.1088/0266-5611/16/3/308 doi: 10.1088/0266-5611/16/3/308
    [13] M. Chafik, A. El Badia, T. Ha-Duong, On some inverse eeg problems, Inverse Probl. Eng. Mech. II, pages 537–544, 2000. https://doi.org/10.1016/B978-008043693-7/50129-X
    [14] T. Nara, S. Ando, A projective method for an inverse source problem of the poisson equation, Inverse Probl., 19 (2003), 355–369. https://doi.org/10.1088/0266-5611/19/2/307 doi: 10.1088/0266-5611/19/2/307
    [15] H. Kang, H. Lee, Identification of simple poles via boundary measurements and an application of eit, Inverse Probl., 20 (2004), 1853. https://doi.org/10.1088/0266-5611/20/6/010 doi: 10.1088/0266-5611/20/6/010
    [16] A. El Badia, Inverse source problem in an anisotropic medium by boundary measurements, Inverse Probl., 21 (2005), 1487. https://doi.org/10.1088/0266-5611/21/5/001 doi: 10.1088/0266-5611/21/5/001
    [17] L. Baratchart, A. Ben Abda, F. Ben Hassen, J. Leblond, Recovery of pointwise sources or small inclusions in 2d domains and rational approximation, Inverse Probl., 21 (2005), 51. https://doi.org/10.1088/0266-5611/21/1/005 doi: 10.1088/0266-5611/21/1/005
    [18] Y. S. Chung, S. Y. Chung, Identification of the combination of monopolar and dipolar sources for elliptic equations, Inverse Probl., 25 (2009), 085006. https://doi.org/10.1088/0266-5611/25/8/085006 doi: 10.1088/0266-5611/25/8/085006
    [19] D. Kandaswamy, T. Blu, D. van de Ville, Analytic sensing: Noniterative retrieval of point sources from boundary measurements, SIAM J. Sci. Comput., 31 (2009), 3179–3194. https://doi.org/10.1137/080712829
    [20] T. Nara, S. Ando, Direct localization of poles of a meromorphic function from measurements on an incomplete boundary, Inverse Probl., 26 (2010), 015011. https://doi.org/10.1088/0266-5611/26/1/015011 doi: 10.1088/0266-5611/26/1/015011
    [21] A. El Badia, T. Nara, An inverse source problem for helmholtz's equation from the cauchy data with a single wave number, Inverse Probl., 27 (2011), 105001. https://doi.org/10.1088/0266-5611/27/10/105001 doi: 10.1088/0266-5611/27/10/105001
    [22] M. Clerc, J. Leblond, J.-P. Marmorat, T. Papadopoulo, Source localization using rational approximation on plane sections, Inverse Probl., 28 (2012), 055018. https://doi.org/10.1088/0266-5611/28/5/055018 doi: 10.1088/0266-5611/28/5/055018
    [23] R. Mdimagh, I.Ben Saad, Stability estimates for point sources identification problem using reciprocity gap concept via the helmholtz equation, Appl. Math. Model., 40 (2016), 7844–7861. https://doi.org/10.1016/j.apm.2016.03.024 doi: 10.1016/j.apm.2016.03.024
    [24] J. Vorwerk, Ü. Aydin, C. H. Wolters, C. R. Butson, Influence of head tissue conductivity uncertainties on eeg dipole reconstruction, Front Neurosci., 13 (2019), PMCID: PMC6558618. https://doi.org/10.3389/fnins.2019.00531
    [25] M. Rubega, M. Carboni, M. Seeber, D. Pascucci, S. Tourbier, G. Toscano, et al., Estimating eeg source dipole orientation based on singular-value decomposition for connectivity analysis, Brain topogr., 32 (2019), 704–719. https://doi.org/10.1007/s10548-018-0691-2 doi: 10.1007/s10548-018-0691-2
    [26] M. S. Hämäläinen, R. J. Ilmoniemi, Interpreting magnetic fields of the brain: minimum norm estimates, Med. Biol. Eng. & Comput., 32 (1994), 35–42. https://doi.org/10.1007/BF02512476
    [27] R. D. Pascual-Marqui, Review of methods for solving the eeg inverse problem, Int. J. Bioelectromagn., 1 (1999), 75–86, 1999.
    [28] S. Baillet, Toward functional brain imaging of cortical electrophysiology markovian models for magneto and electroencephalogram source estimation and experimental assessments, Orsay, France, 1998.
    [29] S. Baillet, J. C. Mosher, R. M. Leahy, Electromagnetic brain mapping, IEEE Signal Process. Mag., 18 (2001), 14–30. https://doi.org/10.1109/79.962275
    [30] J. C. Mosher, P. S. Lewis, R. M. Leahy, Multiple dipole modeling and localization from spatio-temporal MEG data, IEEE Trans. Biomed. Eng., 39 (1992), 541–557. https://doi.org/10.1109/10.141192
    [31] P. A. Muñoz Gutiérrez, E. Giraldo, M. Bueno-López, M. Molinas, Localization of active brain sources from eeg signals using empirical mode decomposition: a comparative study, Front. Integrat. Neurosci., 2 (2018), 55. https://doi.org/10.3389/fnint.2018.00055 doi: 10.3389/fnint.2018.00055
    [32] C. Kaur, P. Singh, S. Sahni, Electroencephalography-based source localization for depression using standardized low resolution brain electromagnetic tomography-variational mode decomposition technique, Eur. Neurol., 81 (2019), 63–75. https://doi.org/10.1159/000500414 doi: 10.1159/000500414
    [33] V. P. Oikonomou, I. Kompatsiaris, A novel bayesian approach for eeg source localization, Comput. Intell. Neurosci., 2020. https://doi.org/10.1155/2020/8837954
    [34] S. Andrieux, A. Ben Abda, H.D. Bui, Reciprocity principle and crack identification, Inverse Probl., 15 (1999), 59. https://doi.org/10.1088/0266-5611/15/1/010 doi: 10.1088/0266-5611/15/1/010
    [35] H. Azizollahi, M. Darbas, M. M. Diallo, A. El Badia, S. Lohrengel, EEG in neonates: Forward modeling and sensitivity analysis with respect to variations of the conductivity, Math. Biosci. Eng., 15 (2018), 905–932. https://doi.org/10.3934/mbe.2018041 doi: 10.3934/mbe.2018041
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1924) PDF downloads(52) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog