JP3237359B2 - Life activity current source estimation method - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a method for estimating the position, orientation, and size of a biological activity current source, and more particularly to a method for estimating a biological activity current source using a minimum norm method.

[0002]

2. Description of the Related Art When a living body is stimulated, the polarization formed across a cell membrane is broken and a living activity current flows. This biological activity current appears in the brain and heart, and is recorded as an electroencephalogram and an electrocardiogram. The magnetic field generated by the biological activity current is recorded as a magnetoencephalogram and a magnetocardiogram.

Recently, as a device for measuring a minute magnetic field in a living body, SQUID (Superconducting Quantum Inter
face Device: A sensor using a superconducting quantum interferometer) has been developed. This sensor is placed outside the head, and the sensor can non-invasively measure a small magnetic field generated by a current dipole (hereinafter simply referred to as a current source), which is a biological activity current source, generated in the brain. From the measured magnetic field data, the position, direction, and size of the current source related to the lesion are estimated, and the estimated current source is displayed on a tomographic image obtained by an X-ray CT or MRI device, and the physics of the affected part is displayed. It is used to identify the target position.

Conventionally, as one of current source estimation methods, there is a method using a minimum norm method (for example, WHKullman
n, KDJandt, K. Rehm, HASchlitt, WJDallas and
WESmith, Advances in Biomagnetism, pp.571-574, P
lenum Pless, New York, 1989).

A conventional current source estimation method using the minimum norm method will be described below with reference to FIG. As shown in FIG. 7, a multi-channel SQUID sensor 1 is provided near the subject M. Multichannel SQUID sensor 1 is accommodated by immersing a number of magnetic sensors (pickup coil) S ₁ to S _m in the refrigerant such as liquid nitrogen in a container called a dewar.

On the other hand, a large number of grid points (1) to (n) are set three-dimensionally in, for example, the brain, which is a diagnosis target area of the subject M, and an unknown current source (current dipole) is set at each grid point. Assuming that each current source is represented by a three-dimensional vector VP _j (j = 1 to n). Then, the magnetic field B _i .about.B _m detected by the magnetic sensor S ₁ to S _m of the SQUID sensor 1 is expressed by the following equation (1).

[0007]

(Equation 1)

In equation (1), VP _j = (P _jx , P _jy , P _jz ) α _ij = (α _ijx, α _ijy, α _ijz ). Note that α _ij is the magnetic sensors S ₁ to S ₁ when a current source having a unit size in the X, Y, and Z directions is placed on a grid point.
Is a known coefficient which represents the strength of the magnetic field detected by each position of the S _m.

[B] = (B ₁ , B ₂ ,..., B _m ) [P] = (P _1x , P _1y , P _1z , P _2x , P _2y , P _{2z ,.}
(P _nx , P _ny , P _nz ), the expression (1) can be rewritten into a linear relational expression like the expression (2). [B] = A [P] (2)

In the equation (2), A is a matrix having 3n × m elements represented by the following equation (3).

[0011]

(Equation 2)

Here, when the inverse matrix of A is represented by A ⁻ ,
[P] is represented by the following equation (4). (P) = A ^- (B) ......... (4)

The simultaneous equation represented by the equation (4) is obtained by comparing the number m of equations (the number of magnetic sensors S _{1 to} S _m , for example, several 10 to several hundred) with the number n of unknowns (for each grid point). Since the number of assumed current sources is, for example, several hundred to several thousand), no solution is obtained. Therefore, a condition is added to minimize the norm | [P] | of the vector [P]. Then, the above equation (4) is expressed as the following equation (5). [P] = A ⁺ [B] (5) where A ⁺ is a generalized inverse matrix expressed by the following equation (6). ^{A + = A t (AA t} ) -1 ......... (6) However, A ^t is the transpose matrix of A.

By solving the above equation (5), the current sources VP _j on each grid point
Are estimated, and the one having the largest value among them is assumed to be close to the true current source. This is the principle of the current source estimation method using the minimum norm method.

