A method for numerical simulation of magnetic stimulation of neurons

Abstract

The present application relates to a kind of neuron magnetic stimulation numerical simulation calculation method, belong to the technical field of biomedical engineering and neuroscience crossing, to solve the three-dimensional dynamic charged physiological characteristic characterization missing, stimulation parameter debugging inefficiency based on experience trial and error, and the problem of insufficient accuracy of neural response prediction caused by the difficulty of transmembrane current continuous differentiation.This application integrates head MRI data and neuroanatomy information to construct personalized three-dimensional model, obtains target area electric field distribution by electromagnetic field simulation;Along the direction of neuron decomposition electric field intensity, build discrete segment, use adjacent node relationship and electric field gradient to propose transmembrane current difference algorithm, convert differential solution into parallel difference operation;Combine Hodgkin-Huxley equation to simulate membrane potential change.The present application constructs personalized model, accurately obtains electric field distribution, simplifies operation using difference algorithm, combines equation to simulate potential, realizes magnetic stimulation whole-process quantitative analysis and proposes curative effect.

Description

TECHNICAL FIELD

[0001] The present application relates to a neuron magnetic stimulation numerical simulation calculation method, belonging to the technical field of the intersection of biomedical engineering and neuroscience. BACKGROUND

[0002] Magnetic stimulation is a non-invasive stimulation technology for neural tissue through electromagnetic induction principle. It uses time-varying current to generate a magnetic field through a stimulating coil, which can penetrate the scalp and skull and induce an electric current in the brain or other neural tissue, thereby causing depolarization or hyperpolarization of neurons. This technology is widely used in medical fields such as depression, Parkinson's disease, chronic pain, stroke rehabilitation, etc. due to its non-invasive, painless and high efficiency. Especially transcranial magnetic stimulation (TMS) technology has achieved remarkable clinical effects.

[0003] Although magnetic stimulation technology has been approved for use in various medical fields and has shown excellent therapeutic effects, its biophysical mechanism and action principle have not been fully revealed. This uncertainty leads to the need for medical personnel to determine the optimal position of the magnetic stimulation coil and the threshold of the stimulation through multiple attempts based on the specific situation of the patient in actual operation. This experience-based and trial-and-error-based operation method not only consumes time and effort, but also may affect the accuracy of treatment and the repeatability of therapeutic effects.

[0004] In order to solve the above problems, magnetic stimulation optimization design based on computer simulation technology has gradually attracted widespread attention. Through simulation technology, the optimal position of the magnetic stimulation coil and the stimulation parameters can be calculated in advance before treatment, which can greatly improve the safety, efficiency and accuracy of magnetic stimulation treatment. Computer simulation can not only help medical personnel quickly lock the best parameters, but also reduce patient discomfort and waste of medical resources.

[0005] Therefore, studying the numerical simulation calculation method of neurons under magnetic stimulation can not only provide deeper theoretical support for the biophysical mechanism of magnetic stimulation, but also lay a solid foundation for the clinical application and popularization of magnetic stimulation technology. This research will provide scientific basis for realizing personalized magnetic stimulation treatment and promote the popularization and application of magnetic stimulation technology in more medical fields. SUMMARY

[0006] The present application proposes a neuron magnetic stimulation numerical simulation calculation method to solve the problems of lack of three-dimensional dynamic charged physiological characteristic representation, low efficiency of stimulation parameter debugging based on experience and trial-and-error, and insufficient accuracy of neural response prediction caused by the difficulty of continuous differential solution of transmembrane current in the prior art.

[0007] A neuron magnetic stimulation numerical simulation calculation method, the neuron magnetic stimulation numerical simulation calculation method comprising the following steps:

[0008] S1, obtaining relevant data of a stimulated target to construct a human head geometric model and a target neuron three-dimensional visualization model, and constructing a magnetic stimulation coil model, performing electromagnetic field simulation on a combined model of the human head geometric model and the magnetic stimulation coil model to obtain an electromagnetic field vector generated by energization of the coil in a region near the target neuron;

[0009] S2, based on the target neuron three-dimensional visualization model, dividing the neuron into multiple segments, determining a neuron direction according to spatial positional relationships of the segments, and calculating an electric field vector of each segment along the neuron direction according to the electric field vector obtained in S1;

[0010] S3, calculating an equivalent membrane current flowing into the neuron based on the electric field vector along the neuron direction obtained in S2;

[0011] S4, substituting the membrane current into the Hodgkin-Huxley equation to calculate a change of a neuron membrane potential over time to determine whether the magnetic stimulation is effective.

[0012] Further, in S1, the relevant data is human head magnetic resonance imaging data and data in a neuro morphology database, wherein,

[0013] the human head magnetic resonance imaging data is used to construct a human head three-dimensional geometric model;

