Software implementation of the multivariate method for the Hodgkin-Huxley model Vasyl Martsenyuka, Zoryana Mayhrukb, Oleksandra Kuchvarab, Oksana Bahrii-Zaiatsb, Igor Andrushchakc a University of Bielsko-Biala, Willowa St. 2, Bielsko-Biala, 43-300, Poland b I. Horbachevsky Ternopil National Medical University, 12 Rus'ka St., Ternopil, 46001, Ukraine c Lutsk National Technical University, Lutsk, Lvivska St. 75, Ukraine Abstract 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 tree induction. The developed method consists of 5 stages. The approach is proven as software in a package of Java classes. Keywords multivariate method, qualitative analysis, Hodgkin-Huxley model, classification rules. 1. Introduction Despite significant achievements in the field of mathematical modeling of medical and biological processes, a number of problems arise when using the approach to studying the electrical activity of cells. The vast majority of studies concern the construction of models of ionic conductivity of cell membranes based on the Hodgkin-Huxley formalism [1-7]. The main practical results concern the establishment of an approximate analytical form of velocity parameters in nonlinear differential equations by statistical methods for different groups of cells. At the same time, obtaining results regarding the qualitative characteristics of cells (for example, classification of types of excitability) in an analytical form (in the form of conditions that include numerical parameters) based on a system of the Hodgkin-Huxley type does not seem possible due to complex nonlinearities in the equations. One possible approach used is to reduce to a two-dimensional Fitzhugh-Nagumo model, which is based on certain restrictions on the heuristic model. Therefore, the work offers a fundamentally different approach to the qualitative analysis of the Hodgkin-Huxley system, which uses a multivariate method of simulation modeling with the subsequent application of data mining algorithms. The criteria obtained in this way can be presented in the form of knowledge structures - classification rules, or decision trees [8]. The results make it possible to obtain conditions for the classification of types of excitability for certain groups of cells, and implementation in the control model can be applied to specific pathologies and certain groups of patients. Thus, the development of a method of qualitative analysis of the electrical activity of the cell for the Hodgkin-Huxley system, its mathematical and computer implementation, is relevant. The purpose of this work is to improve the methods and means of mathematical, computer modeling and computing methods for the integration of these methods and their results into information technologies for the diagnosis of electrophysiological processes of cell excitability. Research methods. Qualitative analysis of the Hodgkin-Huxley model using the multivariate method. The task of the method is to establish the mechanisms of multiparametric effects in the Hodgkin-Huxley model. ITTAP’2022: 2nd International Workshop on Information Technologies: Theoretical and Applied Problems, November 22–24, 2022, Ternopil, Ukraine EMAIL: vmartsenyuk@ath.bielsko.pl, majhruk@tdmu.edu.ua, kuchvara@tdmu.edu.ua, bagrijzayats@tdmu.edu.ua, 9000@lntu.edu.ua ORCID: 0000-0001-5622-1038 (A.1); 0000-0001-8909-0644(A.2); 0000-0002-0248-3224 (A.3); 0000-0002-5533-3561 (A.4); 0000-0002- 8751-4420 (A.5) ©️ 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) 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 C m . Ion channels are represented by electrical conductance g i , 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 . I c  Cm dt The current through a given ion channel will be: I i  g i (V  Vi ) , 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: dV   g K n 4 (V  VK )  g Na m 3 h(V  VNa )  g L (V  VL )  I app , (1) dt 25  V V dm   (1  m) * 0.1* 25V  m * 4 * exp 18 , (2) dt exp 10  1 10  V V dn   (1  n) * 0.1* 10  n * 0.125 * exp 80 , (3) 10V dt exp 10  1 V dh  h  0.07*exp 20 *(1 h)  dt 30 V 1 exp 10 . (4) 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  1  x () t tlimln  t x  0 0 x0 (5) 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 vectors attached to the center of the sphere. The Gram-Schmidt method is used to construct an orthonormal basis. The linearization of the system of nonlinear differential equations (1) – (4) was carried out in the 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 n 4 (V  VK )  g Na m 3 h(V  VNa )  g L (V  VL )  I appu , (6) dt 25  V V dm   (1  m) * 0.1 * 25V  m * 4 * exp ,18 (7) dt exp 10 1 10  V V dn   (1  n) * 0.1 * 10  n * 0.125 * exp , 80 (8) 10V dt exp 10 1 V dh  h  0.07 * exp 20 * (1  h)  30V . (9) dt 1  exp 10 Multiple controls: U  {u (t ) : a  u (t )  b, t 0  t  t f , u (t )  calculable } . Here a, b, t f  0 . 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: tf 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: H  V 2  1 ( g K n 4 (V  VK )  g Na m 3 h(V  V Na )  g L (V  VL )  I appu )  25  V V    2 ((1  m) * 0.1 * 25V  m * 4 * exp 18 )  exp 1 10 10  V V   3 ((1  n) * 0.1 * 10  n * 0.125 * exp )  80 10V exp 10  1 V  h   4 (0.07 * exp 20 * (1  h)  30V ) 1  exp 10 We have a coupled system: d1 H   2V  1 ( g K n 4  g Na m 3 h  g L )  dt V 25V 25V 25  V 1  exp 10  exp 10  V 10 2  2 ((1  m) * 0.1 * 25V  m * * exp 18 ) 9 (exp 10  1) 2 10V 10V , 10  V 1  exp 10 exp 10  V 10 0.125  3 ((1  n) * 0.1 * 10V  n* * exp 80 )  80 (exp 10  1) 2 V 30V 0.07  h  4 ( * exp 20 * (1  h)  30V exp 10 ) (11) 20 10 * (1  exp 10 ) 2 d2 H 25  V V    1 * 3 * g Na m 2 h(V  VNa )  2 (0.1 * 25V  4 * exp 18 ) , dt m exp 10  1 10  V d3 H V    1 g K * 4 * n 3 (V  VK )  3 (0.1 * 10  0.125 * exp 80 ), 10V dt n exp 10  1 d4 H V  1   1 g Na m 3 (V  VNa )  4 (0.07 * exp 20  30V ). dt h 1  exp 10 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: dV   g K n 4 (V  VK )  g Na m 3 h(V  VNa )  g L (V  VL )  I app sgn 1 (t ) , dt 25  V V dm   (1  m) * 0.1 * 25V  m * 4 * exp 18 , dt (12) exp 10 1 10  V V dn   (1  n) * 0.1 * 10  n * 0.125 * exp , 80 10V dt exp 10 1 V dh  h  0.07 * exp 20 * (1  h)  30V . dt 1  exp 10 V (0)  V0 , m(0)  m0 , n(0)  n0 , h(0)  h0 , d1 H   2V  1 ( g K n 4  g Na m 3 h  g L )  dt V 25V 25V 25  V 1  exp 10  exp 10  V 10 2  2 ((1  m) * 0.1 * 25V  m * * exp 18 ) 9 (exp 10  1) 2 10V 10V , 10  V 1  exp 10 exp 10  V 10 0.125  3 ((1  n) * 0.1 * 10V  n* * exp 80 )  80 (exp 10  1) 2 V 30V 0.07  h  4 ( * exp 20 * (1  h)  30V exp 10 ) 20 10 * (1  exp 10 ) 2 d2 H 25  V V    1 * 3 * g Na m 2 h(V  VNa )  2 (0.1 * 25V  4 * exp 18 ) , dt m exp 10  1 10  V d3 H V    1 g K * 4 * n 3 (V  VK )  3 (0.1 * 10  0.125 * exp 80 ), 10V dt n exp 10  1 d4 H V  1   1 g Na m 3 (V  VNa )  4 (0.07 * exp 20  30V ), dt h 1  exp 10 k (t f )  0, k  1,4 . 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 , min g Na  g Na  g Na max , g Lmin  g L  g Lmax , VKmin  VK  VKmax , VNamin  VNa  VNamax , 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 І Figure 1. Neuronal excitability and type І. - type ІІ Figure 2. Neuronal excitability and type ІІ. - type ІІІ Figure 3. Neuronal excitability and 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 g 1K g 1Na g 1L VK1 1 V Na VL1 C m1 x1m x1n x1h    M   V02 m02 n02 h02 g K2 2 g Na g L2 VK2 V Na2 VL2 C m2 x m2 x n2 x h2   R N 14 V N m0N n0N h0N g KN N g Na g LN VKN V NaN VLN C mN x mN x nN x hN   0 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  n 0 m0  h0 n0  h0 g K  g Na gK  gL g Na  g L VK  V Na  1 1 1 1 1 1  x(m0 , n0 ) x(m0 , h0 ) x(n0 , h0 ) p( g 1K , g 1Na ) p( g 1K , g 1L ) p( g 1Na , g 1L ) p(VK1 ,VNa 1 ) D x(m0 , n0 ) x(m0 , h0 ) x(n0 , h01 ) 2 2 1 1 1 p( g 1K , g 1Na ) p( g 1K , g 1L ) p( g 1Na , g 1L ) 1 1 p(VK ,VNa )   x ( m k , n k ) x ( m1 , h 1 ) x ( n 1 , h 1 ) p( g 1K , g 1Na ) p( g 1K , g 1L ) p( g 1Na , g 1L ) p(VK1 ,VNa 1  0 0 0 0 0 0 ) VK  VL VNa  VL Cm  1 xm  xn xm  xh xn  xh C  p(VK1 ,VL1 ) 1 p(VNa ,VL1 ) p(C m1 ,1) p( x1m , x1n ) p( x1m , x1h ) 1 1 p( xn , xh ) C1   R k 14 p(VK1 ,VL1 ) 1 p(VNa ,VL1 ) p(C m1 ,1) p( x1m , x1n ) p( x1m , x1h ) p( x1n , x1h ) C 2   p(VK1 ,VL1 ) 1 p(VNa ,VL1 ) p(C m1 ,1) p( x1m , x1n ) p( x1m , x1h ) p( x1n , x1h ) C k  0, if u  v Here x(u, v)  p(u, v)   1, if u  v , 2, if u  v  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. Figure 4. Package decision_tree.fde.hh Figure 5. UML-class diagram MultiVariateMethod. In the MultiVariateMethod class (Fig. 5.), random parameter values are generated (step 2): M_x0 = dm.getRandomInitialValues(); M_rateConstants = dm.getRandomRateConstants(); Next, the Hodgkin-Huxley system integration class-applet is launched. At the same time, the expert chooses the shape of the obtained trajectory (step 3). After that, the parameter interrelation matrix generation step is launched (step 4). Note that the sequence of steps 2-4 can be performed as many times as desired. At any moment, the user can run the decision tree induction algorithm (step 5): dtDecision_tree_HH = new decision_tree.fde.hh.DecisionTree(dmtnRoot, dataManager_FDE, htAttribute_list, "Information gain"); The hh database used in the package is implemented in the MySQL DBMS. It includes the following tables (Fig. 6): attribute – description of attributes for building a decision tree, i.e. relationships between initial values and between rate constants; categorized_data – training sets used in the classification algorithm (in this case, decision tree induction) and represent the matrix in the fourth step; init_values_values – matrix of randomly generated initial values; initial_values – description of initial values (including minimum and maximum values); rate_constants – description of rate constants (including minimum and maximum values); 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. Figure 7. Decision tree A set of classification rules is built for n  397 the following: IF x_mx_h THEN class=type I IF n_0v_L AND x_m>x_h AND x_n