However, in the above-mentioned minimum norm method, since current sources are assumed at all positions where the grid points are arranged, the degree of freedom of the solution increases spatially, which may cause an erroneous solution to be obtained. I have.

Therefore, based on the assumption that the magnetic field source is composed of current sources existing in several local regions, grid points having a low possibility of having a solution are gathered near higher grid points to obtain a space. The present inventors have proposed a method of obtaining a solution that is closer to the true solution by lowering the degree of freedom of the solution by adding a general constraint (hereinafter referred to as a lattice point moving minimum norm method) (Tomita) Sada, Shigeki Kajiwara, Yasushi Kondo, Keiichi Yoshida, Shinji Yamamoto, Takashi Otsu, Naoichi Yakimaki, Hisashi Kato: Transactions of the 8th Annual Meeting of the Japanese Society of Biomagnetism, Vol. 6, pp. 86-89, 1993
). Hereinafter, the grid point moving minimum norm method will be described.

In the grid point moving minimum norm method, the probability that a solution exists on the j-th grid point is expressed by the following equation (7).

[0018]

(Equation 3)

In the above equation (7), Vr _j and VP _j are the position vector of the j-th grid point and the minimum norm solution on the grid point, and n is the number of grid points. The first term of equation (7) represents the intensity of the current source on the grid point, and the second term represents the density of the current source near the grid point. Further, β and γ are empirical parameters, and are appropriately determined based on, for example, an isomagnetic field map created from data obtained by the above-described magnetic sensors. The second term of equation (7) is that if the individual current sources have small intensities but a large number are dense, those lattice points will be replaced by other lattice points with low density but strong individual intensities. Added to avoid subordination. formula
The higher the accuracy given by (7), the higher the possibility that a solution exists near that grid point.

Next, another grid point is moved to a position near a grid point having a large value of the equation (7). However, there is a possibility that the current source is located at a plurality of positions or is spread spatially. There is also. So, we divide the grid points into several groups,
Move the grid points within each group. In order to group each grid point, a group function represented by the following equation (8) is used.

[0021]

(Equation 4)

The above equation (8) represents the degree of influence of the j-th grid point on the i-th grid point. α is a parameter that determines the shape of the function of equation (8), and is appropriately set empirically, for example, similarly to β and γ described above.

FIG. 8 shows an example in which grid points are grouped. In the figure, the horizontal axis is a position vector of a grid point, the vertical axis is a group function, A, B, and C are grid points, φ _A , φ
_B and φ _C are group functions of the lattice points A, B and C.
When the lattice point B has the largest influence on the lattice point A, the lattice point A is assumed to be subordinate to the lattice point B. Therefore, grid points A and B belong to the same group. on the other hand,
Since the lattice point C itself has the greatest influence on the lattice point C, the lattice point C belongs to a different group from the lattice points A and B. In this way, the grid point group is divided into a plurality of groups.

After the grid points are divided into groups, the other grid points belonging to the same group are moved to the grid point having the highest probability of existence in each group, and the minimum norm solution is obtained again.
The moving distance of each lattice point at this time is a very small distance, and is determined by, for example, applying a predetermined coefficient to the distance between each lattice point. By repeating the above operation, the solution is made to converge on some local regions.

The above is the principle of the grid point moving minimum norm method. FIGS. 9 and 10 schematically show how grid points are grouped. In FIG. 9, N is a lattice point group initially set, and N ₁ and N ₂ are new lattice point groups obtained by moving the lattice point group N in groups. Further, N ₃ , N ₄ , N ₅ , N ₆ , N
Reference numeral ₇ denotes a lattice point group obtained based on further dividing the lattice point groups N ₁ and N ₂ into groups.

[0026]

Although the grid point moving minimum norm method is excellent in that a solution closer to a true solution can be obtained, the following new problem has been clarified. .

That is, in the grid point moving minimum norm method,
The parameters α, β, and γ need to be appropriately set empirically. These parameters are the location of the current source,
Since it depends on the size, direction, etc., it is not easy for an inexperienced operator to appropriately set the values of the respective parameters using the isomagnetic field map as described above. Incorrect settings can result in results that deviate from the true solution.

