A simulation method for electric-thermal coupling of neuron action potential based on temperature effect

Abstract

The application discloses a neuron action potential electric-thermal coupling simulation method based on temperature effect, relates to the technical field of neuron electric activity simulation methods, and comprises the following steps: S1, geometric modeling and mesh division; S2, temperature-dependent improvement is carried out on a Hodgkin-Huxley model; S3, material assignment: material properties are assigned to each structure of the model in Abaqus; S4, a user-defined subroutine is called to realize the regulation of temperature as a continuous variable parameter on the dynamics of a Hodgkin-Huxley ion channel; the user-defined subroutine comprises a USDFLD subroutine, a UMATHT subroutine and a HETVAL subroutine; a neuron axon three-dimensional finite element model is established through the Hodgkin-Huxley model, the temperature-dependent Hodgkin-Huxley dynamics equation is directly embedded through the user-defined subroutine, and electric-thermal-mechanical three-field full coupling is realized; meanwhile, the application also provides an effective simulation tool and a theoretical basis for the research on the influence of a subsequent thermal environment on neural functions.

Description

Technical Field

[0001] This invention relates to the field of neuronal electrical activity simulation methods, and in particular to a neuronal action potential electrothermal coupling simulation method based on temperature effects. Background Technology

[0002] Neuronal electrical activity is fundamental to the function of the nervous system, and action potentials are one of the most basic characteristics of nerve cells. Temperature has a crucial impact on neuronal electrical activity, causing changes in the frequency of neuronal action potentials. Even small temperature fluctuations can affect the activation and inactivation of multiple ion channels, thus influencing the generation of action potentials. This effect is particularly pronounced in extreme high-temperature environments, potentially leading to neuronal overexcitation or functional inhibition, ultimately resulting in severe neurological dysfunction. Furthermore, temperature has been found to affect various biophysical processes, including biomolecular diffusion, enzyme activity, and heat shock-induced gene expression. Therefore, research and simulation of temperature-sensitive neuronal electrical activity based on the classic Heckscher-Hirschman (HH) model have significant theoretical and clinical implications for understanding the mechanisms of neurological dysfunction under high-temperature conditions.

[0003] Methods for exploring the effect of temperature on neuronal excitability using one-dimensional cable theory or simplified HH models are becoming increasingly sophisticated. These studies often employ analytical methods, neglecting the three-dimensional geometric features of nerve cells. While a few three-dimensional finite element models can handle complex geometries or electrothermal equivalent couplings, they are mostly temperature-fixed, studying neuronal action potential changes through electrical, optical, or chemical stimulation. Systematic studies using custom subroutines to implement temperature-sensitive ion channel gating variables and multi-physics coupling are rare. Furthermore, in three-dimensional finite element modeling, the thin-layer approximation of neuronal axons is widely used in neuronal electrophysiological finite element modeling due to its computational simplicity. However, this approximation neglects the radial charge gradient, electrostriction, and piezoelectric effects within the membrane layer, making it difficult to accurately realize action potential-induced membrane thickness changes and mechanical stress. Therefore, this paper proposes a simulation method for electrothermal coupling of neuronal action potentials based on temperature effects. Summary of the Invention

[0004] In view of this, this invention addresses the shortcomings of existing technologies by providing a temperature-based simulation method for the electrothermal coupling of neuronal action potentials. It establishes a three-dimensional finite element model of neuronal axons using the HH biophysical model and directly embeds the temperature-dependent Hodgkin-Huxley dynamic equations through a user-defined subroutine. By increasing the ambient temperature from 37°C to 45°C, the method predicts the generation and propagation of axonal action potentials at different high temperatures, achieving full coupling of the electro-thermal-mechanical three fields. This also provides an effective simulation tool and theoretical basis for subsequent research on the effects of thermal environments on neural function.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A simulation method for electrothermal coupling of neuronal action potentials based on temperature effects includes the following steps:

[0007] S1. Geometric Modeling and Mesh Generation: A finite element model of neuronal axons is constructed based on the Hodgkin-Huxley model. The finite element model of neuronal axons includes a two-dimensional axisymmetric model and a three-dimensional finite element model. After the finite element model of neuronal axons is constructed, mesh generation is performed.

[0008] S2. Temperature-dependent improvements to the Hodgkin-Huxley model: Introducing Q... 10 Temperature correction factor affects the gating rate constant α in the Hodgkin-Huxley model. x and β x Scaling was performed to establish temperature-dependent ion channel kinetic equations;

[0009] S3. Material Assignment: Assign material properties to each structure in the model in Abaqus;

[0010] S4. User-defined subroutine settings: Call user-defined subroutines to regulate the dynamics of Hodgkin-Huxley ion channels as a continuously variable parameter; these user-defined subroutines include USDFLD, UMATHT, and HETVAL subroutines.

[0011] As a preferred embodiment, the core membrane potential dynamics equation of the Hodgkin-Huxley model in step S2 is:

[0012] ;

[0013] The evolution of gated variables satisfies:

[0014] ;

