A method for simulating action potential conduction of neurons based on spectral element method

Abstract

The application discloses a neuron action potential conduction simulation method based on a spectral element method, and belongs to the technical field of bio-electromagnetics. The method comprises the following steps: 1) neuron model input; 2) conduction equation construction: a partial differential equation for action potential propagation in a single nerve fiber is established; 3) conduction equation solving: the specific form of the conduction equation established in the step 2) is changed by using the neuron model input in the step 1), and the resting potential is solved by using the spectral element method; 4) the voltage change of each point of the neuron with time is obtained according to time iteration. The spectral element method can simulate the action potential propagation of various forms of neurons, and only the basic parameters of the neuron model need to be imported and the external input is set, so that the propagation process of the action potential in the nerve fiber can be calculated, that is, the voltage change and distribution in and outside the cell membrane in the action potential propagation process. The method has the advantages of fast calculation speed, accurate result, simple operation and the like.

Description

Technical Field

[0001] The present invention belongs to the field of bioelectromagnetic technology, and in particular relates to a neuron action potential conduction simulation method based on a spectral element method. Background Art

[0002] An action potential in a neuron is a diffusible change in electrical potential from its resting potential that occurs when an excitable cell is stimulated. This change in electrical potential propagates along the cell membrane and is called a nerve impulse. Studying how nerve signals propagate within the human body is crucial for understanding the human nervous system.

[0003] Currently, most studies of neuronal action potentials rely on experimental methods, such as patch clamping, which measures the voltage difference between the inside and outside of the membrane using an external circuit to derive the waveform of the action potential. Although physiological methods like patch clamping offer extremely high temporal resolution and accuracy for studying neural electrical signals, accurately recording action potentials, they are highly disruptive to the sample and suffer from low measurement throughput and spatial resolution, making them difficult to meet the needs of neurobiologists. Alternatively, some studies use biological dissection to establish anatomical models of neurons, construct geometric models, and employ complex algorithms to determine action potential propagation. However, the geometric model construction and solution process require significant resources and take a long time to solve. Furthermore, some experimental methods use multi-channel electrodes to measure the voltage difference between the inside and outside of the membrane. While this can improve spatial resolution, it still suffers from certain disadvantages, as it only records extracellular signals, making it difficult to accurately capture the precise action potential propagation process along dendrites. In theoretical research such as modeling and simulation, for example, in the 2017 publication "Mathematics for Neuroscientists," an analysis of nerve fibers and the integration of the transmembrane Hodgkin-Huxley model were used to generalize the action potential propagation equation for a single nerve fiber. Furthermore, the single nerve fiber model was expanded to develop a matrix representation for multiple branches, which was then solved using a finite difference method. While the equation was cleverly formulated, the finite difference method used to solve it was not very precise, resulting in significant errors between sparse and dense grids in complex neuron models. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for simulating neuronal action potential conduction based on the spectral element method to address the problems in the prior art of low accuracy or slow solution speed when solving complex neurons.

[0005] The present invention comprises the following steps:

[0006] 1) Neuron model input;

[0007] 2) Construction of the conduction equation: Establishing the partial differential equation for the propagation of action potentials in a single nerve fiber;

[0008] 3) Solving the conduction equation: Using the neuron model input in step 1), the specific form of the conduction equation established in step 2) is changed, and the resting potential is solved using the spectral element method;

[0009] 4) Obtain the voltage variation of each neuron point over time based on time iteration;

[0010] In step 1), the neuron model input is to import the neuron morphological parameters, electrical parameters, and radius of each point to be simulated. The neuron morphological parameters include the length of each neuron branch and the connection mode between branches. The electrical parameters include the conductivity of sodium ions, potassium ions, and chloride ions at each location, cytoplasmic coupling resistance, etc.

[0011] In step 2), the specific steps of constructing the conduction equation may be:

[0012] (1) After setting the initial values ​​of the equation and the input waveform of the external stimulus, the stimulus time, position, duration, and stimulus intensity, each branch is constructed as a single nerve fiber model, and the electrical parameters and morphological parameters are introduced. Since the voltages at the connection points of adjacent branches are equal and the partial derivatives of the voltage with respect to time are equal, multiple single nerve fibers can be connected together to form neurons.

[0013] (2) When the neuron is at resting potential, set the initial value of the conductivity at each point of the neuron.

