Studies On The Numerical Methods Of The Inverse Problems In EEG And MEG

EEG and MEG, which involve the most advanced technologies in the 21th century including the areas of biotechnology, electronic engineering, medical engineering and so on, are the latest noninvasive techniques for investigating neuronal activities of the living human brain. At present, EEG and MEG have attracted the attentions of people among the world. This paper just considers the inverse source problem arising from EEG/MEG. We built the mathematical models of the inverse EEG/MEG problems for the static and harmonic electromagnetic fields respectively, and then pay our main attentions to the numerical reconstruction methods for the pointwise sources of which the uniqueness results have been proved by the predecessors. For the static field, we propose two methods-the algebraic direct method and the rational interpolation method. For the harmonic field, we just consider the dipoles and give out a reconstruction method by algebraic considerations.I. Reconstruction of pointwise sources in static fieldsSupposing the simply connected domainΩ(?)Rd(d= 2,3) is bounded, and its bound-aryΓbelongs to C1 class. For the single layer model of EEG/MEG in static field, the potential u satisfies the Cauchy problem Here, F is the source term. If there exist monopoles with strengthλj∈R at points Sj∈Ωfor j = 1,2,...,N1 and dipoles with strength Pk∈Rd at points Ck∈Ωfor k= 1,2,...,Nz, then F is of the form: Here, be the total number of the pointwise sources. Assume that the upper bound of N is known to be M and there existδ> 0, d> 0 such thatAt this time, the inverse EEG/MEG problem of single layer model of static field is:(IP1):Given the Cauchy boundary data g∈H1/2(Γ) and (?)∈H-1/2(F), reconstruct the source term F of the form (3), i.e,reconstruct the number parameters N1 and N2, the location parameters Sj and Ck and the strength parametersλj and Pk.In [15], the uniqueness of the solution of problem (IP1) has been proved. Based on this, we discuss the numerical methods for the static model in this paper, where two methods-algebraic direct method and rational interpolation method are considered.Ⅰ.1. Algebraic direct method Supposing the sequence{an} is of the form We assume that there existδ0> 0, do> 0 such that pk, ck∈C satisfyingTheorem 1 ([50]) Suppose{an} is defined as (5). Let then In particular, Theorem 2 ([50]) Suppose{an} is defined as (5). If l1,...,lN satisfy the equations then c1,c2,..., cN are solutions of The converse is also true.Supposing M is the upper bound of N, Theorem 1 and Theorem 2 implicate a numerical method to reconstruct the parameters N and pk, ck by the values of{an}n=02M-1:1. For n= M, M-1,..., calculate the determinant of Dn in turn. The first integer such that det(Dn)≠0 is identified as N.2. After N has been determined, solve the equations (9) to obtain the values of{lk}k=1N. Substitute them to equation (10), then{ck}k=1N are obtained.3. Using the values of{ck}k=1N obtained by the above step, solve the equations to find p1,p2,...,PN.Given the algebraic theory, we discuss the application of the algebraic method in the inverse EEG/MEG problem.(一). Reconstruction of the non-combined sources in a discConfine that the solving domainΩis the unit disc in R2 with its center at the origin, and denote it by B.Identifying R2 with the complex plane Z, let x= (x1,x2), Sj= (Sj,1,Sj,2), Ck= (Ck,1,Ck,2),Pk＝(Pk,1,Pk,2)andTheorem 3 ([67]) Suppose that u is the solution of the Cauchy problem (1) and (2). Then there exist a function A analytic in B such that u(z)= Re(f(z)), where Here, Re(·) represents the real part of a complex number. Moreover, if fr(?)ds= 0, then Here, for every givenξ0∈Γ,C= i(?)(ξ0) is a constant.From the Cauchy integral formula and Theorem 3, two conclusions are obtained:Theorem 4 Supposing there just exist dipoles in B, i.e. N1= 0. For non-negative integer n, define thenTheorem 5 Suppose that there just exist monopoles in B, i.e, N2= 0. For non-negative integer n, define thenAs a result of Theorem 4 and Theorem 5, we obtain a numerical method to reconstruct the dipoles or monopoles for the unit disc domain:·Supposing g and (?) are given,for n= 0,1,2,...,2M - 1, calculate the values of an orβn by the formulas (12) or (14).·According to the formulas (13) or (15), N2, Ck, Pk or N1, Sj,λj can be reconstruct by the algebraic method.(二). Reconstruction of the non-combined sources in smooth domainsSupposing thatΩ(?)Rd(d= 2,3) is bounded and its boundaryΓis smooth. Based on the variational form of the boundary problem, we construct the algebraic direct method for the non-combined sources for the R2 and R3 cases respectively. Choose the test function v∈H(Ω)={v∈H1(Ω),△v=0),then the variational form of the Cauchy problem(1)and(2)is where,1.The monopole caseSuppose that there just exist monopoles,i.e,N=N1.Case 1.1(The R2 case). For m=0,1,...,2M一1,let vm(1)=Re(zm),vm(2)= Im(zm),and vm:=vm(1)+ivm(2)=zm.Substituting vm to(17),thenThe reconstruction method for the monopoles in smooth domain of R2 case is:.For m=0,1,...,2M一1,choosing vm=zm,compute the values of R(vm)by(18)..Lettingαm=R(vm),reconstruct the parameters N1,Sj,λj by the algebraic method.Case 1.2(The R3 case). Suppose that Sj=(Sj,1,Sj,2,Sj,3)and sj=Sj,1+iSj,2 for j=1,2,...,N1.Choosing vm=(x1+ix2)m and substituting them to(17),we obtain Assuming that there are not two points in {Sj)j=1N1 sharing the same projection the plane xlx2,which means that sj satisfy the condition(6).For sequence(19),adopt the algebraic method,then the values of N1,{Sj,1,Sj,2)j=1N1,{λj}j=1N1 can be obtained.Choosing the other test functioin (?)m=x3(x1+ix2)m,we have Substituting the values of N1,{Sj,1,Sj,2}j=1N1,{λj}j=1N1 to the equations(20),then{Sj,3)j=1N1 are obtained by solving this linear equations.Hence,the reconstruction method for the monopoles in smooth domain of R3 is: ·For m=0,1,...,2M-1,letting vm=(x1+ix2)m,compute the values of R(vm) defined by(17).Then the parameters can be obtained by the algebraic method..For m=1,2,...,N1,letting (?)m=x3(x1+ix2)m,compute the values of R((?)m).Solve (20)to find{Sj,3}j=1N1.2.The dipole caseIf there just exist dipoles,then N=N2.Case 2.1(The R2 case).For m=0,1,...,2M-1,choose vm=zm.We have where pk=Pk,1+iPk,2,ck=Ck,1+iCk,2.Therefore,the reconstruction method for the dipoles in smooth domain of R2 is:.For m=0,1,...,2M-1,choosing vm=zm,calculate the values of R(vm)..Letting reconstruct the parameters by the algebraic method.Case 2.2(The R3 case).For j=1,2,...,N2,let Pk=(Pk,1,Pk,2,Pk,3),Ck= (Ck,1,Ck,2,Ck,3)and pk=Pk,1+iPk,2,ck=Ck,1+iCk,2.Assume that there not exist two dipoles such that the x1x2 projections of their location parameters are the same.Letting vm=(x1+ix2)m,then And then choosing (?)m=x3(x1+ix2)m,we have.For m=0,1,2,…,2M一1,letting vm=(x1+ix2)m,αm. Using the algebraic method,reconstruct the parameters and ·For m= 0,1,...,2N2-1, letting calculateβm. Under the conditions that the parameters have been obtained, solve the algebraic equations (21) to find andFor the purpose of verifying the methods we proposed, we give out several simple examples. For the unit disc domain, comparing the numerical results of the two methods, which the one base on the complex expression of the solution and the other base on the variational form of the boundary problem, we propose that the accuracy of the them are accordant. Moreover, we also show the validity of the variational algebraic method for the other smooth domains.1.2. Rational interpolation methodFor the reconstruction of the combined sources, we construct an effective rational interpolation method. Introducing the inverse difference and the continued fraction, we give a method to determine the degree of the denominator of a rational function, and then invert the combined sources by rational interpolation.1. The R2 caseWe denote by L(P(z)) and d(P(z)) the leading coefficient and the degree of the polynomial function P(z). Let (?)n to be the set of polynomials with degree no more than n. Define the rational function setTheorem 6 Assume that and the preconditions (?)(Q)>(?)(P) and (?)(Q)=n(1≤n≤n0) are known. Given the notes zo,z1,..., z2n-1 which are different from each other such thatκ[z0, z1,…,zl]≠0,∞for l= 0,1,…,2n-1, we conclude that the inverse difference functions k[z0, z1,…,zk,z] are rational functions and belong to & Rn,n for k= 0,1,...,2n-1.To be precise, letting thenFrom the Theorem 6, we conclude that:For even number k, supposing the well posed notes z0, z1,...,zk are given, find the other n0-k/2 points satisfying For k=0,2,4,...,we identity the first number such that is the degree of H(z).IdentifyingΩ(?)R2 as a domain in C.For fixed point z∈C／(?),letting v(x,z)= 1/(z-x),then Letting R(z)=G(z)/H(z),it is obvious that and G(z) and H(z)are relatively-prime.Letting satisfies the preconditions of Theorem 6.From the results of Theorem 6,we can determine the degree of H(z).If the degree of denominator H(z)is known to be N0,given(zk,R(zk))for k= 0,1,...,2No such that k[z0,z1,…,zl]≠0,∞for any l=0,1,...,2N0,then From this,we obtain a rational interpolation method to invert the rational function R(z) by finite discrete values.Here,{sj}j=1N1 and{ck}k=1N2 are the poles of R(z)with first order and second order respectively.2.The R3 caseFor fixed point letting thenIf we suppose that the projections of {Sj}j=1N1 and {Ck}k=1N2 on the x1x2 plane are dif-ferent,the rational function R(ξ)can be inverted using the rational interpolation method which is proposed for the two dimensional case.At this time,we can obtain the parame-ters andChoosing the other test functions we have Suppose that the numerator and denominator of and respectively. From the numerical results of and we can write out directly. Adopting the Pade-type approximation, rational function is obtained. Then compute the residues of at the poles and to obtain Compute the residues of at the poles to obtainThe numerical examples illustrate that the rational interpolation method is applicable to the combined sources for static fields. Comparing to the algebraic method, the stability of rational interpolation method has descended. We point out that we just offer a new idea to reconstruct the combined sources, some other problems such as improving the accuracy of the method or adopting other interpolation method with better stability are waiting for further research.II. Reconstruction of dipoles in harmonic fieldsSuppose that the simply connected domain included the origin is bounded and its boundaryΓbelongs to C1 class. For the harmonic electromagnetic field, the electric field strength E and the magnetic field strength H satisfy where w,ε,μ,σare constants inΩ. If there are finite dipoles inΩ, the source function JP is of the form: where JP∈R3 is a vector function.The inverse source problem of harmonic field we discussed is:(IP2):Given the boundary conditions n×E|r and nx H|Γ, reconstruct the source term JP, i.e, reconstruct the number N, the locations{Cm}mN=1 and the strength{Pm}mN=1.Based on the uniqueness result of (IP2) proposed in [85], we consider the numerical reconstruction method for dipoles. Using the variational form of the problem, we choose suitable test functions to deduce the algebraic relations among the parameters and then solve the problem by algebraic consideration.Introducing the operator where define For e∈N(M,Ω), the variational form of (25)-(27) is whereIn the following discussion we would see that the parameters of dipoles can be re-constructed by the values of R(e). For the convenience of the choice of e, we give a characterization of a subset of N(M,Ω).Theorem 7 If e∈(C2,α((?)))3 satisfies then M(e)=0,i.e, e∈M(M,Ω).(1). Reconstruction of one dipoleAt first, we consider the case that there just exist one dipole in the domain, i.e, we know that N=1 andLetting e1= (0, eikx1,ieikx1) and substituting it to (29), we obtainCase 1.1. R(e1)≠0. Choosing then Suppose that Re(·), Im(·) and arg(·) represent the real part, image part and and amplitude of a complex number respectively. Letη=e-ikCo,1, then we have It is obvious thatRemark 1 In the actual MEG model, the value of k is close to zero and satisfies In view of the stability of the algorithm, we will do not useηto solve C0.1 directly but just take it as a middle value to invert the other parameters.Case 1.2. R(e1)= 0.Choosing thenIf R(e1)=0, which means that P0,2=P0,3=0, then F0,1≠0. Therefore, R(e2)≠0, R(e4)≠0. We have (2). Reconstruction of dipolesIf there are finite dipoles in the homogeneous medium, then Here, the upper bound of the number of the dipoles is known to be M. Moreover, we suppose that there exist constantsγ0> 0, do> 0 such thatStep 1. Reconstruct the parameters N and{Cm,2, Cm,3}mN=1.At first, choose test functions Substituting en to (29), we have where cm andβm satisfy the condition (6), then the parameters N, cm andβm can be obtained by the algebraic method andStep 2. Reconstruct parametersIf the number of the dipoles N has been obtained according to Step 1, letting then we have where Solving the linear algebraic equations (33),βm(1) can be uniquely determined. Letting thenStep 3.Reconstruct parameters For n=0,1,...,2N-1,choosing then where Identifyingβm(2) as the unknown variables and supposing the parameter matrix of(34)is A,we obtain th Here,[N/2]represents the integral part of N/2.case 2.1.det A≠0.It means that cm≠0 for m=1,2,...,N.We haveCase 2.2.det A=0.It means that there exist an cm0=0(1≤m0≤N).We suppose that m0=N,i.e,CN=CN,2+iCN,3=0.For the equations suppose its parameter matrix is B,then det B≠0 which means that the solution of(35) is unique.Hence Choosing e(3)= (eikx2,0,0), thenFrom the numerical results, we conclude that the algebraic method for the recon-struction of dipoles is valid. This algebraic method, which need neither solving forward problem nor any iterative step, is a simple and practicable algorithm. For the improving of the stability of the algebraic method, the regularization strategy can be considered, which is waiting for further research.