BioMed Research International

Volume 2015, Article ID 936295, 13 pages

http://dx.doi.org/10.1155/2015/936295

## Application of Stochastic Automata Networks for Creation of Continuous Time Markov Chain Models of Voltage Gating of Gap Junction Channels

^{1}Department of Mathematical Modelling, Kaunas University of Technology, Studentų Street 50, 51368 Kaunas, Lithuania^{2}Laboratory of Systems Control and Automation, Lithuanian Energy Institute, Breslaujos Street 3, 44403 Kaunas, Lithuania^{3}Department of Applied Informatics, Vytautas Magnus University, Vileikos Street 8-409, 44404 Kaunas, Lithuania^{4}Department of Business Informatics Research in Systems, Kaunas University of Technology, Studentų Street 56, 5142 Kaunas, Lithuania^{5}Department of Anesthesiology, Albert Einstein College of Medicine, 1300 Morris Park Avenue, Bronx, NY 10461, USA^{6}Department of Anesthesiology, New York Hospital Queens, 56-45 Main Street, Flushing, NY 11355, USA^{7}Institute of Cardiology, Lithuanian University of Health Sciences, Sukileliu Street 17, 50009 Kaunas, Lithuania^{8}Department of Neuroscience, Albert Einstein College of Medicine, 1300 Morris Park Avenue, Bronx, NY 10461, USA

Received 4 July 2014; Revised 7 December 2014; Accepted 8 December 2014

Academic Editor: Carlo Cattani

Copyright © 2015 Mindaugas Snipas et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The primary goal of this work was to study advantages of numerical methods used for the creation of continuous time Markov chain models (CTMC) of voltage gating of gap junction (GJ) channels composed of connexin protein. This task was accomplished by describing gating of GJs using the formalism of the stochastic automata networks (SANs), which allowed for very efficient building and storing of infinitesimal generator of the CTMC that allowed to produce matrices of the models containing a distinct block structure. All of that allowed us to develop efficient numerical methods for a steady-state solution of CTMC models. This allowed us to accelerate CPU time, which is necessary to solve CTMC models, ∼20 times.

#### 1. Introduction

Gap-junctional communication plays an important role in many processes, such as impulse propagation in the heart, communication between neurons and glia, organ formation during early development, regulation of cell proliferation, and metabolic exchange between cells of various tissues, including the lens that lack blood circulation. Gap junction (GJ) channels are formed of connexin (Cx) proteins, which belong to a family of integral membrane proteins exhibiting a tissue specific expression pattern. GJs provide a direct pathway for electrical and metabolic signalling between the cells [1]. In humans, twenty-one isoforms of Cxs form GJ channels [2]. Each GJ channel is composed of two hemichannels (HCs), both oligomerized of six Cxs. Cxs have four alpha helical transmembrane domains (M1 to M4), intracellular N- and C-termini (NT and CT), two extracellular loops (E1 and E2), and a cytoplasmic loop (CL) [3]. Docking of two HCs from neighbouring cells leads to formation of the GJ channel composed of 12 Cxs.

However, despite such complexity, all GJ channels share the same common property—sensitivity to transjunctional voltage (), called voltage gating. Junctional conductance () measured under steady-state conditions decays symmetrically in response to of either polarity, which have been explained by the presence of a -sensitive gate in each opposed HC [4]. Such property, being inherently quantitative, is amenable to the investigation by computational methods.

Earlier, we developed stochastic 4- and 16-state models of voltage gating, containing 2 and 4 gates in series in each GJ channel. In order to demonstrate that the proposed -gating model is adequate, it is necessary to compare its output to experimental results. For example, the proposed stochastic 4- and 16-state models of -gating contained a sizable number (>10) of parameters and for their estimation global optimization (GO) algorithms were successfully used [5, 6]. The simulation of -gating was performed with different sets of parameters. However, for the estimation of a global minimum, it typically requires running thousands of iterations, each lasting for up to 10 seconds, and consequently GO takes tens of hours or days. Thus, the reduction of the computation time that is needed for GO of experimental data is a critical task.

In previous work a discrete time Markov chain (DTMC) model was used [7]. This model described the GJ channel containing 12 gates. In such model, differently from the 4- and 16-state models, it was assumed that each Cx of the GJ channel contains the gate. Since all 12 gates operate at the same time, construction of the transition matrix is not a trivial task. Therefore, transition matrix is dense and when direct methods, that is, Gaussian elimination, are applied, the run-time complexity of steady-state probabilities is in the neighbourhood of . Numerical experiments showed that the use of DTMC model, as opposed to simulation, reduced CPU time ~18- and ~7-fold for 4- and 16-state models, respectively.