The present invention has been made in view of such circumstances, and among the above-mentioned parameters α, β, γ, particularly, a parameter (movement parameter) α related to movement on each grid point. It is an object of the present invention to provide a method for estimating a biological activity current source that can easily and accurately estimate a current source without setting.

[0029]

Means for Solving the Problems In order to solve the above problems, the present inventors have made intensive studies and focused on the norm of the solution obtained by the minimum norm method as a criterion for determining the movement parameter α. In the lattice point moving minimum norm method, the minimum norm solution is obtained by moving the lattice point in each iteration, so the norm of the solution changes every time. Therefore, various values were set for the movement parameter α, estimation was performed by simulation, and changes in the norm of the solution were examined. As the sensor, a vector measurement pickup coil schematically shown in FIG. 2 was used. 13 × 3 sensor channels
In the 39 channels, the coils were arranged on a spherical surface having a coil pitch of 37.5 mm and a radius of 117 mm such that the axis of each coil was directed toward the center of the sphere. However, since the tangential component of the magnetic field is affected by the volume current, the magnetic field generated by the current source is calculated using a spherical model and the effect of the volume current is considered (J.Sarvas,
Phys. Med. Biol., Vol.32, pp11-22, 1987).

Assuming that the head is a sphere having a radius of 80 mm, the magnetic sensors are arranged symmetrically about the Z axis (FIG. 11).
(A)). The distance from the surface of the sphere to the sensor is 37m
m. The current source has two current dipoles of position [mm] moment [nAm] (20, 0, 50) (0, 10, 0) (-20, 0, 5) (0, 10, 0). Was.

FIG. 4 shows changes in the norm when the movement parameter α is 0.3 and 0.5. The horizontal axis in FIG. 4 is the number of repetitions, and the vertical axis is the norm of the solution. FIG. 5 shows the estimation results. When α is 0.3, the norm of the solution diverges as shown in FIG. 4A, resulting in a different estimation result as shown in FIG. 5A. On the other hand, when α is 0.5, FIG.
As shown in FIG. 5B, the norm of the solution converges, and a correct estimation result is obtained as shown in FIG.

From this, it can be seen that in each iteration, a grid point movement parameter α that minimizes the norm of the solution is obtained, and the movement of each grid point may be performed using this parameter.

The present invention based on the above findings has the following configuration. That is, the present invention estimates a physical quantity including any of the position, size, and direction of the biological activity current source by detecting a magnetic field generated by the biological activity current with a plurality of magnetic sensors.
A large number of grid points are set in the subject, and the relational expression between the unknown current source on each grid point and the magnetic field data measured by each of the magnetic sensors is a vector having the current source at each grid point as an element. In the biological activity current source estimation method using a method (minimum norm method) of solving by adding a condition that minimizes the norm of
(A) The relationship between an unknown current source on each grid point and magnetic field data obtained by detecting a magnetic field generated by a biological activity current with a plurality of magnetic sensors by setting a large number of grid points in the subject. A first step of solving a formula by the minimum norm method to obtain a current source on each grid point; and (b) calculating each of the grid points based on the physical quantity of the current source at each grid point obtained in the first step. A second step of determining the accuracy of the presence of the current source; and (c) determining the shape of a function (group function) for dividing the grid point group into a plurality of groups based on the accuracy determined in the second step. A third step of optimizing parameters to be determined (moving parameters) according to a condition that a norm of a solution obtained by a minimum norm method is minimized; and (d) the second step.
Using a group function determined by the accuracy obtained in the process and the movement parameter optimized in the third process,
A fourth step of dividing the grid point group into a plurality of groups, and (e) other grid points near the grid point where the current source having the largest size exists in each of the divided grid point groups. A fifth process of rearranging the grid point group by moving
(F) a sixth step of obtaining a current source on each of the rearranged grid points by using the minimum norm method.

[0034]

