Qualitative Analysis of the Hodgkin-Huxley Model of Neuron Excitability based on Classification Rules Vasyl Martsenyuk1[0000-0001-5622-1038], Yurii Rudyak 2[0000-0003-1836-9132], Andriy Sverstiuk2[0000-0001-8644-0776], Zorana Mayhruk2[0000-0001-8909-0644], Andriy Horkunenko2[0000-0002-2021-006X], Mykhailo Kasianchuk3[0000-0002-4469-8055] 1 University of Bielsko-Biala, Department of Іnformatics and Automatics, Poland 2 I. Gorbachevsky Ternopil National Medical University, Ternopil, Ukraine 3 Ternopil National Economic University, Ternopil, Ukraine vmartsenyuk@ath.bielsko.pl, sverstyuk@tdmu.edu.ua, majhruk@tdmu.edu.ua, horkunenkoab@tdmu.edu.ua Abstract. The proposed multivariate method of qualitative analysis of models of electrophysiological processes is an approach that allows to solve problems of classification of excitability of cells which cannot be solved by other tradi- tional methods, for example stability theory, or boundary cycles. As a whole, the method combines the Monte Carlo approach for the formation of training kits and the data mining classification algorithms: the sequential coverage method with the generation of classification rules and the induction method of the decision tree. The advantages of the classification rules, which can be built on the 5th step of the algorithm, are that they correspond to the natural reflec- tion of knowledge in the thinking of people and are more expressive. In addi- tion, the sequential coverage algorithm is easier to implement and debug than the recursive algorithms of decision trees, and its computational complexity is simpler than that of finite state machines. The decision tree induction method, being more complex to implement, allows visualization and a priori probability values for the type of cell excitability based on the relationship between the ini- tial values and the velocity parameters in the Hodgkin-Huxley model. The de- veloped method consists of 5 stages. The approach is proven in the form of software in the Java-class package. Keywords: Qualitative analysis Hodgkin-Huxley model, classification rules. 1 Introduction This paper examines the electrical activity model of a giant squid axon based on the Hodgkin-Huxley mathematical model [1-7]. When studying the excitability of the axon through the construction of classification rules of the type of excitability, the initial conditions are calculated. Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). The systems of ordinary differential equations of mathematical biology use veloci- ty constants and initial values as parameters. In this paper, we present an automated multivariate method that determines the conditions for parameter combinations to obtain a particular type of neuron excitability in the Hodgkin-Huxley model. At work [8] proposed a multivariate approach based on a decision tree that investi- gates only the influence of initial conditions in ODE (systems of ordinary differential equations) at fixed velocity constants. In real models, it is rarely possible to use or directly measure the values of velocity constants in the literature. As a rule, they are estimated based on time series [9] based on solution trajectories [10]. 2 Development of a multivariate method of qualitative analysis of the Hodgkin-Huxley model 2.1 The objective of the method is to establish mechanisms of multiparameter effects in the Hodgkin-Huxley model. General ideas of the method were developed in [8] for the case of the initial condi- tions of the ODE. In this paper, it will also be developed for the velocity constants of ODE. In doing so, 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 apply the sequential coverage algorithm to construct the classification rules. The model of electrical activity of a giant squid axon, proposed in [11], is consid- ered. In the model, each component of the excitatory cell is considered as an electrical element. The lipid layer is represented as a container Cm . Ion channels are represented by electrical conductivity g i , where i is a specific ion channel that depends on both voltage and time. Ion pumps are a current source I . Denote by V the difference app 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 m3h(V  VNa )  g L (V  VL )  I app dt , (1) dm 25  V  V  (1  m) * 0.1* 25V  m * 4 * exp 18 dt exp 10  1 , (2) 10  V V dn   (1  n) * 0.1 * 1010V  n * 0.125 * exp 80 dt exp 10  1 , (3) V dh  h  0.07 * exp * (1  h)  20 30 V dt 1  exp 10 . (4) A system based on (1) - (4) is assumed to exist at initial values and velocity param- eters 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 Namin  g Na  g Namax , 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 . We will then randomly generate the initial and velocity parameters that belong to a practically valid region. For each of the sets of such parameters, the system (1) - (4) is integrated to obtain the corresponding trajectories. The results are further applied by the sequential coverage algorithm to find specific patterns for decision making. So, as a whole, the following five steps are taken. 1. Definition of classes of trajectories of the system. It should be noted that in practi- cal applications they are mostly dealing with much more complex behaviors in or- der to characterize their concepts as “stable-unstable” and accordingly to resort to the analysis of eigenvalues or Lyapunov exponents of the dynamic system. There- fore, 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 [12]: type I, type II, type III. To denote a trajecto- ry 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 І Fig. 1. Neuronal excitability and type І. - type ІІ Fig. 2. Neuronal excitability and type ІІ. - type ІІІ Fig. 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 ini- tial values and velocity parameters is generated based on probability distributions at defined intervals. In this paper, we assume that the initial values and velocity pa- rameters 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 m10 n10 h01 g 1K g 1Na g 1L VK1 1 VNa VL1 Cm1 x1m x1n x1h   M   V02 m02 n02 h02 g K2 2 g Na g L2 VK2 VNa 2 VL2 Cm2 xm2 xn2 xh2   R N 14  N  V0 m0N n0N h0N g KN N g Na g LN VKN N VNa 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 [13]. 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 veloci- ty parameters are assigned to the corresponding class attributes. 4. Construction of a matrix of relationships between initial values and between ve- locity 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 con- structed 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  VNa   x(m0 , n0 ) x(m0 , h0 ) x(n0 , h0 ) p( g K , g Na ) p( g K , g L ) p( g Na , g L ) p(VK ,VNa ) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 D x(m0 , n0 ) x(m0 , h0 ) x(n0 , h0 ) p( g K , g Na ) p( g K , g L ) p( g Na , g L ) p(VK ,VNa1 ) 2 2 1 1 1 1 1 1 1 1 1 1 1  k k 1 1 1 1 1 1 1 1 1 1 1 1  x(m0 , n0 ) x(m0 , h0 ) x(n0 , h0 ) p( g K , g Na ) p( g K , g L ) p( g Na , g L ) p(VK ,VNa ) VK  VL VNa  VL C m  1 xm  xn xm  xh xn  xh C  p (VK1 ,VL1 ) p (VNa1 ,VL1 ) p (C m1 ,1) p ( x1m , x1n ) p ( x1m , x1h ) p ( x1n , x1h ) C1   R k 14 p (VK1 ,VL1 ) p (VNa1 ,VL1 ) p (C m1 ,1) p ( x1m , x1n ) p ( x1m , x1h ) p ( x1n , x1h ) C 2   p (VK1 ,VL1 ) p (VNa1 ,VL1 ) p (C m1 ,1) p ( x1m , x1n ) p ( x1m , x1h ) p ( x1n , x1h ) C k  Here 0, if u  v , C  1,3 – values of the class attribute associated  l x(u, v)  p(u, v)  1, if uv 2, if uv  with the corresponding trajectory forms. Therefore, in this step, the numerical values of the initial values and velocity pa- rameters are transformed into categorical values of the attributes of the training da- ta 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 val- ues 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 . 3 Software implementation of the multivariate method for the Hodgkin-Huxley model A rule.model Java-package was developed to implement the method rule.model. The package includes classes (fig.4): beans-classes Attribute, Attribute_for_list to work with the data in the corresponding tables and Rule to represent the rules. SQL queries to obtain the relevant data are implemented in the AttributeListPeer and TuplesPeer classes [14]. The Rule_set class stores a set of training rules. In addition, this class implements the sequential coverage algorithm directly. The class contains members: the m_dataManager data manager, the hash tables of the m_htTuples training datasets, all the attributes with their possible m_htAtt_vals values, and directly a set of m_htRule_set rules. The Rule_set class constructor constructs the hash tables m_htTuples and m_htAtt_vals, as well as applying the sequential coverage algorithm by calling the Sequential_covering method (m_htTuples, m_htAtt_vals). 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 storing the rule. Using the method public void conjunctCondition(Attribute_for_list attribute, String sAttribute_value) the new condition is conjugated to the rule. Using the method public Rule copy() the new condition is conjugated to the rule. Using the method a "deep" copy of the rule is created. This uses the JOS (Java Object Serialization) protocol. The number of positive and negative training sets is calculated using the methods of the TuplesPeer class. Fig. 4. Package karule.model. In the fde.hh.MultiVariateMethod class (Fig. 5), random parameter values are gener- ated (step 2):  M_x0 = dm.getRandomInitialValues();  M_rateConstants = dm.getRandomRateConstants(); Next, the Hodgkin-Huxley system integration applet is launched. The expert chooses the shape of the obtained trajectory (step 3). After that, the step of generating the pa- rameter matrix is started (step 4). Note that the sequence of steps 2-4 can be per- formed as many times as desired. At any time, the user can run the sequential cover- age algorithm (step 5): rule.model.Rule_set rule_set = new rule.model.Rule_set(rule_dataManager, sql0, rules_file_url) Fig. 5. UML-class diagram MultiVariateMethod. The fde database used in the package is implemented in MySQL DBMS. It includes the following tables (Fig. 6): attribute – a description of the attributes used to construct the classification rules, that is, the relationship between the initial values and the velocity constants; categorized_data – training kits used in the classification algorithm (in this case, the sequential coverage algorithm) and represent the matrix in the fourth step; init_values_values – matrix of randomly generated initial values:  V01 m01 n01 h01   2   V0 m02 n02 h02   R N 4 V N mN n0N h0N   0 0 initial_values – description of initial values (including minimum and maximum values); parameter_kind – type of parameter; rate_constants – description of velocity constants (including minimum and maxi- mum values); rate_constants_values – matrix of randomly generated velocity constants:  g 1K g 1Na g 1L VK1 VNa 1 VL1 Cm1 x1m x1n x1h    g K2 g Na 2 g L2 VK2 VNa 2 VL2 Cm2 xm2 xn2 xh2   R N 10    g KN g Na N g LN VKN VNaN VLN CmN xmN xnN xhN  3.1 Features of software implementation of the Hodgkin-Huxley model The rule.model package can be used for a wide range of functional-differential equa- tion systems. For this purpose, the model based on functional-differential equations should be implemented in the form of a corresponding package of Java classes. An example of such a package in the case of the Hodgkin-Huxley model is the medbioin- vestigations.hodgkin_hyxley package (Fig. 7). To integrate with the decision_tree.fde.hh package in the Hodgkin_HuxleyGraph class, which renders the model graphically, a new constructor was added to the exist- ing one, using the instance reference to the MultiVariateMethod instance. This con- structor additionally creates a JComboBox class object that lets you choose the shape of the trajectory and run the method step 4: String[] classStrings = {"subclini- cal","chronic","acute","lethal"}; JComboBox m_jcbClassName = new JComboBox(classStrings); m_jcbClassName.addActionListener(new ActionListener() { public void actionPerformed(ActionEvent e) { JComboBox jcbClass = (JComboBox) e.getSource(); m_sClassName=(String)jcbClass.getSelectedItem(); mvm.m_sClassName = m_sClassName; ((AdvancedFrame)getParent()).dispose(); mvm.run4thStep(); } }); You should also make appropriate changes to the hh database table in the following order: ─ describe all initial values and velocity constants in the tables initial_values and rate_constants respectively; ─ describe the dependencies between the initial values and the velocity constants to be examined in the attribute table; ─ in the categorised_data table, create fields according to the attribute table data. Fig. 6. Database Tables hh. 3.2 Estimation of complexity of algorithm execution From the analysis of the sequential coverage algorithm, we see that the computational complexity is determined by the product of the number of possible values of the class attribute K  3 (the number of iterations of the outer loop) and the computational complexity of the Get_one_ rule (D, Att_vals, c) procedure that is performed inside each cycle. Fig. 7. Packet medbioinvestigations.hodgkin_huxley. The Get_One_Rule (D, Att_vals, c) procedure involves performing p  18 iterations. At each iteration for a particular attribute Ai , a measure is calculated FOIL _ Gain for each of the K i  2 attribute values. That is, the inner body of the loop in the p Get_one_ rule (D, Att_vals, c) procedure  K  36 is executed once. The measure i 1 i FOIL _ Gain is calculated as a result of 3 – SQL-queries that can be evaluated with complexity O(log( N )) (see MySQL 5.0 documentation - http://dev.mysql.com/doc/refman/5.0/en/select-speed.html). Therefore, in general, the Get_one_Regulation (D, Att_vals, c) procedure has computational complexi- ty O K  log( N )   O(36 * log( N )) . p  i   i 1  To sum up, we have the computational complexity of the entire sequential order algo- rithm O K   K i  log( N )   O(108 * log( N ))  O(log( N )) p  i 1  (5) A numerical experiment was conducted on the basis of model (1) - (4) to find out the consistency of the time of construction of a set of classification rules with an esti- mate (6). A system based on the Celeron (R) Dual-Core CPU T3300 @ 2.00 GHz and 2 GB RAM was used. The results are presented in Fig. 8. Fig. 8. Difficulty of sequential covering algorithm. Figure 8 shows a significant deviation from the estimation time of the construction of classification rules N  350 , which is associated with an increase in computing re- sources. A set of classification rules for n  397 below. IF x_mx_h THEN class=type I IF n_0v_L AND x_m>x_h AND x_n