[0014] (3) Enter the equation initialization and construct a unit matrix for each solving unit, including the constant term in the equation, the differential terms of each order of voltage, and the stimulation term on the right end; after obtaining the required unit matrix, assemble each unit matrix together to form a global matrix, substitute the model parameters that change with position and construct the matrix equation to be solved.

[0015] (4) The initial value of the solution of the matrix equation is the voltage at each point when there is no stimulation, that is, the value of the resting potential. Next, the spectral element method matrix equation form of the equation is applied with time central difference to obtain the time iteration equation, that is, the left side of the equation is the related term of the voltage at the next moment, and the right side is the related term of the voltage at the current moment. By performing time iteration, the change of the voltage at each point of the neuron over time can be obtained.

[0016] In step 3), the basic equation for action potential conduction is first constructed, which is a second-order partial differential equation for voltage with respect to time and position. The electrical and geometric parameters in the equation are then set according to the imported external model. The second-order partial differential equation is solved using the spectral element method. Specifically, the resting potential in the absence of stimulation is first solved. Then, using the resting potential as the initial time value and combining it with the time iteration method, the voltage at each point on the neuron is iteratively solved at all times to obtain the change in time.

[0017] The voltage of the action potential is solved using the spectral element method, and the Gauss-Lobatto-Legendre (GLL) point interpolation method inside the unit is applied to make the radius change smoother and closer to the actual situation. In the assignment of the radius of the GLL point inside the unit, after the actual model is sampled, the radius obtained by the actual sampling is added to the two endpoints of the unit. According to the position distance between the GLL point inside the unit and the unit endpoint, the size of the radius at the GLL point is obtained by linear interpolation, so that the radius has a relatively smooth change inside the unit.

[0018] The spectral element method simulates the propagation of action potentials in neurons. By importing the geometric and electrical parameters of the neuron model, it solves the propagation of the action potential generated by the applied stimulus in this neuron. Compared with the finite difference solution method, the solution error of the spectral element method decreases exponentially with the increase of the number of nodes N, while the finite difference method generally decreases linearly. That is, under the same conditions, the spectral element method has higher accuracy. In terms of grid division, the spectral element method can freely change the grid density according to the actual situation, using a denser grid in areas with large radius changes, thereby improving accuracy while ensuring speed. In practical applications, when solving complex real neuron models, it uses fewer resources, has higher accuracy, and is faster to solve. It can also use fewer units and higher orders to obtain more accurate results.

[0019] Compared with the prior art, the outstanding technical effects and advantages of the present invention are:

[0020] The present invention is applied in the field of bioelectromagnetism and has a very significant role and value in understanding the working mode of biological nervous systems and the transmission process of neural signals. The spectral element method can be used to simulate the action potential propagation of neurons of various forms. It only needs to import the basic parameters of the neuron model and set the external input to calculate the propagation process of the action potential in the nerve fiber, that is, the change and distribution of the voltage inside and outside the cell membrane during the action potential propagation. Compared with other solution methods, such as the finite difference solution method, the spectral element method has the advantages of fast calculation speed, accurate results, and simple operation. In addition, when importing a neuron model with a low discrete radius sampling rate, the simulation results are closer to the actual situation. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 This is a diagram illustrating GLL integration point interpolation.

[0022] Figure 2 This is a diagram of neuron morphology.

[0023] Figure 3 It is the waveform diagram of external input stimulation.

[0024] Figure 4 It is the voltage feedback waveform of the neuron soma.

[0025] Figure 5 It is the unit convergence curve of the spectral element to calculate the real neuron voltage and time.

[0026] Figure 6 It is the order convergence curve of the spectral element calculation of the real neuron voltage and time.

[0027] Figure 7 This is the voltage response diagram of each branch of the bifurcated nerve fiber at 7.2ms. DETAILED DESCRIPTION

[0028] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the following embodiments will be further described in conjunction with the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0029] The present invention first imports the neuron model parameters to be solved, such as the length of each branch of the neuron, the connection method between the branches, the radius of each location, the conductivity of sodium ions, potassium ions, chloride ions, cytoplasmic coupling resistance and other parameters. Then set the input waveform of the external stimulation, the stimulation time, position, duration, stimulation intensity, etc., and in this link, you can customize the initial value of each neuron conductivity. Then enter the equation initialization link, and construct a unit matrix for each solving unit, mainly including the constant term in the equation, the differential terms of each order of voltage, and the right-hand end stimulation term. After obtaining the required unit matrix, each unit matrix is ​​assembled together to form a global matrix, and the model parameters that change with position are substituted to construct the matrix equation to be solved. The specific solution process of the equation is as follows.

[0030] The present invention uses the spectral element method to solve the second-order partial differential equation for the conduction of action potentials in neurons. This equation is obtained by superimposing the axial propagation equation and the Hodgkin-Huxley equation. The main form of the conduction equation is as follows:

[0031]

[0032] The voltage V in this equation is not only related to the position, but also changes with time due to the influence of ion concentration.

[0033] In formula (1), V(x,t) is the voltage value of each neuron position at time t, is the time constant, where C m is the capacitance per unit area, g cl is the conductivity per unit area. where R a is the cytoplasmic resistance, and the right end is the input stimulus. Here, the stimulus applied at a certain point changes with time, I stim (t) is the current magnitude that changes with time, and the stimulus term can be replaced by the synaptic stimulus in a real neuron.

[0034] Among them, m t =α m (u)(1-m)-β m (u)m,h t =α h (u)(1-h)-β h (u)h,n t =α n (V)(1-n)-β n (V)n,α m ,β m ,α h ,β h ,α n ,β n is a function related only to voltage, are the maximum conductivity of each ion.

[0035] The spectral element method is used to solve the above equations, using the Gauss-Lobatto-Legendre (GLL) orthogonal polynomial basis function. The N-order GLL basis function in a standard reference unit ξ∈[-1,1] is as follows:

[0036]

[0037] Where j = 0, 1, ... N, L N (ξ) is an Nth-order Lagrange polynomial, which converts the nodes {ξ∈[-1,1] j ,j=0,1…N} as GLL integration points, which are equations (1-ξ 2 )L′ N The N roots with (ξ)=0 are called GLL points because each node is expanded using the GLL basis function.

[0038] Multiplying both ends of equation (1) by the test function w(x) and integrating them yields equation (2):

[0039]

[0040] Arrange and transform equation (2) and include G a Integrate by parts for the term (x). Considering that the end of the neuron branch here is not connected to other parts, the current propagates to this point and stops, that is, the partial derivative of the voltage with respect to the displacement at the end point is equal to zero. Therefore, the boundary condition is the Newman boundary condition, and equation (3) is obtained:

[0041] Expand the voltage V on x using basis functions and apply the Galerkin method, i.e. the test basis function is the same as the expansion basis function, both of which are taken as N(x). Substituting into equation (3), we get The matrix equation form is , where A, B, C, and D are the matrix including the voltage term to be solved, the right-hand stimulation matrix, the coefficient matrix related to the time partial derivative, and the matrix related to the conductivity, respectively.

[0042] When at resting potential, there is no external stimulation and no current in the cell membrane. Therefore, in the form expressed by the matrix equation A, B, C, D, the external stimulation term B is 0 and the term multiplied by C is 0, that is, the matrix equation becomes AV re =D, that is, V re =A -1 D, the resting potential V can be obtained re size.

[0043] Get the resting potential V re Then, enter the time iteration phase.

[0044] The resulting matrix equation Perform time difference (time difference method), use the central difference format, and take the first square of the voltage as the midpoint of the two moments Partial derivative of voltage with respect to time becomes the time difference divided by the time step, that is, The numerator is the voltage difference between adjacent moments, and the denominator is the time step. The stimulus term B on the right is the average of the two moments. D is a constant term and remains unchanged. Substituting it into equation (4) yields equation (4).

[0045]

[0046] u n+1 Put one item on the left side of the equation, u n Moving to the right side of the equation, we get the time iteration form of the matrix equation

[0047]

[0048] The resting potential V obtained by the spectral element method without external stimulation is re As the initial value u1 of the time iteration, and through the time iteration form of the action potential conduction equation, substitute it into the u on the right side of equation (5)n In the equation, we get u n+1 is the voltage value u2 of each point on the neuron at the second moment, and then substitute the voltage value u2 into the right side of the equation u n In the equation, we get u n+1 is the voltage value u3 of each point on the neuron at the third moment... and so on, the voltage value at each moment can be obtained.

[0049] In the sampling of the actual model, due to the interval of sampling, the radius will jump between sampling points. In order to make the calculation results more consistent with the actual situation, the radius at the Gauss-Lobatto_Legendre (GLL) point inside the unit is interpolated, such as Figure 1 As shown. The specific operation method is to use linear or nonlinear interpolation methods based on the position of the GLL point in the unit. Using the radius a1 and a2 of the left and right end points of the unit, as well as the distance between the GLL points and the distance to the two end points, the radius a3 and a4 at the GLL point can be interpolated, thereby achieving a smoother change in the radius within the unit. A relatively advantageous feature of this method is that, by leveraging the advantages of the spectral element method, the voltage value at the multi-branch connection point can be simultaneously determined by the GLL points within the nearby units, which makes the result closer to the actual situation and has a smoother change process.

[0050] In the External Stimulation section, select the number of stimulation points and the location of the stimulation point. Additionally, you can adjust the values ​​of various parameters individually. For synaptic input, you can adjust the synaptic voltage to increase or decrease the stimulation. You can also adjust the duration of the stimulation, the time step dt, and various conductivities to simulate the stimulation process of various neurons. In the Output section, the default setting is to output the somatic cell feedback waveform over time.

[0051] The present invention plays a very important role in simulating the conduction process of action potentials in nerve fibers. It can greatly shorten the resources and time required for solving the problem, and can significantly improve the accuracy and authenticity. It plays a particularly important role in understanding the biological nervous system.

[0052] This method simulates the propagation of action potentials in neurons. The neuron model used is a text file that records morphological parameters such as the length of each neuron's branches, bifurcation connections, and radius. Input parameters are then set, such as the stimulation location, stimulation time, duration, stimulation intensity, spatial sampling interval, and temporal sampling interval.

[0053] A specific embodiment is given below.

[0054] Example 1:

[0055] Step one, determine the model and input parameters. Morphological model as shown in Figure 2 The applied stimulus is located at Figure 2 the red mark, the stimulation time is at 30 milliseconds, the spectral element method basis function is selected as 3 orders, dx takes 2um, and dt takes 0.02 milliseconds. Figure 3 The waveform diagram of the applied stimulus.

[0056] Step two, model initialization and establishment of conduction equation. According to the imported model, determine the connection mode between neuron branches, add materials to the divided cells according to the actual situation of the branches, including sodium ion conductivity, potassium ion conductivity, chloride ion conductivity, etc., and then according to the radius size of the cell end point and the number and position of Gauss-Lobatto_Legendre (GLL) points inside each cell, attach the radius size to the GLL points inside the cell to make a smooth transition in the cell radius. Then initialize the model, including the construction of each cell matrix, the establishment of the cell matrix, and finally the process of assembling the cell matrix into a global matrix. Since the ends of the neurons are not connected to other branches at this time, the boundary conditions here are uniform as the Neumann boundary condition, i.e. the natural boundary condition, which does not require additional processing of the matrix.

[0057] Step three, equation solving. The first main part of the solution is the solution of the resting potential, i.e. the stable membrane potential difference formed by the ion concentration difference itself without external stimulation. After setting the right end term to zero and obtaining the required resting potential, enter the time iteration link. We add the synaptic stimulus to the right end term, and at each time we judge whether the right end term at the current time is zero. When the right end term is not zero, i.e. at the time of stimulation, the right end term changes from zero to the stimulation term, and the action potential begins to propagate. Save the voltage change of each branch of the neuron with time in the final result, and draw the feedback waveform diagram of the body cell, which is the end point of all cells, i.e. the intersection of each main stem. The feedback obtained by the stimulus on the body cell is shown in Figure 4 .

[0058] To compare the algorithm and the accuracy of the results, calculate the error convergence curve of the real neuron, as shown in Figure 5 and 6 , the cell convergence and order convergence have reached a relatively accurate level.

[0059] In addition, in order to verify the accuracy of the simulation results, the results of the simulation of the neuron action potential using the finite difference method, which is widely used in the field of action potential simulation at this stage, are compared. Use the same neuron model as shown in Figure 2 , and add the same stimulus at the same position, i.e. the red mark in Figure 2 .

[0060] First, we compare the convergence of the two methods, and also increase the element density of the finite difference solution method, and compare it with its own error. The results are shown in Table 1. It can be seen that the convergence speed of the spectral element method is faster than that of the spectral element method. Figure 4 , which is significantly better than the finite difference method. Secondly, the spatial step size of the spectral element method is taken as 2 microns, and the spatial step size of the finite difference method is taken as 0.02 microns. The final simulation results of the two are compared, and the error size is 0.0019. In terms of the time used for solution, the solution time of the spectral element method is 30.677 seconds, while that of the finite difference method is 532.401 seconds. While ensuring that the errors of the two are very small, the method used in the present invention requires fewer grids, occupies less solution resources, and has a faster solution speed. It can be seen that the present invention is superior to the finite difference method of the same type in the modeling and solution method of simulating action potentials. It not only has a fast solution speed, but also high accuracy. Moreover, due to the different logic of the model and solution, the present invention applies the GLL integral point linear interpolation method to make the final result closer to the real situation.

[0061] Table 1

[0062] Reference cell number Error size 764 vs 15255 0.0668 1526 vs 15255 0.0387 3053 vs 15255 0.0144 6100 vs 15255 0.0075 15256 vs 15255 0.0029 30506 vs 15255 0.0013 61014 vs 15255 4.27e-4

[0063] In order to compare the simulation results of the two and reflect the calculation results in the figure, we took the bifurcated nerve fibers for comparison again. We also took the same nerve fibers and added stimulation at the same position. The final results are as follows: Figure 7 As shown, it can be seen that the results of the two at the same moment are in good agreement.

[0064] As the number of branches increases, the model tends to be complex, and the advantages of the present invention become more obvious. It can be seen that the present invention has advantages in solving both simple and complex neurons, and the more complex the model is and the more realistic it is, the more obvious the advantage is.

Claims

1. A method for simulating neuronal action potential conduction based on spectral element method, characterized in that The following steps are involved: 1) Neuron model input; 2) Conduction equation construction: Establish the partial differential equation for the propagation of action potential in a single nerve fiber; The specific steps of constructing the conduction equation are: (1) After setting the initial values ​​of the equation and the input waveform of the external stimulus, the stimulus time, position, duration, and stimulus intensity, each branch is constructed as a single nerve fiber model, and the electrical parameters and morphological parameters are introduced. Since the voltages at the connection points of adjacent branches are equal and the partial derivatives of the voltage with respect to time are equal, multiple single nerve fibers are connected together to form a neuron; (2) When the neuron is at resting potential, set the initial value of the conductivity at each point of the neuron; (3) Enter the equation initialization and construct a unit matrix for each solving unit, including the constant term in the equation, the differential terms of each order of voltage, and the stimulation term on the right side; after obtaining the required unit matrix, assemble each unit matrix together to form a global matrix, substitute the model parameters that change with position and construct the matrix equation to be solved; (4) The initial value of the solution to the matrix equation is the voltage at each point when there is no stimulation, that is, the value of the resting potential. Next, the spectral element matrix equation form of the equation is applied with time central difference to obtain the time iteration equation, that is, the left side of the equation is the related term of the voltage at the next moment, and the right side is the related term of the voltage at the current moment. By performing time iteration, the change of the voltage at each point of the neuron over time is obtained; 3) Solving the conduction equation: Using the neuron model input in step 1), the specific form of the conduction equation established in step 2) is modified, and the resting potential is solved using the spectral element method. To solve the conduction equation, the basic equation for action potential conduction is first constructed, which is a second-order partial differential equation of voltage with respect to time and position, and then the electrical and geometric parameters in the equation are set according to the imported external model; The spectral element method is used to solve the second-order partial differential equation. Specifically, the resting potential in the absence of stimulation is first solved. Then, the resting potential is used as the initial time value and combined with the time iteration method, the voltage at each point on the neuron is iteratively solved at all times to obtain the change of time. The voltage of the action potential is solved by the spectral element method, and the GLL point interpolation method inside the unit is applied to make the radius change smoother and closer to the actual situation. In assigning the radius of the GLL point inside the cell, after actual model sampling, the radius obtained by actual sampling is added to the two end points of the cell. According to the position distance between the GLL point inside the cell and the cell end points, the radius at the GLL point is obtained by linear interpolation, so that the radius has a smooth change inside the cell. 4) Obtain the change of voltage at each neuron point over time based on time iteration.

2. A method for simulating neuronal action potential conduction based on the spectral element method as claimed in claim 1, characterized in that In step 1), the neuron model input is to import the neuron morphological parameters, electrical parameters, and radius of each point to be simulated. The neuron morphological parameters include the length of each neuron branch and the connection mode between branches. The electrical parameters include the conductivity of sodium ions, potassium ions, and chloride ions at each location and the cytoplasmic coupling resistance.

Images