The operation of the present invention is as follows. The current source estimated in the first step is not a true current source, but a current source close thereto. Therefore, if another grid point group is moved to the vicinity of the grid point where the estimated current source is located, the appropriate grid point arrangement is reconstructed and the minimum norm solution is obtained, without increasing the number of grid points. In addition, the estimation accuracy of the current source can be improved. However, when there are a plurality of true current sources, it is important to determine which grid point is closer to each grid point. Therefore, in the second process, the accuracy of the existence of the current source at each grid point is determined. In the next third step, a movement parameter for determining the shape of a group function for dividing and moving a grid point is defined as:
Optimization is performed according to the condition that the norm of the minimum norm solution is minimized. In a fourth step, the group of grid points is divided into a plurality of groups using a group function determined by the optimized movement parameters. Then, in the fifth step, for each divided grid point group, another grid point is moved to the vicinity of the grid point where the current source having the maximum size exists. In the sixth step, a current source on each of the rearranged grid points is obtained by the minimum norm method. By repeatedly performing the above processing, even if there are a plurality of true current sources, each current source is easily and accurately estimated.

[0035]

An embodiment of the present invention will be described below with reference to the flowchart of FIG.

First, as described with reference to FIG.
, The orthogonal three-directional components of the minute magnetic field from the subject M are simultaneously measured (vector measurement) by the multi-channel SQID sensor 1 (step N1).
However, the magnetic sensors (pickup coils) S _{1 to} S _{m of the} multi-channel SQID sensor 1 used here.
Is composed of three pickup coils each having detection sensitivity in three orthogonal directions. As this type of pickup coil, for example, there is a triaxial gradiometer. In the gradiometer, the pickup coil is divided into two oppositely wound coils, thereby canceling a uniform magnetic field and detecting only a magnetic field having a gradient.

FIG. 2 schematically shows a configuration of a triaxial gradiometer. The pickup coils L _X , L _Y , L _Z are:
The magnetic field components in the X, Y, and Z directions are detected. The configuration of the triaxial gradiometer is not particularly limited. For example, a triaxial gradiometer having a triaxial coil mounted so as to be orthogonal to six surfaces of a cubic core material, and Japanese Unexamined Patent Publication No. Hei.
No. 1, the low temperature resistant flexible material is wound in a cylindrical shape, and the flexible material is interconnected on the surface of the flexible material at a predetermined width in parallel with each other and in reverse winding. The superconducting thin-film coil pairs are set at different angles from each other.
A triaxial gradiometer or the like formed as a pair is used.

Next, similarly to the conventional example shown in FIG. 7, a three-dimensional lattice point group N is set evenly in, for example, the brain, which is a diagnostic symmetric region (step N2).

Then, using the minimum norm method described above,
A current source (minimum norm solution) at each grid point is obtained (step N3).

Next, the accuracy of the existence of the current source on each grid point is determined by the following equation (9) (step N4). Q (Vr _j ) = | VP _j | ……… (9)

In the present embodiment, in step N1, three orthogonal components of the magnetic field are simultaneously vector-measured, and in step N4, the accuracy of the presence of the current source on each grid point is calculated by equation (7).
The reason obtained by the equation (9) that does not include the second term is as follows.

Generally, as shown in FIG. 11 (a), each magnetic sensor for detecting a magnetic field generated by a current source in a living body has a structure in which the coil axis of each magnetic sensor S is a sphere of the subject M. They are arranged radially.
As shown in FIG. 11B, a pair of coils L ₁ and L ₂ are provided in each magnetic sensor S in the radial direction of the sphere (FIG. 11B).
In the Z direction), the detected magnetic field data is only the Z direction component. In other words, 3
Since a magnetic field having components in the directions X, Y, and Z is detected only in the Z direction, the detected magnetic field data has low mutual independence and low spatial resolution. result,
It is considered that the second term in the equation (7), that is, the density of the current source near the lattice point has a large effect on the existence probability Q of the solution on each lattice point.