When using Markov chain models one needs to build the probability transition matrix and to estimate steady-state probabilities at different s at ~1000 different time moments during a single iteration. Typically, GO of experimental data require to use 500–5000 iterations to find a global minimum. Altogether, this will require to perform ~2 500 000 simulations using Markov model. Thus, it is evident that modeling requires fast construction of the matrix of transition probabilities (transition rates) and fast solution of the steady-state probabilities because the amount of central processing unit (CPU) time is high even at relatively small number of states. In our prior studies [8], we already used continuous time Markov chain (CTMC) model of GJs gating. A transformation of the transition probabilities into transition rates is necessary to generate CTMC model with the same steady-state probabilities, but infinitesimal generator matrix of CTMC model is sparse. For example, infinitesimal generator of GJ model presented in [8] was a tridiagonal matrix. Therefore, generation of matrix and modeling of GJs under steady-states conditions using the CTMC model require smaller amount of CPU time compared to using the DTMC model.

In the present study, we used CTMC for the modeling of GJs containing two voltage-sensitive gates, each of which is composed of six subgates attributed to each Cx; in stochastic 4- and 16-state gating models each gate is regarded as one unit. We also used a stochastic automata network (SAN) formalism for the Markov model specification. SAN formalism allowed accelerating generation of the transition-rate matrices.

In previous study [8] we used piece linear aggregate (PLA) formalism for CTMC model creation. PLA formalism allows building and storing infinitesimal generator of Markov chain model automatically, but matrix structure cannot be easily deduced, especially in more complex models. On the other hand, SAN formalism is a method that is based on using tensor algebra matrix operations. Consequentially, infinitesimal generators have the distinct structure allowing for very efficient application of the numerical methods. Here, we present SAN description of CTMC models of three different GJ models and analyze the efficiency of numerical solution. Since CPU time depends on software and its implementation, we focused on more universal evaluation of a complexity of algorithms. It is based on the exact number of mathematical operations, which is necessary to perform steady-state probabilities calculation.

Our studies showed that if a proper numerical method is applied, then a steady-state solution of the proposed CTMC model of GJs requires 10–20-fold less CPU time compared to DTMC models. We suggest that the use of iterative methods might be especially beneficial in estimation of gating parameters, since it requires repetitive simulations with different sets of parameters. We showed that using a previous solution for evaluating continuous one, when changes are small (<1 mV), allows reducing the number of iterations for at least 30 percent.

#### 2. Methods

##### 2.1. Markov Chain Modeling

We assume the stationary analysis of a homogenous irreducible Markov chain with a finite number of system states, denoted by . Markov chain modeling consists of two stages: (1) construction of transition-rates matrix called an infinitesimal generator and (2) calculation of steady-state probabilities.

The first stage is a model specification. For a GJ gating model this could mean defining the states of a single gate and possible transitions among them, the number of gates in the GJ, and so forth. Basically, it results in formation of a transition matrix , if one assumes a DTMC, or infinitesimal generator matrix , if one assumes CTMC. In this paper we consider mainly the CTMC models.

Formation of and can be performed manually if the size of state space () is relatively small. For larger models the special software can be used, for example, methods based on events language [9], Petri nets [10], and stochastic automata networks [11].

One of the main problems in Markov chain modeling of real systems is a rapid growth of the number of states. The number of states of the Markov chain grows exponentially, when the number of system components grows linearly. Therefore the use of efficient model creation tools and numerical methods is crucial [12].

##### 2.2. Calculation of Steady-State Probabilities

The most difficult and time consuming part of Markov chain modeling is calculation of steady-state probabilities.

Computation of steady-state probabilities of DTMC, which are stored as a row-vector , is the solution of a system of linear equations:

Similarly, computation of steady-state probabilities of DTMC, which are stored as a row-vector **,** is the solution of a system of linear equations:
where denotes a zero row-vector of size .

Equations (1)-(2) can also be interpreted as computations of left side eigenvector corresponding to eigenvalue 1 (in case of DTMC) or eigenvalue 0 (in case of CTMC). Since and are singular, an additional condition is used in all cases.

There are three big classes of algorithms allowing evaluation of steady state probabilities: direct methods, iteration methods, and projection methods. More about numerical methods for solution of general linear systems can be found in [13]; more about application of numerical methods specifically for Markov chains is in [14].

##### 2.3. Stochastic Automata Networks

SAN formalism [15] is one of the most efficient methods used to solve state-space explosion problem, which is very detrimental in Markov chain modeling. SAN allows very efficient construction and storage of infinitesimal generator by using tensor (Kronecker) algebra operations.

Though SAN formalism originally was developed for the modeling of computer networks and communication systems [16, 17], there are multiple examples of SAN use in biology. For example, DeRemigio et al. [18] and Hao et al. [19] used SANs formalism to model calcium channels; Wolf [20] used SANs to describe kinetics of biochemical reactions.

The main idea of SAN formalism is based on the division of a system into smaller subsystems, which can interact among themselves. Those subsystems are described by different stochastic automata. A single automaton is represented by a Markov chain, that is, by the set of subsystem states and possible transitions among them. If two (or more) automata somehow interact among themselves, then transition in one automaton may depend on the state of another one.

The state of the whole system, so called global state, is a compositional state of all automata. An infinitesimal generator matrix of the whole system, so called global generator matrix, can be expressed by infinitesimal generators of individual automata and Kronecker algebra operations. We recall basic definitions below.

Kronecker product of two matrices and is given by

Kronecker sum of two squared matrices and is given by where , are identity matrices of sizes and , respectively. The Kronecker sum of more than two square matrices is also well defined [21].

If the network consists of independent stochastic automata , each governed by infinitesimal generators , , then the global infinitesimal generator can be expressed as a tensor sum:

Expression (5) is also called the SAN descriptor of the system. For a SAN of independent automata steady-state probability, the vector of the whole system is given by where is a steady-state probability vector of an individual automaton .

If we are to consider a network, describing two independent gates that operate between open and closed states with transition rates and , then automata and model each hemichannel/gate, and each of them can be described using the infinitesimal generator:

Thus, the SAN descriptor of the network of two independent gates is given by

If automata in SAN are not completely independent, the interaction among them can also be expressed by the use of Kronecker algebra operations. Plateau expressed two different ways [15] to describe the interaction among automata.

*(1) Functional Transition Rates*. A transition rate in a single automaton may depend on the state of the other automata, that is, on the global system state. Transition rates, which are independent on the global system state, are called constant transition rates.

If we are to consider a network composed of 2 automata and suppose that the transition rates of each connexin depend on the number of connexins/subgates in the open state (denoted by ), then transition rates are functions: and . Consequently, the SAN descriptor of the system is as follows:where denotes the generalized Kronecker product, which deals with functional transition rates [22].

*(2) Synchronizing Events*. Transition in one automaton can cause transitions in other automata. Transition rates are called local, if they are not transition rates of synchronizing events. Synchronizing transitions may also be functional. In this paper, we do not use synchronizing events for the creation of gap junction models.

Plateau and Atif showed that the SAN descriptor of a network consisting of automata and having synchronizing events can be expressed as follows:

The use of SAN formalism basically solves matrix construction problems, since even large matrices can be built and stored (assuming there is enough storage space in operative memory) in a very short time.

The main problem that arises when dealing with SANs of interacting automata is that the steady-state solution cannot be expressed as a simple product form (6). In this case, steady-state probabilities can be found either from solving (2) after is built from the SAN descriptor or directly from the descriptor. That is, building and storing of is not necessary, if special numerical methods are applied. It is possible not to build , since vector SAN descriptor product can be implemented efficiently, for example, by using shuffle algorithm. These problems are considered in detail in [23].

#### 3. Results and Discussion

In this section, we present CTMC models of GJs, created by using SAN formalism. The structure of infinitesimal generators and efficient application of numerical methods for steady-state solutions are considered in detail.

##### 3.1. CTMC Model of the GJ Channel Containing the 12 Two-State Subgates

Gap junctions form clusters (junctional plaques) of individual channels arranged in parallel in the junctional membrane of two adjacent cells. The GJ channel is composed of 2 hemichannels (left and right) arranged in series. Each hemichannel is composed/oligomerized from six Cxs forming a hexamer with the pore inside. We envision that each hemichannel forms the gate, which is composed of six subgates arranged in parallel; that is, to each connexin the subgate is attributed and the GJ channel ultimately contains two gates composed of 12 subgates (see Figure 1).