[0015] in, V is the stimulation current; C is the membrane potential; V is the membrane capacitance. , and It is the equilibrium potential of each ion; , α and β are the maximum conductivities of each ion; m and h are the sodium ion channel gating variables; n is the potassium ion channel gating variable; α and β are rate functions related to membrane potential, and the specific equations are as follows:

[0016] ;

[0017] ;

[0018] ;

[0019] ;

[0020] ;

[0021] ;

[0022] Introducing Q 10 Temperature correction factor on gating rate constant and Scaling is performed, the Q 10 The expression for the temperature correction factor is:

[0023] ;

[0024] Where: T is the current temperature, and the reference temperature is 6.3℃.

[0025] As a preferred embodiment: in step S4, the USDFLD subroutine uses Q... 10 The rate constant of the gated variable is corrected in real time by the temperature coefficient; and the Na, K and L ion currents are calculated based on the equivalent membrane potential of the current temperature field, and the calculation results are stored in the state variables.

[0026] As a preferred embodiment, the specific steps of the USDFLD subroutine in step S4 are as follows: First, write the USDFLD subroutine in Fortran; second, calculate Q. 10 The correction factor is used to calculate the evolution of the gated variable, calculate the ion current, and finally use the calculated gated variable or current as the field variable for the material properties to depend on.

[0027] As a preferred embodiment: in step S4, the UMATHT subroutine defines heat transfer properties, including specific heat capacity and thermal conductivity updated with state variables, thereby realizing dynamic updating of electrical conductivity.

[0028] As a preferred option, the specific steps of the UMATHT subroutine in step S4 are as follows: first, write the UMATHT subroutine interface; second, input the current temperature, temperature gradient, and state variables; and finally, output the heat flux and associate it with the thermal conductivity in the material definition.

[0029] As a preferred embodiment: in step S4, the HETVAL subroutine applies an external stimulation current and calculates the equivalent heat flux generated by the ion current.

[0030] As a preferred embodiment, the specific steps of the HETVAL subroutine in step S4 are as follows: first, write the HETVAL subroutine interface; second, input the current temperature, time, and state variables; and finally, calculate the equivalent heat generation. The HETVAL subroutine directly converts the calculated temperature-dependent ion current into the force driving the evolution of the temperature field, forming a closed-loop feedback of temperature-gated-current-temperature.

[0031] As a preferred embodiment, the meshing in step S1 is specifically as follows: in the two-dimensional axisymmetric model, axisymmetric quadrilateral elements are used for structured meshing; in the three-dimensional finite element model, hexahedral dominant meshes are used, combined with wedge elements to handle the transition region, and for the cell membrane layer with a thickness of only 3 nm, the minimum mesh size is set to be less than 3 nm to fully analyze its physical properties.

[0032] As a preferred embodiment, the temperature-effect-based neuronal action potential electrothermal coupling simulation method further includes step S5: verifying HH dynamics using a two-dimensional axisymmetric model; verifying electrothermal equivalence using a three-dimensional finite element model; specifically, the electrothermal equivalence verification involves applying a voltage boundary to the outer surface of the extracellular medium, and in the axisymmetric verification, verifying the correctness of the electrothermal equivalence by comparing the membrane potential decay law along the axial direction with the analytical solution of the passive cable equation; specifically, the HH dynamics verification involves verifying the correctness of embedding temperature-dependent ion channel dynamics in the custom subroutine by using the action potential waveform, gating variable, and ion current calculated by Abaqus after applying a stimulation current;

[0033] Compared with the prior art, the present invention has significant advantages and beneficial effects. Specifically, as can be seen from the above technical solution, this application is based on the classic Hodgkin-Huxley model, combined with Q... 10 A temperature correction mechanism was established using a three-dimensional finite element model of a temperature-sensitive neuronal axon with finite membrane thickness in the Abaqus platform, and temperature-dependent ion channel dynamics were embedded through a subroutine. The model showed good agreement with the classical HH model and original experimental results at 6.3℃, demonstrating the effectiveness of the established model. The established model can effectively characterize the electrophysiological response of neuronal axons under high-temperature conditions, providing a theoretical basis and simulation tool for thermal injury risk assessment, research on the mechanisms of abnormal neuronal function under high-temperature conditions, and optimization of neurothermotherapy protocols.

[0034] To more clearly illustrate the structural features and effects of the present invention, a detailed description is provided below in conjunction with the accompanying drawings and specific embodiments. Attached Figure Description

[0035] Figure 1 This is a schematic diagram of the neuron axon model of the present invention;

[0036] Figure 2 This is a schematic diagram illustrating the electrothermal equivalence verification of the model of the present invention;

[0037] Figure 3 This is a schematic diagram of the HH kinetic verification results of the present invention;

[0038] Figure 4 This is a schematic diagram illustrating the waveform changes of neuronal action potentials at different temperatures according to the present invention.

[0039] Figure 5 This is a schematic diagram illustrating the dynamic changes of neuronal ion channel gating variables (m, h, n) at different temperatures according to the present invention.

[0040] Figure 6 This is a schematic diagram illustrating the changes in neuronal ionic conductivity at different temperatures according to the present invention.