On the other hand, when the magnetic field generated by the current source of the living body is vector-measured, three orthogonal magnetic fields from the subject are obtained.
Since the directional components X, Y, and Z are detected, the independence between the measured magnetic field data is increased, and as a result, the spatial resolution is improved. Therefore, even without considering the second term of equation (7),
It is considered that the existence probability Q of the solution on each grid point can be accurately obtained by only the first term.

As in the present embodiment, the magnetic field from the subject is vector-measured and the accuracy of the presence of the current source on each grid point is determined by equation (9), which is included in the second term of the equation. There is no need to set empirical parameters β and γ, and the current source can be easily and accurately estimated. Under the same conditions as the simulation described in FIG. 5, the current source was estimated by simulation by applying Equation (9) to each magnetic field data of vector measurement and radial measurement.

FIG. 3 shows the results. FIG. 9A shows the estimation result of the current source based on the magnetic field data measured by the vector.
FIG. 7B shows an estimation result based on the measurement in the radial direction. In the figure, the circles indicate the set positions of the current sources, and the arrows indicate the estimated current sources. As is apparent from FIG. 3, according to the embodiment method of estimating the current source by applying the above-described equation (9) to the vector-measured magnetic field data, the current source is correctly estimated and measured in the radial direction. Equation for the magnetic field data
When (9) is applied, the current source is not correctly estimated.

However, the present invention is not necessarily limited to such a method, and a single-axis type gradiometer is used in step N1 as in the previously proposed minimum grid point moving norm method. Only the radial component of the magnetic field may be measured, and in step N4, the accuracy of the existence of the current source on each grid point may be obtained using equation (7). According to this method, empirical setting of the parameters β and γ is inevitable, but by using the automatic adjustment of the movement parameter α described later, the estimation of the current source can be made easier and the accuracy can be improved accordingly. Can be.