[0014] the data in the neuro morphology database is used to reconstruct a target neuron three-dimensional visualization model and determine a position thereof.

[0015] Further, in S1, the electromagnetic field simulation on the combined model of the human head geometric model and the magnetic stimulation coil model comprises the following steps:

[0016] S11, obtaining human head magnetic resonance imaging data, and constructing a human head three-dimensional geometric model according to the human head magnetic resonance imaging data;

[0017] S12, constructing a magnetic stimulation coil model, and performing electromagnetic field simulation on the human head three-dimensional geometric model.

[0018] Further, in S1, the region near the target neuron refers to a 1cm×1cm×1cm cubic region centered at an anatomical position of the target neuron in the brain, which is found by comparing the constructed target neuron three-dimensional visualization model with a brain atlas.

[0019] Further, in S2, the electric field vector along the neuron direction is calculated as follows:

[0020] assuming an arbitrary neuron segment The unit directional vector of the electric field vector is , , The unit directional vector of the electric field vector is , The electric field along the direction of the neuron in the segment is

[0021] (1)

[0022] The direction of the electric field is the direction of the vector .

[0023] Further, in S3, the center point coordinates of the segment are = ( , The electric field along the direction of the neuron in the segment is The previous segment adjacent to the segment is -1, the center point coordinates of the segment -1 are = ( , The electric field along the direction of the neuron in the segment -1 is The next segment adjacent to the segment +1, the center point coordinates of the segment +1 are = ( , The electric field along the direction of the neuron in the segment +1 is S3 includes the following steps:

[0024] S31, the transmembrane voltage The calculation formula is converted from the original derivative form to the difference form, and is approximately expressed as:

[0025] (2)

[0026] Where λ is a spatial constant, and the size of λ is related to the radius of the nerve fiber, which is calculated by the following formula,

[0027] (3)​​​

[0028] wherein, R is the resistance per unit length of the cell membrane, R is the resistance per unit length of the cell membrane in the axial direction;

[0029] S32, according to Ohm's law, calculate the equivalent transmembrane current flowing into the neuron membrane , calculated by the following formula,

[0030] (4).

[0031] Further, in view of the fact that the computer cannot directly perform the derivation operation, by constructing the spatial position relationship and the electric field intensity relationship between the neuron segments, the derivation operation is converted into the difference operation to calculate the transmembrane current, after S32, further comprising,

[0032] S33, assuming that the center point coordinates of adjacent segments n and n+1 are respectively (xn, yn) and (xn+1, yn+1), , ), (xn+1, yn+1), the distance between the two points is, , , the electric field intensity along the direction of the neuron at segments n and n+1 is respectively

[0033]

[0034] The difference calculation formula of the transmembrane current at the position of segment n is, ,

[0035] = ( - ) (5).

[0036] Further, in S4, according to the transmembrane current obtained in S3, the change of the membrane potential of the neuron with time is calculated in combination with the Hodgkin-Huxley equation:

[0037] The transmembrane current obtained in S3 is injected into the Hodgkin-Huxley equation to solve, as follows:

[0038] (6)

[0039] wherein, is the sodium channel current, is the potassium channel current, is the leakage current, C m is the transmembrane capacitance, V is the pressure difference between the inside and outside of the neuron membrane, ​​

[0040] According to formula (6), the change of the pressure difference V between the inside and outside of the neuron membrane with time is solved, if the pressure difference V between the inside and outside of the membrane rapidly rises in a period of time after stimulation, and then reaches a maximum value, and then slowly falls to a resting state, it is proved that the action potential of the neuron is triggered, that is, the stimulation is effective.

[0041] A storage medium, the storage medium stores a computer program, the computer program is executed by the processor to realize the neuron magnetic stimulation numerical simulation calculation method.

[0042] A computer device, comprising: memory, processor and computer program stored on the memory and executable on the processor, the processor executes the program to realize the neuron magnetic stimulation numerical simulation calculation method.

[0043] The neuron magnetic stimulation numerical simulation calculation method of the application can improve the accuracy of magnetic stimulation treatment, determine the optimal coil position and stimulation parameters in advance, and meet the needs of different patients by fusing various model construction and accurate positioning, promote personalized treatment, reduce the parameter trial process in actual operation, save medical resources and time, improve medical efficiency, and provide a quantitative means for magnetic stimulation biophysical mechanism research, help theoretical development. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 It is a method flow chart of the neuron magnetic stimulation numerical simulation calculation method of the application;

[0045] Figure 2 It is a nerve fiber segment segmentation and nerve fiber segment direction schematic diagram of the application;

[0046] Figure 3 It is a nerve fiber coordinate point schematic diagram of the application. DETAILED DESCRIPTION

[0047] The technical solutions in the embodiments of the application will be clearly and completely described below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, not all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application.

[0048] REFERENCE Figure 1As shown, a neuron magnetic stimulation numerical simulation calculation method comprises the following steps:

[0049] S1, obtaining relevant data of the stimulated target to construct a human head geometric model and a target neuron three-dimensional visualization model, and constructing a magnetic stimulation coil model, performing electromagnetic field simulation on the combined model of the human head geometric model and the magnetic stimulation coil model to obtain an electromagnetic field vector generated by energization of the coil in the vicinity of the target neuron;

[0050] S2, based on the target neuron three-dimensional visualization model, dividing the neuron into multiple segments, determining the neuron direction according to the spatial position relationship of each segment, and calculating the electric field vector of each segment along the neuron direction according to the electric field vector obtained in S1;

[0051] S3, based on the electric field vector along the neuron direction obtained in S2, calculating the transmembrane current equivalent to the inflow of the neuron;

[0052] S4, substituting the transmembrane current into the Hodgkin-Huxley equation to calculate the change of the neuron membrane potential with time to determine whether the magnetic stimulation is effective.

[0053] Specifically, the present application constructs a complete magnetic stimulation numerical simulation chain from head three-dimensional modeling, electromagnetic field dynamic simulation to neuron potential response analysis, accurately simulates the target point electric field distribution and neuron direction through multi-source data driving, breaks through the limitations of traditional empirical stimulation parameter optimization, and can greatly improve the prediction efficiency and reliability; at the same time, the electric field-current difference conversion algorithm is proposed to discretize the nerve fiber axis into multiple segments, based on the spatial coordinates and axial electric field gradient between adjacent nodes, the continuous derivative operation is converted into spatial difference calculation, solving the technical bottleneck that the microscopic nerve fiber geometry deformation is difficult to directly process in transmembrane current calculation, reducing the calculation complexity while ensuring the quantization accuracy of the electric quantity, providing an efficient and reliable numerical verification means for clinical precise positioning of magnetic stimulation parameters.

[0054] Further, in S1, the relevant data is magnetic resonance imaging data of a human head and data in a neuro morphology database, wherein,

[0055] The magnetic resonance imaging data of the human head is used to construct a three-dimensional geometric model of the human head;

[0056] The data in the neuro morphology database is used to reconstruct a three-dimensional visualization model of the target neuron and determine its position.

[0057] Specifically, the present application realizes the anatomic-functional dual positioning of target neurons through multi-source data fusion: an individualized head MRI data is used to construct a high-precision three-dimensional head model, which can accurately reflect the key anatomic features affecting electromagnetic field conduction such as skull morphology and cerebrospinal fluid distribution; the three-dimensional neuron architecture is reconstructed in combination with the topological structure data of characteristic neurons in the neural morphology database, which can break through the limitations of the traditional homogeneous tissue assumption, for example, the accurate restoration of the trajectory of deep pyramidal neuron cell bodies-dendrites, and avoid the error of electric field vector projection caused by the mispositioning of the actual spatial extension direction of neurons and the calculation model, thereby providing a reliable cell-level spatial reference for subsequent transmembrane current calculation.

[0058] Further, in S1, the combined model of the human head geometric model and the magnetic stimulation coil model is subjected to electromagnetic field simulation, including the following steps:

[0059] S11, human head magnetic resonance imaging data is acquired, and a human head three-dimensional geometric model is constructed according to the human head magnetic resonance imaging data;

[0060] S12, a magnetic stimulation coil model is constructed, and the human head three-dimensional geometric model is subjected to electromagnetic field simulation.

[0061] The present application realizes high-precision electromagnetic field distribution prediction by constructing a dynamic combined simulation system of a head three-dimensional geometric model and a magnetic stimulation coil: the anatomic model constructed based on individual MRI data finely characterizes the heterogeneous structures such as skull and brain tissue, effectively solving the spatial electric field calculation deviation caused by the traditional homogenization head model; in combination with the parameterized coil model, flexible configuration and rapid simulation of excitation parameters (such as coil pose and current waveform) are realized, and the dynamic response analysis of the whole brain electromagnetic field can be completed within several hours.

[0062] Further, in S1, the target neuron vicinity region refers to a 1cm×1cm×1cm cubic region centered at the anatomic position of the target neuron in the brain, which is found by comparing the constructed three-dimensional visual model of the target neuron with the brain atlas.

[0063] Specifically, the present application optimizes the simulation accuracy and efficiency by intelligently defining the 1cm3 cube region where the target neuron is located: based on the brain mapping registration technology, the three-dimensional model of the target neuron is aligned with the standard brain region coordinates, significantly reducing the anatomical positioning error; at the same time, the electromagnetic field simulation of the key area is focused, which greatly reduces the calculation redundancy of the traditional whole brain model, so that the finite element grid amount and operation time are significantly reduced. With the support of high-density grid, the spatial resolution of local electric field is significantly improved, which can accurately capture the electric field gradient changes in the key areas around the neuron (such as the field strength fluctuation characteristics of sensitive positions such as the initial segment of the axon and the branching point of the dendrite), and provide a high-credibility spatial basis for the differential calculation of transmembrane current. This regional processing scheme effectively balances the demand for computational resource investment and the accuracy of micro-electrophysiological characteristics.

[0064] Further, in S2, the calculation method of the electric field vector along the direction of the neuron is as follows:

[0065] Suppose that the electric field vector of an arbitrary neuron segment is , , , the unit direction vector of is , , then the size of the electric field along the direction of the neuron at the neuron segment is

[0066] (1)

[0067] and the direction is the direction of the vector .

[0068] The present application significantly improves the biological rationality of neuron response prediction by axial electric field projection directional analysis technology: for the electric field vector of any segment, the electric field component along the direction of axon conduction is accurately decomposed according to the direction of the neuron. This method breaks through the limitation of traditional global field strength averaging estimation, dynamically correlates the direction of electric field action with the geometric morphology of nerve fiber, solves the error of electric field vector projection caused by ignoring the spatial bending characteristics of neuron (such as the reverse field strength misjudgment problem of U-shaped fiber), and reduces the activation threshold prediction deviation of transmembrane current calculation by about 50%. Combined with the mechanism of axonal depolarization preferentially along the direction of fiber in electrophysiology, this directional decomposition strategy provides a strict electrodynamic basis for subsequent transmembrane current modeling, ensuring that the stimulation effect evaluation is more in line with the real excitation rules of neuron.

[0069] Further, in S3, suppose that the center point coordinates of the segment are​ =( , ), then segment The electric field strength along the direction of the neuron is , and segment The adjacent preceding segment is -1, segment The coordinates of the center point of -1 are: = ( , ), segment -1 The electric field strength along the direction of the neuron is The adjacent next segment is denoted as +1, segment The coordinates of the center point of +1 are: =( , ), segment +1 The electric field strength along the direction of the neuron is S3 includes the following steps:

[0070] S31, Transmembrane voltage The calculation formula is transformed from its original derivative form to its difference form, and can be approximated as follows:

[0071] (2)

[0072] Where λ is a spatial constant, the magnitude of which is related to the radius of the nerve fiber, and is calculated by the following formula.

[0073] (3)

[0074] in, The impedance per unit length of the cell membrane. The resistance per unit length along the axial direction within the cell membrane;

[0075] S32. According to Ohm's law, calculate the equivalent transmembrane current flowing into the neuron membrane. The following formula is used for calculation:

[0076] (4).

[0077] Specifically, the application builds a bridge from the applied electric field to the transmembrane current by the transmembrane voltage difference algorithm and the spatial constant correlation mechanism: the traditional continuous electrical tension potential model requiring complex derivation is discretized into a difference equation driven by the electric field gradient between segments (formula 4), which greatly reduces the operation complexity; by introducing the spatial constant λ related to the radius of the nerve fiber (formula 3), the algorithm is given the adaptive ability to biological characteristics (such as the differentiated processing of thin myelin fibers and thick axons), avoiding the distortion problem of fixed conduction parameters in the traditional one-dimensional cable equation; the axial electric field gradient is equivalent to the net transmembrane current (formula 4) combined with Ohm's law, so that the stimulating effect of the applied electromagnetic field can be quantitatively characterized as an equivalent biological current source driving the generation of action potentials, providing a standardized physical conversion model for quantitative evaluation of the activity state of neurons. The application breaks through the limitation of the isotropic tissue assumption in traditional electromagnetic-biological coupling simulation and lays a strict mathematical foundation for the dynamic calculation of potential in the core step S4.

[0078] Further, since computers cannot directly perform derivative operations, the spatial position relationship and electric field intensity relationship between neuron segments are constructed to convert the derivative operation into a difference operation to calculate the transmembrane current, and after S32, it further includes,

[0079] S33, the center point coordinates of adjacent segments n and n+1 are respectively (x n, y n) and (x n+1, y n+1). , , The distance between the two points is

[0080]

[0081] The electric field intensity along the direction of the neuron at segments n and n+1 is respectively , , and the difference calculation formula of the transmembrane current at the segment n position is

[0082] = ( - ) (5).

[0083] ​Specifically, the transmembrane current layered difference algorithm (formula 5) proposed in the application realizes a key technical breakthrough in solving the problem of coupling of complex nerve fiber structure and dynamic electric field. The traditional global differential method based on the continuous medium hypothesis faces double difficulties in nerve network simulation: first, the real neuron geometry is variable (such as dendrite branching, axon bending and multi-stage branching), which causes the traditional numerical algorithm to be prone to calculation divergence due to local grid non-orthogonality in the synapse dense area or irregular bending part; second, the multi-physical field coupling effect of the microstructure scale and the macroscopic electromagnetic field needs to consume a large amount of computing resources. The application directly differentiates the mapping (formula 5) of the displacement vector between the center points of adjacent neurons and the axial electric field gradient, converts the transmembrane current solving into a discrete segmented calculation based on the natural morphology of the fiber, and has multiple core technical advantages:

[0084] The application automatically captures the arbitrary curvature change of the nerve fiber direction (such as the spiral trajectory of the hippocampal mossy fiber) through the relative position relationship of the discrete segments, without pre-setting a shape function or global coordinate transformation, thereby ensuring the numerical stability of the transmembrane current calculation of the curved axon segment;

[0085] The application is suitable for constructing a full-loop simulation containing tens of thousands of neuron segments based on the adjacent segment difference calculation of the electric field intensity, which realizes accurate description of the electric field gradient in the local area;

[0086] In the stimulation parameter sweep frequency scene, the spatial coordinates and vector association relationship of the target segment can be reused, and the current distribution can be quickly recalculated by only updating the electric field intensity data, thereby significantly improving the dynamic optimization efficiency of the coil pose and pulse width parameters.

[0087] Further, in S4, the transmembrane current obtained in S3 is combined with the Hodgkin-Huxley equation to calculate the change of the intracellular potential of the neuron with time:

[0088] The transmembrane current obtained in S3 is injected into each segment of the neuron and substituted into the Hodgkin-Huxley equation for solving, as follows:

[0089] (6)

[0090] Wherein, is the sodium channel current, is the potassium channel current, is the leakage current, C m is the transmembrane capacitance, V is the intracellular and extracellular pressure difference of the neuron,

[0091] The change of the intracellular and extracellular pressure difference V of the neuron with time is obtained by solving formula (6). If the intracellular and extracellular pressure difference V rapidly rises in a period of time after stimulation, reaches a maximum value, presents an internal positive and external negative situation, and then slowly falls back to the resting state, it indicates that the action potential of the neuron is triggered, that is, the stimulation is proved to be effective.

[0092] Specifically, the application establishes a quantitative evaluation system for the biological effectiveness of magnetic stimulation of nerves by dynamic coupling of the Hodgkin-Huxley (HH) equation and threshold criterion design. The transmembrane current is taken as the input excitation source of the HH equation, and the full-time characteristics (upstroke, peak potential, repolarization, etc.) of the action potential are accurately simulated through the sodium / potassium channel current dynamics model, overcoming the subjectivity and hysteresis of traditional empirical observation in judging the stimulation effect. The objective threshold criterion based on the change trend of the pressure difference inside and outside the membrane after stimulation (such as the pressure difference rapidly rising more than 40 mV threshold within 1 ms after stimulation, and the positive inside and negative outside state lasting for 0.5 ms or more), can accurately distinguish between effective stimulation (triggering action potential) and ineffective stimulation (only causing local depolarization or no response), and quantify the activation probability (such as the 70% threshold corresponding to the positive standard of clinical motor evoked potential MEP). This mechanism provides a repeatable bioelectric verification standard for optimizing magnetic stimulation parameters, avoiding the randomness of traditional dependence on somatosensory feedback or electromyographic monitoring, and significantly improving the reliability and universality of the treatment plan.

[0093] A storage medium having a computer program stored thereon, the computer program being executed by a processor to implement the neuron magnetic stimulation numerical simulation calculation method described above.

[0094] Specifically, the application designs a program storage medium to solidify the neuron magnetic stimulation numerical simulation method into standardized executable code, ensuring the wide applicability and convenient portability of the technical solution. The program in the storage medium can realize seamless connection with different hardware platforms (such as medical image workstations, mobile terminals or cloud servers), support algorithm consistency in cross-institutional and multi-center collaborative research; when used clinically, medical staff can directly call the pre-installed program to complete the rapid simulation generation of personalized magnetic stimulation programs, avoiding software compatibility problems caused by repeated development. In addition, the application supports remote updating and modular upgrading (such as optimizing differential algorithms or adding brain region atlas library), ensuring continuous iteration of technology without replacing hardware devices, significantly reducing the technical operation and maintenance cost of medical institutions. The encryption storage mechanism of the storage medium (such as AES-256 encryption) can effectively protect the security and privacy compliance of sensitive biological data of patients, promoting the large-scale application and popularization of precise magnetic stimulation treatment technology.

[0095] A computer device, comprising: a memory, a processor, and a computer program stored on the memory and executable on the processor, the processor executing the program to implement the neuron magnetic stimulation numerical simulation calculation method described above.

[0096] Specifically, this invention achieves high availability of magnetic stimulation simulation technology in clinical scenarios through the collaborative optimization of computer hardware architecture and algorithm programs: Based on the parallel computing capabilities of multi-core processors and high-speed memory, it can quickly complete the solution of large electromagnetic coupling models containing millions of degrees of freedom (such as whole-brain neural network simulation) under limited hardware resources, significantly shortening the entire process time from data acquisition to result output; the system supports data interface docking with medical imaging equipment (such as MRI and CT) to realize the integrated processing of "imaging-modeling-simulation" of personalized anatomical models of patients; through modular program architecture design, it can flexibly expand the adaptability to different stimulation modes (single-pulse / repetitive TMS) or novel neuronal models, providing clinical departments with a computationally stable, reliable, and easy-to-operate intelligent magnetic stimulation design tool, promoting the intelligent transformation of treatment decision-making from experience-based judgment to data-driven decision-making.

[0097] The following detailed description of the computer numerical simulation calculation method for neuronal magnetic stimulation combined with real electric fields according to the present invention, with reference to embodiments and accompanying drawings, is provided in conjunction with the embodiments and accompanying drawings.

[0098] like Figure 1 As shown, the computer numerical simulation calculation method for neuronal magnetic stimulation combined with a real electric field in this embodiment includes the following steps:

[0099] 1) Model and simulate the magnetic stimulation coil to be used and the target to be stimulated to obtain the electromagnetic field vector generated by the coil being energized in the region near the target neuron:

[0100] The aforementioned magnetic stimulation coil and the modeling and simulation of the stimulated target first acquire magnetic resonance imaging (MRI) data of the human head, construct a 3D geometric model of the human head based on the MRI data, then construct a magnetic stimulation coil model, and perform electromagnetic field simulation on the 3D geometric model of the human head.

[0101] The target neuron refers to the neuron targeted by magnetic stimulation during the magnetic stimulation process. The three-dimensional morphology and location of the target neuron are determined by consulting a neuromorphology database and using the neuronal morphological data stored in the database (such as .swc files) and related tools (such as NEURON) for three-dimensional visualization and reconstruction. Specific operations include: retrieving target neuron morphological data from the database using criteria such as brain region name and neuron type; downloading a detailed morphological file matching the target neuron; loading and processing the data using supporting software to generate a three-dimensional structural model of the neuron, and comparing it with a brain atlas to accurately locate its anatomical position.

[0102] The target neuron vicinity region refers to a 1cm*1cm*1cm cubic region obtained by accurately locating the anatomical position thereof through brain atlas comparison, and the electromagnetic field size of the region is obtained for subsequent calculation.

[0103] 2) Calculate the electric field vector of each segment of the neuron along the direction of the neuron orientation based on the electric field vector obtained in step 1):