[0041] Figure 7 This is a schematic diagram illustrating the dynamic changes of neuronal ion currents at different temperatures according to the present invention.

[0042] Explanation of reference numerals in the attached diagram:

[0043] Figure 2 In the diagram: a is a schematic diagram of the spatial variation of the membrane voltage; b is a schematic diagram of the axial voltage distribution of the membrane layer.

[0044] Figure 3 In the diagram: a represents the action potential time waveform; b represents the curves of the gating variables (m, h, n) over time; c represents the changes in Na and K ion conductivities over time; d represents the changes in Na, K, and L ion currents over time.

[0045] Figure 5 In the figure: the black curve represents the sodium channel activation variable m, the red curve represents the potassium channel activation variable n, and the blue curve represents the sodium channel inactivation variable h; in the figure, a, b, c, d, and e correspond to the gating variable responses at 37℃, 39℃, 41℃, 43℃, and 45℃, respectively.

[0046] Figure 6 In the figure: the black curve represents sodium conductivity (gNa), and the red curve represents potassium conductivity (gK); a, b, c, d, and e in the figure correspond to the gated variable responses at 37℃, 39℃, 41℃, 43℃, and 45℃, respectively.

[0047] Figure 7 In the diagram: the black curve represents the inward sodium current INa, the red curve represents the outward potassium current IK, the blue curve represents the leakage current IL, and the green curve represents the capacitance current IC; in the diagram, a, b, c, d, and e correspond to the schematic diagrams of operating conditions at 37℃, 39℃, 41℃, 43℃, and 45℃, respectively. Detailed Implementation

[0048] The present invention is as follows Figure 1As shown in Figure 7, a simulation method for electrothermal coupling of neuronal action potentials based on temperature effects includes the following steps:

[0049] S1. Geometric Modeling and Mesh Generation: A finite element model of neuronal axons is constructed based on the Hodgkin-Huxley model. The finite element model of neuronal axons includes a two-dimensional axisymmetric model and a three-dimensional finite element model. After the finite element model of neuronal axons is constructed, mesh generation is performed.

[0050] S2. Temperature-dependent improvements to the Hodgkin-Huxley model: Introducing Q... 10 Temperature correction factor affects the gating rate constant α in the Hodgkin-Huxley model. x and β x Scaling was performed to establish temperature-dependent ion channel kinetic equations;

[0051] S3. Material Assignment: Assign material properties to each structure in the model in Abaqus;

[0052] S4. User-defined subroutine settings: Call user-defined subroutines to regulate the dynamics of Hodgkin-Huxley ion channels as a continuously variable parameter; these user-defined subroutines include USDFLD, UMATHT, and HETVAL subroutines.

[0053] The core membrane potential dynamics equation of the Hodgkin-Huxley model in step S2 is:

[0054] ;

[0055] The evolution of gated variables satisfies:

[0056] ;

[0057] in, V is the stimulation current; C is the membrane potential; V is the membrane capacitance. , and It is the equilibrium potential of each ion; , α and β are the maximum conductivities of each ion; m and h are the sodium ion channel gating variables; n is the potassium ion channel gating variable; α and β are rate functions related to membrane potential, and the specific equations are as follows:

[0058] ;

[0059] ;

[0060] ;

[0061] ;

[0062] ;

[0063] ;

[0064] Introducing Q 10 Temperature correction factor on gating rate constant and Scaling, the Q 10 The expression for the temperature correction factor is:

[0065] ;

[0066] Where: T is the current temperature, and the reference temperature is 6.3℃.

[0067] In step S4, the USDFLD subroutine uses Q. 10 The rate constant of the gated variable is corrected in real time by the temperature coefficient; and the Na, K and L ion currents are calculated based on the equivalent membrane potential of the current temperature field, and the calculation results are stored in the state variables.

[0068] The specific steps of the USDFLD subroutine in step S4 are as follows: First, write the USDFLD subroutine in Fortran; second, calculate Q. 10 The correction factor is used to calculate the evolution of the gated variable, calculate the ion current, and finally use the calculated gated variable or current as the field variable for the material properties to depend on.

[0069] In step S4, the UMATHT subroutine defines heat transfer properties, including specific heat capacity and thermal conductivity updated with state variables, thereby achieving dynamic updates of electrical conductivity.

[0070] The specific steps of the UMATHT subroutine in step S4 are as follows: First, write the UMATHT subroutine interface; second, input the current temperature, temperature gradient, and state variables; and finally, output the heat flux and associate it with the thermal conductivity in the material definition.

[0071] In step S4, the HETVAL subroutine applies an external stimulation current and calculates the equivalent heat flux generated by the ion current.

[0072] The specific steps of the HETVAL subroutine in step S4 are as follows: First, write the HETVAL subroutine interface; second, input the current temperature, time, and state variables; and finally, calculate the equivalent heat generation. The HETVAL subroutine directly converts the calculated temperature-dependent ion current into the force driving the evolution of the temperature field, forming a closed-loop feedback of temperature-gated-current-temperature.

[0073] The mesh generation in step S1 is as follows: In the two-dimensional axisymmetric model, axisymmetric quadrilateral elements are used for structured mesh generation; in the three-dimensional finite element model, hexahedral dominant mesh is used, combined with wedge elements to handle the transition region. For the cell membrane layer with a thickness of only 3 nm, the minimum mesh size is set to be less than 3 nm to fully analyze its physical properties.

[0074] The temperature-effect-based neuronal action potential electrothermal coupling simulation method also includes step S5: verifying HH dynamics using a two-dimensional axisymmetric model; verifying electrothermal equivalence using a three-dimensional finite element model; specifically, the electrothermal equivalence verification involves applying a voltage boundary to the outer surface of the extracellular medium, and in the axisymmetric verification, verifying the correctness of the electrothermal equivalence by comparing the membrane potential decay law along the axial direction with the analytical solution of the passive cable equation; specifically, the HH dynamics verification involves verifying the correctness of embedding temperature-dependent ion channel dynamics in the custom subroutine by using the action potential waveform, gating variable, and ion current calculated by Abaqus after applying a stimulation current;

[0075] Example: A simulation method for electrothermal coupling of neuronal action potentials based on temperature effects

[0076] Geometric modeling and mesh generation:

[0077] Based on research into the physiological structure of single-neuron axons, which are characterized by their elongated shape, smooth surface, and few branches, this application simplifies the axon geometry to a standard three-layer coaxial cylinder, consisting of intracellular mediator, cell membrane, and extracellular mediator from the inside out. Specific geometric dimensions are shown in Table 1. All simulations were performed using Abaqus CAE 2020 software. The specific modeling process is as follows:

[0078] First, a two-dimensional axisymmetric analysis platform was used to draw a two-dimensional sketch. Strictly adhering to the geometric dimensional parameters in Table 1, the radial range and length of the three-layer axons were defined. After completing the sketch dimensional constraints, a three-dimensional solid model was generated from the two-dimensional cross-section through rotation operations to ensure the coaxiality and symmetry of each layer. Since the motion signal propagates longitudinally in the axons, the two-dimensional axisymmetric model can significantly reduce computational costs while maintaining computational accuracy; therefore, it was used for HH dynamics verification. The three-dimensional model was used for electrothermal equivalence verification to recreate the physical field distribution in real three-dimensional space.

[0079] After the model was completed, meshing was performed. In the two-dimensional axisymmetric model, axisymmetric quadrilateral elements were used for structured meshing, and the final model contained 15,652 nodes and 15,300 elements. In the three-dimensional model, hexahedral dominant mesh was used, combined with wedge elements to handle the transition region. For the cell membrane layer with a thickness of only 3 nm, the minimum mesh size was set to be less than 3 nm to fully analyze its physical properties. The final three-dimensional model contained 1,276,922 nodes and 1,426,524 elements, which ensured both computational accuracy and computational efficiency (see Figure 1).

[0080] Table 1: Geometric Model Parameters

[0081]

[0082] HH model and temperature-dependent improvements

[0083] The HH model explains the electrical activity of neurons at the ionic level, dividing the total transmembrane current into capacitive current and three types of ionic currents (sodium current, potassium current, and leakage current). It uses three gating variables (m, h, n) to describe the activation and deactivation dynamics of voltage-gated ion channels. Its core membrane potential dynamic equation is:

[0084] ; (1)

[0085] The evolution of gated variables satisfies:

[0086] ; (2)

[0087] in, V is the stimulation current; C is the membrane potential; V is the membrane capacitance. , and It is the equilibrium potential of each ion; , α and β are the maximum conductivities of each ion; m and h are the sodium ion channel gating variables; n is the potassium ion channel gating variable; α and β are rate functions related to membrane potential, and the specific equations are as follows:

[0088] ; (3)

[0089] ; (4)

[0090] ; (5)

[0091] ; (6)

[0092] ; (7)

[0093] ; (8)

[0094] The parameters of the classic HH model are based on experimental data from the voltage clamp test of giant squid axons at 6.3℃ (see Table 2). Since the original model cannot accurately simulate action potential conduction above 30℃, this application introduces Q... 10 Temperature correction factor on gating rate constant and Scaling is applied to adjust the temperature range from 37℃ to 45℃. This temperature-dependent improvement forms the core foundation for subsequent subroutine implementations.

[0095] ; (9)

[0096] Table 2: HH Model Parameter Values

[0097]

[0098] Material Assignment: Each structure in the model was assigned its own material properties in Abaqus CAE 6.13-3 (see Table 3). To verify the electrothermal equivalence of the model, the film capacitance and conductivity were equivalent to specific heat capacity and thermal conductivity, respectively. The coefficient of thermal expansion was dynamically coupled to the voltage field through a subroutine, thus transforming the electrophysiological process into an equivalent heat transfer problem. This equivalent method provides a multiphysics solution framework for embedding temperature-sensitive HH dynamics into user subroutines.

[0099] Table 3: Material properties of each layer of the neuron

[0100]

[0101] User subroutine implementation:

[0102] Based on the HH model equations (Equations (1)–(3)) and the parameters in Table 2, this application developed and invoked three user-defined subroutines (USDFLD, UMATHT, and HETVAL) to achieve precise control of Hodgkin-Huxley ion channel dynamics by using temperature as a continuously variable parameter. Specifically, the USDFLD subroutine corrects the rate constants of the gate variables (m, h, n) in real time using the Q10 temperature coefficient, calculates Na, K, and L ion currents based on the equivalent membrane potential of the current temperature field (TEMP), and stores the calculation results in the state variables; the UMATHT subroutine defines the heat transfer properties (specific heat capacity and thermal conductivity updated with the state variables); and the HETVAL subroutine applies an external stimulation current and calculates the equivalent heat flux generated by the ion current. The three subroutines collaboratively solve the membrane potential evolution equation shown in Equation (1), supporting the simulation of action potential propagation along the axon longitudinal direction and reserving an interface for electro-thermal-mechanical three-field coupling.

[0103] Implementation principle: The core is the HH equation.

[0104] ;

[0105] Where current = conductance × voltage difference, conductance is related to the gating variables, and the dynamics of the gating variables (m, h, n) are regulated by voltage and temperature. Therefore, changing the temperature can change the current, thereby changing the membrane potential.

[0106] In Abaqus, electrical equations cannot be solved directly, but thermal analysis can simulate the diffusion process. Therefore, by mapping electrical quantities to thermal quantities through electrothermal equivalence, the electrical equations can be solved in Abaqus.

[0107] Specific implementation steps:

[0108] USDFLD Subroutine - Real-time Calculation of Temperature-Dependent Gated Variables and Ion Current

[0109] After modeling is completed, material properties are assigned to it. The USDFLD subroutine is used to declare that some of the defined material properties (electrical conductivity, thermal conductivity, coefficient of thermal expansion) change with temperature.

[0110] Implementation principle: It reads the current temperature field (TEMP) and applies Q... 10 The gating rate constants (α, β) are corrected, the evolution of the gating variables (m, h, n) is calculated, and then Na⁺, K⁺, and leakage current are calculated. The results are stored in the state variables for use by other subroutines.

[0111] Steps: First, write the USDFLD subroutine in Fortran, then calculate Q. 10 The correction factor is used to calculate the evolution of the gated variable and the ion current. Finally, the calculated gated variable or current is used as a field variable for material property dependence. USDFLD makes temperature a continuously variable parameter, correcting the rate constant in real time at each time increment and each integration point, thus achieving Q... 10 Precise temperature-dependent control of the dynamics of the gated variables (m, h, n). At higher temperatures, the rate increases → the action potential duration shortens and repolarization accelerates; simultaneously, the conductance / current amplitude is adjusted → the peak value decreases.

[0112] UMATHT subroutine – Define / update heat transfer properties (dynamic material behavior)

[0113] Implementation principle: UMATHT is used to define custom heat transfer constitutive behavior (specific heat capacity, thermal conductivity). It is used to define heat transfer properties (specific heat capacity and thermal conductivity updated with state variables) to achieve dynamic updates of electrical conductivity.

[0114] The specific steps involve first writing the UMATHT subroutine interface, then inputting the current temperature, temperature gradient, and state variables (gated variables / currents passed from USDFLD), and finally outputting the heat flux and associating it with thermal conductivity in the material definition. This subroutine, together with USDFLD, implements the dependence of material properties on temperature / gated state.

[0115] HETVAL subroutine – Applying an equivalent heat source

[0116] Implementation principle: HETVAL defines the internal heat generation rate, which corresponds to the driving force of ion current generation. The HETVAL subroutine applies an external stimulation current and calculates the equivalent heat flux generated by the ion current.

[0117] Steps: First, write the HETVAL subroutine interface; second, input the current temperature, time, and state variables (ion current values ​​obtained from USDFLD); and finally, calculate the equivalent heat generation.

[0118] HETVAL directly converts the calculated temperature-dependent ion current into the force driving the evolution of the temperature field (i.e., membrane potential), thus closing the entire system loop: Temperature → USDFLD update gating → Calculated current → HETVAL drive → New temperature field

[0119] Model validity verification

[0120] To verify the effectiveness of the constructed neuronal axon finite element model, the axon model was validated by both the electrothermal equivalent model and the classical HH model (both were performed at the reference temperature of 6.3℃).

[0121] (1) Verification of electrothermal equivalence: A voltage boundary is applied to the outer surface of the extracellular medium (Gaussian distribution is used in the axial direction to simulate local electrode stimulation). A circumferential mode can be further introduced into the three-dimensional model to examine the potential distribution under asymmetric stimulation. In the axisymmetric verification, the correctness of the electrothermal equivalence is verified by comparing the membrane potential decay law along the axial direction with the analytical solution of the passive cable equation.

[0122] (2) HH kinetic verification: applied size 10 After a stimulation current lasting 0.5 ms, the action potential waveform, gating variable, and ion current calculated by Abaqus were compared with the theoretical predictions from directly solving the classical HH equation and the original data from Hodgkin-Huxley, verifying the correctness of embedding the HH ion channel dynamics into the custom subroutine.

[0123] Discussion of Results:

[0124] Model validity results verification:

[0125] Electrothermal equivalence verification:

[0126] The model was validated for electrothermal equivalence at a reference temperature of 6.3℃. Figure 2 (a) The spatial distribution cloud map of the membrane voltage in the entire three-dimensional axon model after applying voltage boundary conditions is presented, which intuitively shows the continuous propagation characteristics of the potential along the axon length direction; Figure 2 (b) A quantitative comparison of the membrane potential decay curves along the axial direction was then performed, where the solid line represents the Abaqus finite element simulation results and the dashed line represents the analytical solution of the passive cable equation. The results show that the two are highly consistent.

[0127] HH kinetic verification

[0128] The model was validated for HH kinetics at a reference temperature of 6.3℃. Figure 3 (a) shows the time waveform of the action potential. Figure 3 (b) Presents the curves of the changes in sodium activation-gated variable m, potassium activation-gated variable n, and sodium inactivation-gated variable h over time. Figure 3 (c) shows the dynamic response of sodium conductivity and potassium conductivity. Figure 3 (d) shows the variation of Na, K, and L ion currents, as well as the action potential waveform, gating variables m, h, and n variation curves, and current changes obtained from Abaqus simulation.

[0129] The results show that the peak value, peak width, rise and fall slopes of the action potential, as well as the steady-state values ​​and arrival times of the gating variables obtained from the Abaqus simulation, are in high agreement with the theoretical predictions of the classical HH model and the original experimental data from Hodgkin-Huxley. The peak amplitude, time history, and polarity of each ion current are also consistent with the theoretical values. This verification result fully demonstrates that the user subroutines such as USDFLD developed in this application can correctly realize the temperature-dependent HH ion channel dynamics, providing a reliable guarantee for subsequent simulation analysis under high-temperature conditions.

[0130] Changes in neuronal electrical activity at different temperatures:

[0131] The effect of temperature on membrane potential:

[0132] Figure 4 The dynamic changes of neuronal action potentials are shown within a temperature range of 37℃ to 45℃. As can be seen from the figure, the peak amplitude of the action potential progressively decreases as the temperature gradually increases from 37℃ to 45℃. At higher temperatures (close to 45℃), the peak potential decreases significantly, with some simulation results showing the peak value dropping from approximately 105 mV at 37℃ to even lower levels. Simultaneously, the action potential duration shortens significantly, particularly with accelerated repolarization, resulting in a narrower and sharper waveform.

[0133] Within the intermediate temperature range (around 39°C), the frequency of action potential firing increases slightly. However, as the temperature rises further, excitability turns into inhibition. At 45°C, the action potential waveform tends to simplify, with a significant decrease in peak value, even exhibiting a near-plateau pattern. These characteristics are consistent with experimental observations in the literature showing reduced action potential amplitude and shortened duration under high-temperature conditions. Miller et al. further confirmed through thermodynamic models that increased temperature accelerates ion channel gating dynamics and reduces channel availability, directly leading to a decrease in action potential peak value and accelerated repolarization. These changes reflect the combined effect of high temperature on action potential morphology by accelerating ion channel gating dynamics while simultaneously reducing channel availability and maximum conductivity.

[0134] Overall, within the temperature range of 37℃–45℃, action potentials exhibit characteristics of decreasing amplitude, shortening duration, and simplified waveform as temperature increases, which may have a significant impact on nerve signal transmission and the stability of cardiac rhythm.

[0135] The effect of temperature on gated variables:

[0136] To further elucidate the intrinsic mechanism of temperature-regulated action potential morphology, this application analyzed the kinetic response of ion channel gating variables at different temperatures, and the results are as follows: Figure 3 As shown, temperature changes can modulate the activation and deactivation dynamics of sodium channels, affecting the amplitude and duration of action potentials. Similarly, temperature fluctuations affect potassium channel activity, altering the repolarization phase of action potentials and thus impacting the efficiency of neural signal transmission. As the temperature increases from 37°C, the dynamics of the sodium channel activation-gated variable (m), deactivation-gated variable (h), and potassium channel activation-gated variable (n) all accelerate significantly. The rise phase and peak time of the m variable appear significantly earlier, while the peak amplitude decreases slightly at higher temperatures; the deactivation process of the h variable accelerates, resulting in a steeper fall phase and a slower recovery process. At the highest temperature (close to 45°C), the oscillation amplitudes of the three gating variables gradually decrease, and the waveforms tend to flatten, especially in the recovery phase of the later action potential. This indicates that high temperatures significantly accelerate the activation and deactivation dynamics of the gating variables, while reducing their dynamic range, ultimately leading to a narrower and sharper trend in the action potential.

[0137] The effect of temperature on ionic conductivity:

[0138] Figure 4 The response characteristics of neuronal sodium and potassium ion conductance at different temperatures are described. Increased temperature leads to a higher peak sodium conductance (…). ) and peak potassium conductivity ( Both exhibited a progressive decreasing trend. The peak sodium conductivity decreased rapidly from a relatively high level at lower temperatures, and at higher temperatures, the peak value decreased significantly and the duration shortened markedly. Potassium conductivity also showed a decrease in peak amplitude and a relatively shortened activation delay. At around 45℃, the peak values ​​of both major ion conductivities decreased significantly, and the waveforms became flat, reflecting that the maximum conductivity of voltage-gated ion channels was suppressed under high-temperature conditions.

[0139] The effect of temperature on ion current:

[0140] like Figure 7 As shown, the inward sodium current (I Na ) and outward potassium current (I K The temperature dependence of ion current changes significantly with increasing temperature. Combining the changes in the gating variables (Figure 5) and ion conductance (Figure 6) described earlier, the temperature dependence of ion current is fully explained by the physical mechanism. As the peak amplitude of sodium current progressively decreases, both its activation and deactivation processes accelerate, resulting in a narrower and sharper current peak. At higher temperatures, the residual component of the later stages of sodium current also weakens significantly. The peak value of potassium current also decreases, and the activation kinetics accelerate, but the overall amplitude decay is more significant. Under the highest temperature conditions, the oscillation amplitude of both ion currents decreases significantly, and in some records, the current almost plateaus or approaches the baseline level, indicating that high temperatures strongly suppress the intensity of transmembrane ion flow.

[0141] Overall, during the temperature increase from 37℃ to 45℃, the kinetics of the gating variable accelerated significantly, the maximum value of ionic conductance decreased, while the peak amplitude of ionic current progressively decreased and its duration shortened. These changes collectively reflect the combined effects of temperature on the conformational kinetics of ion channel proteins, channel opening probability, and conductivity, ultimately leading to a decrease in action potential amplitude and a shortening of its duration. This phenomenon is related to the accelerated channel gating rate at high temperatures (…). (Effect) But it is closely related to the mechanism of reducing channel availability and single-channel conductance, and has an important inhibitory effect on cell excitability and signal transduction under high temperature environment.

[0142] The effects of temperature on neuronal electrical activity are multi-layered. This application, using a three-dimensional finite element model within an electrothermal equivalence framework, systematically reveals the variations in action potential, gating variables, ionic conductance, and ionic current within the temperature range of 37℃–45℃. Temperature primarily affects neuronal electrical activity through Q... 10 Temperature coefficient modulates ion channel kinetics. In this model, Q is introduced. 10 The factor (the gating rate constant is usually taken as Q) 10 ≈3, maximum conductance is taken as Q 10Temperature corrections to the rate constants (α and β) by approximately 1–1.5 result in a synchronous acceleration of all gating processes (rapid activation of m, slow inactivation of h, and activation of n). This is consistent with the original Hodgkin-Huxley model and subsequent temperature-dependent studies: high temperatures make action potentials "narrow and sharp," accelerating repolarization and thus shortening the action potential duration and theoretically increasing neural conduction speed.

[0143] However, the decrease in peak action potential amplitude and excitatory inhibition observed at higher temperatures (>41°C) cannot be explained solely by gating acceleration. Our results show that high temperatures simultaneously lead to a decrease in peak ionic conductance and a reduction in ionic current intensity, which is related to decreased conformational stability of channel proteins, reduced single-channel conductance, and decreased channel availability. Increased temperature can cause a mild change in maximum channel conductance (Q0). 10 (1.2–1.5), but when a certain threshold is exceeded, the competitive effect between channel dynamics and conductance leads to an initial increase followed by a decrease in excitability, and even conduction block. The near-plateau waveform at 45°C in this model is highly consistent with the manifestations of abnormal nerve function under physiological conditions such as heatstroke or high fever, suggesting that high temperature may impair signal transduction by inhibiting sodium current-dominated depolarization.

[0144] Compared with the traditional one-dimensional HH model, the electro-thermal-mechanical three-field coupling model implemented in Abaqus in this application has significant advantages: it can not only handle the three-dimensional geometric features of nerve axons and the radial gradient of the thin film layer (3 nm thick), but can also be further extended to mechanical stress analysis (such as action potential-induced membrane thickness changes or electrostriction effects). This provides a simulation platform that is closer to physiology for exploring the mechanism of thermal damage.

[0145] The results show that within the temperature range of 37℃–45℃, increased temperature significantly affects neuronal electrical activity. With increasing temperature, the action potential duration shortens significantly, the peak amplitude gradually decreases, and the waveform tends to simplify; the activation and deactivation mechanics of gating variables m, h, and n are significantly accelerated; and the peak conductance of sodium and potassium ions and the corresponding ion current amplitude both show a decreasing trend. This indicates that high temperature, on the one hand, affects the neuronal electrical activity through Q... 10 The effect accelerates the ion channel gating process, but on the other hand, it reduces channel availability and transmembrane ion flow intensity, ultimately causing neuronal excitability to gradually shift from moderate enhancement to inhibition.

[0146] The key design focus of this invention is:

[0147] This application is based on the classic HH model, combined with Q 10A temperature correction mechanism was established using a three-dimensional finite element model of a temperature-sensitive neuronal axon with finite membrane thickness in the Abaqus platform, and temperature-dependent ion channel dynamics were embedded through a subroutine. The model showed good agreement with the classical HH model and original experimental results at 6.3℃, demonstrating the effectiveness of the established model. The established model can effectively characterize the electrophysiological response of neuronal axons under high-temperature conditions, providing a theoretical basis and simulation tool for thermal injury risk assessment, research on the mechanisms of abnormal neuronal function under high-temperature conditions, and optimization of neurothermotherapy protocols.

[0148] The above description is merely a preferred embodiment of the present invention and does not constitute any limitation on the technical scope of the present invention. Therefore, any minor modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention shall still fall within the scope of the technical solution of the present invention.

Claims

1. A simulation method for electrothermal coupling of neuronal action potentials based on temperature effects, characterized in that: Includes the following steps: S1. Geometric Modeling and Mesh Generation: A finite element model of neuronal axons is constructed based on the Hodgkin-Huxley model. The finite element model of neuronal axons includes a two-dimensional axisymmetric model and a three-dimensional finite element model. After the finite element model of neuronal axons is constructed, mesh generation is performed. S2. Temperature-dependent improvements to the Hodgkin-Huxley model: Introducing Q... 10 Temperature correction factor affects the gating rate constant α in the Hodgkin-Huxley model. x and β x Scaling was performed to establish temperature-dependent ion channel kinetic equations; S3. Material Assignment: Assign material properties to each structure in the model in Abaqus; S4. User-defined subroutine settings: Call user-defined subroutines to regulate the dynamics of Hodgkin-Huxley ion channels as a continuously variable parameter; these user-defined subroutines include USDFLD, UMATHT, and HETVAL subroutines.

2. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: The core membrane potential dynamics equation of the Hodgkin-Huxley model in step S2 is: ； The evolution of gated variables satisfies: ； in, V is the stimulation current; C is the membrane potential; V is the membrane capacitance. , and It is the equilibrium potential of each ion; , α and β are the maximum conductivities of each ion; m and h are the sodium ion channel gating variables; n is the potassium ion channel gating variable; α and β are rate functions related to membrane potential, and the specific equations are as follows: ； ； ； ； ； ； Introducing Q 10 Temperature correction factor on gating rate constant and Scaling is performed, the Q 10 The expression for the temperature correction factor is: ； Where: T is the current temperature, and the reference temperature is 6.3℃.

3. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: In step S4, the USDFLD subroutine uses Q. 10 The rate constant of the temperature coefficient-corrected gating variable in real time; The Na, K, and L ion currents are calculated based on the equivalent membrane potential of the current temperature field, and the calculation results are stored in the state variables.

4. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 3, characterized in that: The specific steps of the USDFLD subroutine in step S4 are as follows: First, write the USDFLD subroutine in Fortran; second, calculate Q. 10 The correction factor is used to calculate the evolution of the gated variable, calculate the ion current, and finally use the calculated gated variable or current as the field variable for the material properties to depend on.

5. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: In step S4, the UMATHT subroutine defines heat transfer properties, including specific heat capacity and thermal conductivity updated with state variables, thereby achieving dynamic updates of electrical conductivity.

6. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 5, characterized in that: The specific steps of the UMATHT subroutine in step S4 are as follows: First, write the UMATHT subroutine interface; second, input the current temperature, temperature gradient, and state variables; and finally, output the heat flux and associate it with the thermal conductivity in the material definition.

7. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: In step S4, the HETVAL subroutine applies an external stimulation current and calculates the equivalent heat flux generated by the ion current.

8. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 7, characterized in that: The specific steps of the HETVAL subroutine in step S4 are as follows: First, write the HETVAL subroutine interface; second, input the current temperature, time, and state variables; and finally, calculate the equivalent heat generation. The HETVAL subroutine directly converts the calculated temperature-dependent ion current into the force driving the evolution of the temperature field, forming a closed-loop feedback of temperature-gated-current-temperature.

9. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: The meshing in step S1 is specifically as follows: in the two-dimensional axisymmetric model, axisymmetric quadrilateral elements are used for structured meshing; in the three-dimensional finite element model, hexahedral dominant meshes are used, combined with wedge elements to handle the transition region, and for the cell membrane layer with a thickness of only 3 nm, the minimum mesh size is set to be less than 3 nm to fully analyze its physical properties.

10. The method for simulating the electrothermal coupling of neuronal action potentials based on temperature effects according to claim 1, characterized in that: It also includes step S5: verifying HH kinetics using a two-dimensional axisymmetric model; verifying electrothermal equivalence using a three-dimensional finite element model; the electrothermal equivalence verification specifically involves applying a voltage boundary to the outer surface of the extracellular medium, and in the axisymmetric verification, verifying the correctness of the electrothermal equivalence by comparing the membrane potential decay law along the axial direction with the analytical solution of the passive cable equation; the HH kinetic verification specifically involves verifying the correctness of embedding temperature-dependent ion channel kinetics in the custom subroutine by using the action potential waveform, gating variable, and ion current calculated by Abaqus after applying a stimulation current;

Images