When the accuracy of the existence of the current source on each grid point is determined in step N4, the process proceeds to step N5, where the movement parameter α in the group function defined by equation (8) is used to divide and move the group of grid points. To optimize. Here, the following equation (1
The optimization of the movement parameter α is performed using the evaluation function f represented by (0).

[0048]

(Equation 5)

In equation (10), VP _j (α) is a solution obtained by moving the lattice points using the movement parameter α and obtained by the minimum norm method. Therefore, the evaluation function f is the norm of the solution. is there. N represents the number of grid points.

In steps N5 to N7, several parameters α ₁ , α ₂ , α ₃ ,.
By temporarily moving the grid points using these parameters, the respective minimum norm solutions are obtained. Then, these minimum norm solutions are given to equation (10), and the evaluation functions f (α ₁ ), f
The values of (α ₂ ), f (α ₃ ),... Are determined, and the parameter having the minimum value is adopted as the movement parameter α.
Accordingly, the minimum norm solution obtained by moving the lattice point based on the adopted movement parameter α is adopted, and the minimum norm solution based on other parameters is discarded.

In step N8, the norm of the solution (value f _L-1 of the evaluation function) obtained by the previous processing of steps N4 to N7 and the norm of the current solution (value f _{L of the} evaluation function) Is obtained. If the amount of change Δf of the norm of this solution is equal to or less than a predetermined value,
If not, the process returns to step N4 to repeatedly execute the processes of steps N4 to N7 until the change amount Δf of the norm of the solution becomes equal to or less than a predetermined value.

According to the method of this embodiment, not only the empirical parameters α, β, and γ are not required, but also the norm of the minimum norm solution is used as a condition for stopping the repetitive processing (step N8). The current source estimation process can be stopped at an appropriate number of repetitions.

<Simulation> A simulation was performed to confirm the effectiveness of the method of automatically adjusting the movement parameter α. As the pickup coil, the vector measurement coil shown in FIG. 2 was set at 19 points on the spherical surface. Therefore, the number of channels is set to 57 channels of 19 × 3.
This was set on a spherical surface with a coil pitch of 25 mm and a radius of 131 mm. The head was a sphere with a radius of 80 mm, and the calculation of the magnetic field was performed in consideration of the effect of the volume current. Sensor is Z
They were arranged symmetrically about the axis, and the distance from the sphere on the head to the sensor was 36 mm. The current source set multiple current dipoles.

Assuming the cerebral cortex in the head sphere, three current dipoles were set thereon. The position and moment of each current dipole are as follows. Position [mm] Moment [nAm] (-27.08, 4.78, 47.63) (-8.53, 1.50, -5.00) (-27.08, -4.78, 47.63) (-8.53, -1.50, -5.00) (4.91, -56.08, 32.50) (0.44, -4.98, -8.66)

FIG. 6 shows the result of current source estimation by the grid point moving minimum norm method while automatically adjusting the moving parameter α.
Shown in As shown in the figure, the current source is estimated near the correct position as set. When current source estimation is performed by this method, a solution can be obtained as a current dipole distribution even if the true current source is a single current dipole, but when the estimated current dipole moment is integrated, The moment value of the set current dipole was almost the same, confirming the effectiveness of this method.

[0056]

As is apparent from the above description, according to the present invention, in the grid point moving minimum norm method, the grid point group is divided into a plurality of grid points on the basis of the probability that a current source exists on each grid point. Since the movement parameters that determine the shape of the group function for dividing into groups are automatically adjusted according to the condition that the norm of the minimum norm solution is minimized, the movement parameters are set empirically as before. It is not necessary, and the current source can be easily and accurately estimated.

[Brief description of the drawings]

FIG. 1 is a flowchart of an embodiment.

FIG. 2 is a schematic diagram of a vector measurement pickup coil.

FIG. 3 is a diagram showing simulation results of vector measurement and radial measurement.

FIG. 4 is a diagram illustrating a change in a norm of a solution according to a value of a movement parameter.

FIG. 5 is a diagram illustrating a result of a simulation according to a value of a movement parameter.

FIG. 6 is a diagram illustrating a result of a simulation of a process according to the embodiment;

FIG. 7 is an explanatory diagram of a conventional minimum norm method.

FIG. 8 is an explanatory diagram of grouping of grid points in the grid point moving minimum norm method.

FIG. 9 is a diagram showing a state where a grid point is divided and moved.

FIG. 10 is a diagram showing a state where a grid point is further divided and moved.

FIG. 11 is a diagram showing an arrangement of a general magnetic sensor and its configuration.

[Explanation of symbols]

1 ... Multichannel SQID sensor S ₁ to S _m ... magnetic sensor N ... first grid point group N ₁ to N ₇ ... splitting grating point group _{L X, L Y, L Z} ... vector measurement pickup coil

Claims (1)

(57) [Claims] 1. A method for estimating a physical quantity including any one of a position, a size, and a direction of the biological activity current source by detecting a magnetic field generated by the biological activity current with a plurality of magnetic sensors. A large number of grid points are set, and the relational expression between the unknown current source on each grid point and the magnetic field data measured by each of the magnetic sensors is minimized with the norm of the vector having the current source at each grid point as an element. In the biological activity current source estimating method using a method (minimum norm method) for obtaining a physical quantity of each current source by adding a condition of (a), a number of grid points are set in the subject. An unknown current source on each grid point and a relational expression between magnetic field data obtained by detecting a magnetic field generated by a biological activity current with a plurality of magnetic sensors by solving the minimum norm method on each grid point. Current source (B) a second step of determining the accuracy of the current source at each grid point based on the physical quantity of the current source at each grid point obtained in the first step; and (c) Based on the accuracy obtained in the second step, a parameter (moving parameter) for determining a shape of a function (group function) for dividing the grid point group into a plurality of groups was obtained by a minimum norm method. A third step of optimizing according to the condition of minimizing the norm of the solution, and (d) using a group function determined by the accuracy obtained in the second step and the movement parameter optimized in the third step A fourth step of dividing the grid point group into a plurality of groups, and (e) another grid near the grid point where the current source having the largest size exists in each of the divided grid point groups. Move the point to the case A fifth step of repositioning the point cloud, a sixth process of obtaining by using the minimum norm method a current source on the grid points which are the rearranged (f),
A living activity current source estimating method, comprising: