Software implementation of the multivariate method for the Hodgkin-Huxley model

Abstract

tree induction. The paper developed and investigated estimates of computational complexity for multivariate methods of qualitative analysis of the Hodgkin-Huxley system for the purpose of classifying different types of nerve cell excitability using algorithms of sequential coverage and decision The developed method consists of 5 stages. The approach is proven as software in a package of multivariate method, qualitative analysis, Hodgkin-Huxley model, classification rules.

1. Introduction

Ukraine

2022 Copyright for this paper by its authors.

The general ideas of the method were developed in work [8] for the case of the initial conditions of the ODE. In this work, it will be developed also for the speed constants of the ODE. At the same time, we use the Monte Carlo approach, which consists in the random generation of parameters and the construction of the ODE model based on them. Next, we will apply the sequential coverage algorithm to build classification rules

The model of electrical activity of the giant squid axon proposed in [9]. In the model, each component of the excitable cell is considered as an electrical element. The lipid layer is represented as a container Cm . Ion channels are represented by electrical conductance gi , where i is a specific ion channel that depends on both voltage and time. Ion pumps are represented by a current source I app . Denote by V the difference between the membrane potential and the residual potential.

The current through the bilipid layer will be:

dV .

Ic  Cm dt

Ii  gi (V  Vi ) , The current through a given ion channel will be: where Vi ˗ is the equilibrium potential of the i -th ion channel.

For cells with potassium, sodium and chlorine channels, the total current through the membrane I will be:

I  I c  I K  I Na  I L .

The final typical Hodgkin-Huxley model is: ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) dV dt dm dt dn dt  g K n4 (V VK )  g Na m3h(V VNa )  g L (V VL )  I app ,  (1 m) *0.1*  (1 n) * 0.1* dh dt 25 V

25V exp 10 1 10 V 10 10V exp 10 1 

V 0.07*exp 20 *(1h)

V  m * 4*exp 18 ,

V  n * 0.125* exp 80 ,

1 x(t) limln

tx00t x0

A method of identifying model parameters based on the quadratic quality criterion is proposed. An algorithm for studying the chaotic nature of the attractor using the Lyapunov exponent method has been developed. At the same time, the Lyapunov exponential λ is the average speed of divergence of two trajectories, which is determined from the ratio

x(t)etx0 as

The numerical determination of the exponents of the Lyapunov system of differential equations in this paper is based on the methodology proposed in the works of J. Argyris, G. Faust, and M. Haase. At the same time, the Lyapunov exponents are determined by the transition along the main axis from the center of the infinitesimal sphere. The center of the sphere is obtained on the basis of nonlinear differential equations under certain initial conditions. The trajectories of points on the surface of the sphere are determined on the basis of linearized differential equations at points infinitely far from the center of the sphere. The main axis is determined by linearized equations and a set of orthonormal dm dt dn dt dh dt dt neighborhood of a certain stationary state (V * , m* , n* , h* ) .

The existence of complex (even chaotic) behavior in the Hodgkin-Huxley system indicates the need to consider control in the model - primarily due to the external applied current.

Therefore, the problem of optimal bifurcation control in the Hodgkin-Huxley system is considered.

We have a management system: dV

 g K n4 (V  VK )  g Na m3h(V  VNa )  g L (V  VL )  I appu , ( 6 ) Multiple controls: Here a,b,t f  0 .

U  {u(t) : a  u(t)  b, t0  t  t f , u(t)  calculable } . and bounded, the right-hand side of the control system is continuous x with respect to and only measurable with respect to t for a fixed x . Therefore, the solutions of the system are absolutely continuous functions that satisfy ( 10 ) almost everywhere. Under such conditions, the existence of a solution is proven.

The problem of optimal control contains the criterion of the quality J[u] of the view: t f J[u]   L(t, x, u)dt   (x(t f )) ,

t0 where L – is a given really significant function and  – is a continuously differentiable really significant function. The goal is to find a control such that

J[u* ]  inf J[u] . ( 10 )

uU

The maximum principle for the Hodgkin-Huxley system is formulated. Optimal control in problem ( 6-10 ) exists if the integrand expression in the quality criterion is a convex function. At the same time, the trajectory of the system belongs to space.

The necessary optimality conditions have been obtained. The Hamilton-Pontryagin function has the form: ( 7 ) ( 8 ) ( 9 ) 25V (exp 10  1)2

30V exp 10 ) 25  V

25V ,

( 11 ) ( 12 ) According to the Pontryagin maximum principle, we have

u* (t)  sgn 1 (t) .

Therefore, the optimal control in problem ( 6 – 10 ) can be constructed as a result of solving the following boundary value problem:

 g K n4 (V  VK )  g Na m3h(V  VNa )  g L (V  VL )  I app sgn 1 (t) , dV dt dm dt dn dt

V  m * 4 * exp 18 ,

 V  n * 0.125 * exp 80 , dh dt h V (0)  V0 , m(0)  m0 , n(0)  n0 , h(0)  h0 , d1   H dt V

 2V  1 ( g K n4  g Na m3h  g L )  d4   H dt h 30V exp 10 ) ,

Theorem 1. For a sufficiently small value, the solution t f of system ( 12 ) is unique.

The constructed control as a result of the integration of the system ( 12 ), which is an external applied electric current, can be practically implemented. By looking at different groups of cells, this stabilizing control of bifurcation in the Hodgkin-Huxley system may have important clinical applications for patients suffering from Alzheimer's disease, epilepsy or arrhythmia.

A system based on ( 1 ) - ( 4 ) is assumed to exist at initial values and velocity parameters at specified intervals.

Options p  P  {( g K , g Na , g L ,VK ,VNa ,VL , Cm , xm , xn , xh ) : g Kmin  g K  g Kmax , g Nmain  g Na  g Nmaax , g Lmin  g L  g Lmax ,VKmin  VK  VKmax ,VNmain  VNa  VNmaax , VLmin  VL  VLmax , Cmmin  Cm  Cmmax , xmmin  xm  xmmax , xnmin  xn  xnmax , xhmin  xh  xhmax }  R10 , and the initial conditions (V0 , m0 , n0 , h0 )  X 0  {(V0 , m0 , n0 , h0 ) : V0min  V0  V0max , m0min  m0  m0max , n0min  n0  n0max , h0min  h0  h0max }  R 4 .

The method consists in randomly generating initial values and values of speed parameters that would belong to a practically reasonable region. For each of the sets of such parameters, the system ( 1 ) – ( 4 ) is integrated to obtain the corresponding trajectories. The algorithm of data mining technology (decision tree induction, sequential coverage method, etc.) is then applied to the obtained results in order to obtain certain knowledge structures for decision-making.

So, the approach includes the following five steps.

1. Definition of classes of trajectories of the system. It should be noted that in practical applications they are mostly dealing with much more complex behaviors in order to characterize their concepts as “stable-unstable” and accordingly to resort to the analysis of eigenvalues or Lyapunov exponents of the dynamic system. Therefore, it is advisable to transfer the definition of qualitative forms of trajectories to the competence of expert physiologists. In this case, we will use classes related to the types of neuronal excitability [10]: type I, type II, type III. To denote a trajectory class C , a class attribute is introduced that takes one of 3 discrete values C 1,3 . Figures 1-3 show typical representations of 3 classes of trajectories – types of neuron excitability when the applied current increases: - type І - type ІІ - type ІІІ

2. Generation of a matrix of random initial values and velocity parameters. In order to explore the entire space of initial values and velocity parameters with respect to the generation of trajectory classes defined in the first step, a matrix of random initial values and velocity parameters is generated based on probability distributions at defined intervals. In this paper, we assume that the initial values and velocity parameters are evenly distributed at intervals. Each column corresponds to a set of values of one parameter – either the initial value or the speed parameter. Each row is a set of initial values and velocity parameters for one run of the model based on the ODE:

 V01 m01 n01 h01 g1K g1Na g1L VK1 VN1a VL1 Cm1 x1m x1n x1h  M   V02 m02 n02 h02 g K2 g N2a g L2 VK2 VN2a VL2 Cm2 x m2 xn2 xh2   R N14 V0N m0N n0N h0N g KN g NNa g LN VKN VNNa VLN CmN xmN xnN xhN 

3. Run the model and classify the input set. Each set of initial values and velocity parameters generated in the second step is used as input for the Hodgkin-Huxley model. The numerical integration of the equations is carried out using the Adams method [11]. Output trajectories are classified on the basis of the criteria proposed in the first step. Based on the classification results, sets of initial values and velocity parameters are assigned to the corresponding class attributes.

4. Construction of a matrix of relationships between initial values and between velocity parameters. The method assumes that for the shape of the trajectories of the system, the relationship between the initial values and between the velocity values is much more important than their absolute values. Therefore, a matrix is constructed that includes information in categorized coded form about the relationship between the initial values and between the velocity parameters generated in step 2:  m0  n0 m0  h0 n0  h0 g K  g Na g K  g L g Na  g L VK VNa  x(m01 , n01 ) x(m01 , h01 ) x(n01 , h01 ) p(g1K , g1Na ) p(g1K , g1L ) p(g1Na , g1L ) p(VK1 ,VN1a ) D   x(m02 , n02 ) x(m01 , h01 ) x(n01 , h01 ) p(g1K , g1Na ) p(g1K , g1L ) p(g1Na , g1L ) p(VK1 ,VN1a )

 x(m0k , n0k ) x(m01 , h01 ) x(n01 , h01 ) p(g1K , g1Na ) p(g1K , g1L ) p(g1Na , g1L ) p(VK1 ,VN1a ) VK VL p(VK1 ,VL1 ) p(VK1 ,VL1 ) p(VK1 ,VL1 )

VNa VL p(VN1a ,VL1 ) p(VN1a ,VL1 ) p(VN1a ,VL1 ) 0, if  Here x(u, v)  p(u, v)  1, if

Cm 1 p(Cm1 ,1) p(Cm1 ,1) p(Cm1 ,1) 2, if xm  xn p(x1m , x1n ) p(x1m , x1n ) p(x1m , x1n ) u  v u  v , u  v xm  xh p(x1m , x1h ) p(x1m , x1h ) p(x1m , x1h ) xn  xh p(x1n , x1h ) p(x1n , x1h ) p(x1n , x1h )

C 

 C1  C   R k14

2  Ck  C 1,3 – values of the class attribute associated with the corresponding trajectory forms.

Therefore, in this step, the numerical values of the initial values and velocity parameters are transformed into categorical values of the attributes of the training data sets. Since the probability of equality of random numbers is zero, the matrix D looks like a "binarization" of the relations between the initial values and between the velocity parameters. That is, the matrix D will include only the values 0 and 2.

5. Application of sequential coverage algorithm to correlation between initial values and between velocity parameters. The binary ratio matrix D constructed in step 4 is the training data set for the sequential coverage algorithm. A built-in set of classification rules will contain a check of the relation between the initial values and the velocity parameters in antecedents. As a consequence of the rules will be the trajectory classes of the model C 1,3 .

The proposed multivariate method of qualitative analysis of models of electrophysiological processes is an approach that allows solving the problems of classification of cell excitability, which cannot be solved by other traditional methods, such as stability theory or limit cycles.

In general, the method combines the Monte Carlo approach for forming training sets and data mining classification algorithms: the sequential coverage method with the generation of classification rules and the decision tree induction method.

The advantages of the classification rules that can be built at the 5th step of the algorithm are that they correspond to the natural reflection of knowledge in people's thinking and are more expressive. In addition, the sequential coverage algorithm is easier to implement and debug compared to recursive decision tree algorithms, and its computational complexity is simpler compared to finite automata.

The decision tree induction method, being more difficult to implement, allows visual visualization and a priori values of probabilities for the type of cell excitability based on the relationships between initial values and speed parameters in the Hodgkin-Huxley model.

Software implementation of the multivariate method for the Hodgkin-Huxley model A software environment for the study of electrophysiological processes based on the Hodgkin-Huxley system has been developed in the form of a library of Java classes. In the future, such a software environment should be focused on the construction and research of electrophysiological models in all aspects of physical impact on biological tissue. Approbation of the built mathematical model was carried out by comparing numerical calculations and experimentally obtained data.

To implement the method, a package of Java classes rule.model has been developed. The package includes classes [12] (Fig. 4):

DataManager – class – data manager for obtaining information from the database through the mediation of the corresponding servant classes;

MultiVariateMethod – a class for implementing the multivariate method presented in the work – the main class of the package;

TuplesPeer is a servant class for forming and processing training sets that will be used in the classification algorithm. rate_constants_values – matrix of randomly generated rate constants.

To implement the method of sequential coverage with the construction of classification rules, a package of Java-classes rule.model has been developed. The package includes classes: beans-classes Attribute, Attribute_for_list for working with the data of the corresponding tables, and Rule – for presenting rules. SQL – queries for obtaining relevant data are implemented in the AttributeListPeer and TuplesPeer classes.

The Rule_set class stores a set of training rules. In addition, this class directly implements the sequential coverage algorithm. The class contains members: data manager m_dataManager, hash tables of training data sets m_htTuples, all attributes with their possible values m_htAtt_vals and directly a set of rules m_htRule_set.

In the constructor of the Rule_set class, the m_htTuples and m_htAtt_vals hash tables are built, as well as the sequential coverage algorithm is applied by calling the Sequential_covering (m_htTuples, m_htAtt_vals) method. The resulting set of rules is output to a text file.

The Rule class is designed to store individual rules. Its class members are two hash tables: m_htAntecedent – for storing the antecedent of the rule and m_htConsequent – for the consequent. Using the method

public void conjunctCondition(Attribute_for_list attribute, String sAttribute_value) the conjunction of the new condition to the rule is carried out. Using the method

public Rule copy() a "deep" copy of the rule is created. In this case, the JOS (Java Object Serialization) protocol is used.

Counting the number of positive and negative training sets is carried out in the methods of the TuplesPeer class.

The built decision tree for n  17 is shown in Figure 7. The time to build the decision tree is 1402 ms.

Conclusions

The results show that in the future, such a software environment should be focused on the construction and research of electrophysiological models in all aspects of physical impact on biological tissue and further research can use the obtained results to study cyber physical systems, as an example the obtained results can be used to study the cyber physical systems [12, 13]. Approbation of the built mathematical model was carried out by comparing numerical calculations and experimentally obtained data.

References

[1]. B. Hassard , Bifurcation of periodic solutions of the Hodgkin-Huxley model for the squid giant axon , J. Theoret. Biol. (71) , pp. 401 - 420 , 1978 .

[2] . H. Fukai , T. Nomura , S. Doi , and S. Sato , Hopf bifurcations in multiple-parameter space of the Hodgkin-Huxley equations, I, II, Biol . Cybern., ( 82 ), pp. 215 - 222 ; 223 - 229 , 2000 .

[3] . J. Guckenheimer and I. S. Labouriau , Bifurcation of the Hodgkin-Huxley equations: A new twist , Bull. Math. Biol., (55) , pp. 937 - 952 , 1993 .

[4] . B. Hassard and L. -J. Shiau , Asp ecial point of Z2-codimension three Hopf bifurcation in the Hodgkin-Huxley model , Appl. Math. Lett., (9) , pp. 31 - 34 , ( 1996 ).

[5]. S. Labouriau, Degenerate Hopf bifurcation and nerve impulse II, SIAM J . Math. Anal., ( 20 ), pp. 1 - 12 , 1989 .

[6]. Maihruk Z.V. Software implementation of the numerical method of optimal bifurcation control in the Hodgkin-Huxley model . / Maihruk Z. // Bulletin of the Khmelnytskyi National University No. 1, pp. 186 - 194 , 2014 .

[7]. Martsenyuk V. P. , Mayhruk Z. V. : Construction of optimal bifurcation control in the HodgkinHuxley model based on the maximum principle. Scientific journal of the V.M. Glushkov Institute of Cybernetics of the National Academy of Sciences of Ukraine Kamianets-Podilsky National University named after Ivan Ogienko Mathematical and Computer Modeling Series: Technical Sciences (9) , 78 - 90 , 2014 .

[8]. Yvonne Koch, Thomas Wolf, Peter K. Sorger , Roland Eils, Benedikt Brars, Decision-Tree Based Model Analysis for Efficient Identification of Parameter Relations Leading to Different Signaling States, www .plosone.org, December 2013 , Volume 8 , Issue 12 , e82593 .

[9] . Hodgkin , A. L. ; Huxley , A. F. : A quantitative description of membrane current and its application to conduction and excitation in nerve . The Journal of physiology 117 (4): PMC 1392413. PMID 12991237 , рр. 500 - 544 , ( 1952 ).

[10]. J. R. Clay et al. A simple modification of the Hodgkin and Huxley equations explains type 3 excitability in squid giant axons , J. R. Soc. Interface 5 , 1421 - 1428 , ( 2008 ).

[11]. Hairer E., Nersett S. , Vanner G.: Solution of ordinary differential equations . Tough tasks . - M.: Mir, p. 512 , 1990 .

[12]. V. Martsenyuk , Z. Mayhruk , M. Karpinski , N. Milian , I. Andrushchak , O. Veselska , On implementation of decision tree induction in cloud platforms , in: Advances in Science and Engineering Technology International Conferences (ASET) , Dubai, United Arab Emirates, 2019 , P. 1- 6 . doi: 10 .1109/ICASET. 2019 . 8714450 .

[13]. V. Martsenyuk , A. Sverstiuk , O. Bahrii-Zaiats , Yu. Rudyak, B. Shelestovsky , Software complex in the study of the mathematical model of cyber-physical systems , ICT&ES 2nd International Workshop Information-Communication Technologies & Embedded Systems , 12 November, 2020 Mykolaiv, Ukraine. p. 87 - 97 .

[14]. V. Martsenyuk , A. Klos-Witkowska , A. Sverstiuk , O. Bahrii-Zaiats , M. Bernas , K. Witos , Intelligent Big Data system based on scientific machine learning of cyber-physical systems of medical and biological processes , in: The Fourth International Workshop on Computer Modeling and Intelligent Systems (CMIS-2021). 27 of April 2021 , Zaporizhzhia, Ukraine. pp. 34 - 48 . doi: 10 .32782/cmis/2864-4.