=Paper=
{{Paper
|id=Vol-3309/paper9
|storemode=property
|title=Software implementation of the multivariate method for the Hodgkin-Huxley model
|pdfUrl=https://ceur-ws.org/Vol-3309/paper9.pdf
|volume=Vol-3309
|authors=Vasyl Martsenyuk,Zoryana Mayhruk,Oleksandra Kuchvara,Oksana Bahrii-Zaiats,Igor Andrushchak
|dblpUrl=https://dblp.org/rec/conf/ittap/MartsenyukMKBA22
}}
==Software implementation of the multivariate method for the Hodgkin-Huxley model==
https://ceur-ws.org/Vol-3309/paper9.pdf
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* 25V
m * 4 * exp 18 , (2)
dt
exp 10 1
10 V
V
dn
(1 n) * 0.1* 10 n * 0.125 * exp 80 , (3)
10V
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
tlimln
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 * 25V
m * 4 * exp ,18
(7)
dt
exp 10
1
10 V
V
dn
(1 n) * 0.1 * 10 n * 0.125 * exp , 80
(8)
10V
dt
exp 10
1
V
dh h
0.07 * exp 20 * (1 h) 30V
. (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)
uU
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 * 25V
m * 4 * exp 18 )
exp 1 10
10 V
V
3 ((1 n) * 0.1 * 10 n * 0.125 * exp )
80
10V
exp 10 1
V
h
4 (0.07 * exp 20
* (1 h) 30V
)
1 exp 10
We have a coupled system:
d1 H
2V 1 ( g K n 4 g Na m 3 h g L )
dt V
25V 25V
25 V
1 exp 10
exp 10
V
10 2
2 ((1 m) * 0.1 * 25V
m * * exp 18
)
9
(exp 10 1) 2
10V 10V ,
10 V
1 exp 10
exp 10
V
10 0.125
3 ((1 n) * 0.1 * 10V
n* * exp 80 )
80
(exp 10 1) 2
V 30V
0.07 h
4 ( * exp 20 * (1 h) 30V
exp 10
) (11)
20
10 * (1 exp 10
) 2
d2 H 25 V
V
1 * 3 * g Na m 2 h(V VNa ) 2 (0.1 * 25V
4 * exp 18 ) ,
dt m
exp 10 1
10 V
d3 H
V
1 g K * 4 * n 3 (V VK ) 3 (0.1 * 10 0.125 * exp 80
),
10V
dt n
exp 10 1
d4 H
V
1
1 g Na m 3 (V VNa ) 4 (0.07 * exp 20 30V
).
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 * 25V
m * 4 * exp 18
,
dt (12)
exp 10
1
10 V
V
dn
(1 n) * 0.1 * 10 n * 0.125 * exp , 80
10V
dt
exp 10
1
V
dh h
0.07 * exp 20 * (1 h) 30V
.
dt
1 exp 10
V (0) V0 , m(0) m0 , n(0) n0 , h(0) h0 ,
d1 H
2V 1 ( g K n 4 g Na m 3 h g L )
dt V
25V 25V
25 V
1 exp 10
exp 10
V
10 2
2 ((1 m) * 0.1 * 25V
m * * exp 18
)
9
(exp 10 1) 2
10V 10V ,
10 V
1 exp 10
exp 10
V
10 0.125
3 ((1 n) * 0.1 * 10V
n* * exp 80 )
80
(exp 10 1) 2
V 30V
0.07 h
4 ( * exp 20 * (1 h) 30V
exp 10
)
20
10 * (1 exp 10
) 2
d2 H 25 V
V
1 * 3 * g Na m 2 h(V VNa ) 2 (0.1 * 25V
4 * exp 18 ) ,
dt m
exp 10 1
10 V
d3 H
V
1 g K * 4 * n 3 (V VK ) 3 (0.1 * 10 0.125 * exp 80
),
10V
dt n
exp 10 1
d4 H
V
1
1 g Na m 3 (V VNa ) 4 (0.07 * exp 20 30V
),
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