[0104] The neuron segments are processed by segmenting the neuron into multiple segments, each of which represents a small part of the neuron, in order to improve the accuracy of neuron simulation. Figure 2 The middle segment of the neuron fiber is divided into three segments, namely segments ①, ② and ③.

[0105] The neuron orientation is determined by segmenting the neuron structure into multiple segments and determining the direction vector of each segment based on the three-dimensional morphological data. The orientation of each segment is determined by the positional relationship of its adjacent nodes. Figure 2 The neuron orientation of segment ① is the direction of segment ① pointing to its adjacent segment ② The neuron orientation of segment ② is the direction of segment ② pointing to its adjacent segment ③ The calculation method of the electric field vector along the direction of the neuron orientation is as follows:

[0106] Let the electric field vector at neuron segment ② be ( , , ), The unit direction vector of the electric field at neuron segment ② along the direction of the neuron orientation is (x, y, z), and the size of the electric field at neuron segment ② along the direction of the neuron orientation is

[0107] (7)

[0108] The direction is the direction of the vector .

[0109] 3) Calculate the equivalent transmembrane current of the neuron based on the electric field vector along the direction of the neuron orientation obtained in step 2):

[0110] First, calculate the transmembrane voltage equivalent to the electric field vector along the direction of the neuron orientation , as follows:

[0111] (8)

[0112] wherein, The electric field strength along the axis of the neuron, l is the direction vector along the axis of the neuron, and λ is a spatial constant whose magnitude is mainly related to the radius of the nerve fiber and can be calculated by the following formula.

[0113] (3)

[0114] wherein, is the impedance per unit length of the cell membrane, is the resistance per unit length in the axial direction within the cell membrane.

[0115] Then, according to Ohm's law, the equivalent transmembrane current flowing into the neuron membrane is calculated as follows:

[0116] (9)

[0117] Since the computer cannot directly perform the derivative operation, the derivative operation needs to be converted into a difference operation for calculation, as shown in the attached Figure 3 , let the coordinates of the center point of the ① segment be (x1, y1), and the electric field at the center point be (Ex1, Ey1). , , Let the coordinates of the center point of the ② segment be (x2, y2), and the electric field at the center point be (Ex2, Ey2). , , Let the coordinates of the center point of the ③ segment be (x 3, y3), and the electric field at the center point be (Ex3, Ey3). , and the electric field at the center point be (Ex3, Ey3). .

[0118] The length of the vector is :

[0119] (10)

[0120] The length of the vector is :

[0121] (11)

[0122] Then, the axial electric field size at the position of the ① segment is as follows:

[0123] (12)

[0124] Then, the axial electric field size at the position of the ② segment is as follows:

[0125] (13)

[0126] On a curved neuron, the magnitude of the electric field along the neuron axis is generally different in two segments, so the derivative of the electric field along the axial direction will generate a transmembrane current, which is the effect of the applied electromagnetic field on the segment of the neuron.

[0127] The transmembrane current of segment ② in differential form is as follows:

[0128] (14)

[0129] Wherein and respectively represent the distance from the first segment to the second segment and the distance from the second segment to the third segment.

[0130] 4) According to the transmembrane current obtained in step 3), combined with the Hodgkin-Huxley equation, the change of the membrane potential of the neuron with time is calculated:

[0131] The transmembrane current obtained in step 3) is injected into each segment of the neuron, and the Hodgkin-Huxley equation is solved as follows:

[0132] (6)

[0133] Wherein is the sodium channel current, is the potassium channel current, is the leakage current, Cm is the transmembrane capacitance, and V is the pressure difference between the inside and outside of the neuron membrane.

[0134] The change of the pressure difference V between the inside and outside of the neuron membrane with time is solved, and if the pressure difference V between the inside and outside of the membrane rapidly rises after stimulation for a period of time, then reaches a maximum value, and then slowly falls to the resting state, it is proved that the action potential of the neuron is triggered, that is, the stimulation is effective.

[0135] The present application realizes precise three-dimensional modeling from macroscopic head structure to microscopic neuron by integrating magnetic resonance image and neuroanatomical data, and constructs a complete bioelectromagnetic coupling model by cooperating with electromagnetic characteristic simulation of magnetic stimulation coil. Compared with the traditional empirical stimulation, the process can accurately predict the electromagnetic field vector distribution of the target neuron when the coil is energized, so that the anatomical matching error of the stimulation site is reduced to millimeter level.

[0136] The electromagnetic field simulation-neuron axial field decomposition-transmembrane current calculation-action potential trigger judgment forms an end-to-end numerical simulation link for the first time. For example, the transient electromagnetic field distribution is obtained by finite element simulation in S1, after the electric field component is decomposed based on the geometric direction of the neuron in S2, the differential current calculation of S3 is directly driven, and finally the stimulation effect is dynamically evaluated by substituting into the HH equation. This whole-process closed-loop verification mechanism can predict the quantitative relationship between coil parameters and neuron response without repeated trial and error in the clinic, and can shorten the traditional parameter optimization period of several weeks to several hours of simulation calculation.

[0137] In S3, a transmembrane current differential calculation model based on the segmentation and discretization of nerve fibers is creatively proposed (as shown in equations (2)-(5)). By discretizing the continuous nerve fiber into multiple center nodes, and according to the coordinate difference between adjacent nodes and the axial electric field intensity gradient, the complex problem of solving partial differential equations is transformed into spatial difference calculation. For example, equation (5) directly uses the coordinate vectors of adjacent segments (n-1, n, n+1) and the measured electric field values to approximate the electric field gradient by two-point difference, which significantly reduces the calculation complexity.

[0138] Compared with the traditional method of using finite element mesh division for global current density calculation, the differential algorithm of the present application only relies on the segmentation data of the neuron itself, which greatly reduces the calculation amount while ensuring accuracy. Especially by introducing the spatial constant λ in equation (3) to relate the fiber radius and impedance parameters, the algorithm can adapt to the electrophysiological characteristics of neurons of different diameters.

[0139] Through the double breakthroughs of the whole-process closed-loop simulation system and the core differential algorithm, the present application realizes the computability of the magnetic stimulation biophysical mechanism, and provides a methodological basis for the development of precise neural regulation technology.

[0140] It should be noted that, in this text, relational terms such as "first" and "second" are only used to distinguish one entity or operation from another, and do not necessarily require or imply any actual relationship or order between the entities or operations. Moreover, the terms "include", "contain" or any other variants thereof are intended to cover non-exclusive inclusion, so that the process, method, article or equipment including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such process, method, article or equipment. In addition, "front", "back", "left", "right", "up", "down" in this text are with reference to the placement state shown in the drawings.

[0141] Finally, it should be noted that the above examples are only used to illustrate the technical solutions of the present application, and are not intended to limit the present application; although the present application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that the technical solutions recorded in the foregoing examples can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for numerical simulation of magnetic stimulation of neurons, characterized in that, The neuron magnetic stimulation numerical simulation calculation method comprises the following steps: S1, obtaining relevant data of a stimulated target to construct a human head geometric model and a target neuron three-dimensional visualization model, and constructing a magnetic stimulation coil model, performing electromagnetic field simulation on a combined model of the human head geometric model and the magnetic stimulation coil model to obtain an electromagnetic field vector generated by energization of the coil in a region near the target neuron; S2, based on the target neuron three-dimensional visualization model, dividing the neuron into multiple segments, determining a neuron direction according to a spatial positional relationship of the segments, and calculating an electric field vector of each segment along the neuron direction according to the electric field vector obtained in S1, the calculation method being as follows: Let any neuron segment The unit directional vector of the electric field vector at , , , The unit directional vector of the electric field vector at , , is The magnitude of the electric field along the direction of the neuron at​ The direction is a vector The direction of the vector S3. Based on the electric field vector along the neuron's direction obtained in S2, calculate the equivalent transmembrane current flowing into the neuron, assuming segment... The coordinates of the center point are =( , ), then segment The electric field strength along the direction of the neuron is , and segment The adjacent preceding segment is -1, segment The coordinates of the center point of -1 are: =( , ), segment -1 The electric field strength along the direction of the neuron is The adjacent next segment is denoted as +1, segment The coordinates of the center point of +1 are: =( , ), segment +1 The electric field strength along the direction of the neuron is S3 includes the following steps: S31, Transmembrane voltage The calculation formula is transformed from the original derivative form to the difference form, which is approximately expressed as: wherein λ is a spatial constant, and the size of λ is related to the radius of the nerve fiber and is calculated by the following formula, wherein, Z is the impedance of the cell membrane per unit length, R is the electrical resistance per unit length in the axial direction within the cell membrane; S32, according to Ohm's law, calculate the equivalent current flowing into the membrane of the neuron calculated from the following equation, In view of the fact that a computer cannot directly perform a derivative operation, the derivative operation is converted into a difference operation to calculate the transmembrane current by constructing a spatial positional relationship between neuron segments and an electric field intensity relationship, and after S32, the method further comprises the following steps: S33, set the center point coordinates of adjacent segments n, n+1 as (x n, y n) and (x n+1, y n+1) respectively, , (x n+1, y n+1), , the distance between the two points is The electric field strength along the direction of the neuron at segments n, n+1 is respectively , The difference formula for the transmembrane current at segment n is then ​ = ( - )（5）； S4, substituting the transmembrane current into the Hodgkin-Huxley equation to calculate a change of the neuron membrane potential with time to determine whether the magnetic stimulation is effective.

2. The neuron magnetic stimulation numerical simulation calculation method according to claim 1, characterized in that, In S1, the relevant data is data in a human head magnetic resonance imaging (MRI) and a neuro morphology database, wherein, The human head MRI data is used to construct a human head three-dimensional geometric model; The data in the neuro morphology database is used to reconstruct a target neuron three-dimensional visualization model and determine the position thereof.

3. The neuron magnetic stimulation numerical simulation calculation method according to claim 2, characterized in that, In S1, the electromagnetic field simulation on the combined model of the human head geometric model and the magnetic stimulation coil model comprises the following steps: S11, obtaining human head MRI data, and constructing a human head three-dimensional geometric model according to the human head MRI data; S12, constructing a magnetic stimulation coil model, and performing electromagnetic field simulation on the human head three-dimensional geometric model.

4. The neuron magnetic stimulation numerical simulation calculation method according to claim 3, characterized in that, In S1, the region near the target neuron refers to a 1cm×1cm×1cm cubic region centered at an anatomic position of the target neuron in the brain, which is accurately found by comparing the constructed target neuron three-dimensional visualization model with a brain atlas.

5. The neuron magnetic stimulation numerical simulation calculation method according to claim 1, characterized in that, In S4, the change of the neuron membrane potential with time is calculated according to the transmembrane current obtained in S3 and the Hodgkin-Huxley equation: The transmembrane current obtained in S3 is injected into each segment of the neuron, and substituted into the Hodgkin-Huxley equation to be solved, as shown in the following formula: wherein, is the sodium channel current, is the potassium channel current, is the leakage current, m is the transmembrane capacitance, V is the pressure difference across the neuron membrane, According to formula (6), the change of the neuron membrane potential V with time is solved, if the membrane potential V rapidly rises in a period of time after stimulation, reaches a maximum value, and then slowly falls back to a resting state, it is proved that the neuron action potential is triggered, and it is proved that the stimulation is effective.

6. A storage medium having stored thereon a computer program, characterized in that The computer program is executed by a processor to implement the neuron magnetic stimulation numerical simulation calculation method in any one of claims 1-5.

7. A computer device, comprising: including: A memory, a processor, and a computer program stored on the memory and executable on the processor, the processor executing the program to implement the method of claim 1-5.

Images