Biological mechanism and identifiability of a class of stationary conductance model for Voltage-gated Ion channels
Febe Francis, Míriam R. García, Oliver Mason, Richard H. Middleton
BBiological mechanism and identifiability of a class of stationaryconductance model for Voltage-gated Ion channels
Febe Francis , Míriam R. García ∗ , Oliver Mason , and Richard H. Middleton Bioprocess Engineering Group, IIM-CSIC, Vigo, Spain Complex Dynamic Systems & Control, The University of Newcastle, Australia.
November 11, 2018
The physiology of voltage gated ion channels is complex and insightsinto their gating mechanism is incomplete. Their function is best repre-sented by Markov models with relatively large number of distinct statesthat are connected by thermodynamically feasible transitions. On theother hand, popular models such as the one of Hodgkin and Huxleyhave empirical assumptions that are generally unrealistic. Experimentalprotocols often dictate the number of states in proposed Markov mod-els, thus creating disagreements between various observations on thesame channel. Here we aim to propose a limit to the minimum num-ber of states required to model ion channels by employing a paradigmto define stationary conductance in a class of ion-channels. A simpleexpression is generated using concepts in elementary thermodynamicsapplied to protein conformational transitions. Further, it matches wellmany published channel current-voltage characteristics and parametersof the model are found to be identifiable and easily determined fromusual experimental protocols.
1. Introduction
The electrical activity of a living system is a dynamic function of theionic transport across biological membranes. Ion channels are pore-forming protein ensembles that are responsible for the task of regulat-ing ion flows. Gating arises as conformational changes in the proteinsthat comprise the channel. These conformational changes are driven bychanges in the electric field or by molecules (ligands) that bind to them.For this reason, ion-channels are often classified into voltage-gated andligand-gated categories.Voltage-gated ion channels have charged domains that make theirstructure sensitive to variations in the external electric field. For a par-ticular range of membrane potentials, they adopt a conformation with acentral hole: forming a channel for the free movement of ions. Such an‘open’ state is further defined by certain ‘selectivity filters’ (often aminoacids) that would render specificity for the protein [6]. At other mem-brane potentials, the flow of ionic current is blocked as a result of the‘closed’ or ‘inactive’ conformations that the protein adopts. A channelprotein can thus adopt various conformational states with varying de-grees of conductance, and they can spontaneously switch between thesestates. The steady distribution and dynamics of such switches is centralto any study that involves an ion channel.Equation based kinetic models are useful to interpret the behaviourof a channel in a given situation. Starting with the model of Hodgkinand Huxley [15], several researchers have developed theoretical frame-works that partially explain observations made on channel activity. TheHodgkin-Huxley formalism relies on an underlying model of average ∗ E-mail: [email protected] channel conductance and their equations describe the changes of ionicpermeability with membrane potential. The model makes use of hypo-thetical gating particles to bring about the channel’s function, by forc-ing their motion with respect to the electric field across the membrane.Although this model has been used in many instances, hypothetical gat-ing particles do not appear to be consistent with underlying molecularmechanisms.Further developments in the study of ion-channels made an attemptto give a mechanistic description of the gating phenomena. Such mod-els consist of nonlinear ordinary differential equations (ODEs), includ-ing a current balance equation and the dynamics of conformational tran-sitions incorporated as a ‘gating variable’ that corresponds to the stateof the ion channels. Discrete-state Markov models have been used todescribe the different states of the ion channel [20] and have the advan-tage of providing a mechanistic description for the otherwise abstractHodgkin-Huxley formalism based models.Markov chain models are developed on the assumption that ion chan-nels exist in a finite number of significant energy states, with time-homogeneous rates of transition between them. The model consists of atopology of allowed transitions between these states, together with therates for these transitions. Fitting single-channel recording data witha Markovian kinetic scheme has been standard in neurophysiology forquite some time [7, 20, 26, 32]. However, for good agreement with ex-perimental data, frequently the number of closed states needed varieswith the experimental protocol. A significant benefit of Markov mod-els compared to the Hodgkin-Huxley formalism is that the large degreeof freedom in the model stucture allows it to fit more closely with ex-perimental observations [12]. Thermodynamic models are a two stateMarkovian description of channel flipping, the rate kinetics of whichare described by concepts in thermodynamics [10, 29].Fractal models [22, 23] of ion channel gating provide a different de-scription of the underlying mechanism compared to Markov models.Such models are characterised by equations having continuous ratherthan discrete states. The Diffusion models introduced by Millhauser etal. [27], justify Fractal models at a microscopic level. Statistical anal-ysis, however has often favored Markov models over Fractal models[32, 28].A major hurdle in modelling ion-channel gating using a Markov-jump scheme is in the determination of an appropriate number of closedstates. As the topology space expands with the number of states, gener-ating appropriate kinetic schemes can give rise to ambiguity because itmay be possible to come up with multiple schemes that are consistentwith a given set of data. Moreover, numerical simulation of such mod-els becomes time-consuming. Thus the model has limitations from thepoint of view of parameter estimation and further in its use for multi-cellular simulations [12]. The increasing complexity of the topology1 a r X i v : . [ q - b i o . S C ] D ec enerates a need for model reduction. Kienker [20] talks about theexistence of equivalence in topologies of models that are identifiablewithin the same data set. This would imply that models with a largernumber of states are reducible. Keener [19] illustrates the possibility ofreducing the complexity to stable invariant manifolds. This approachwould reduce the dimension of the system without suffering from largeapproximation errors. Furthermore, the time scale of the Markoviantransitions are much faster than the main time scales involved in theaggregate cell voltage and ion behaviour [1, 14], which in turn offersoptions for model reduction.Another difficulty in the acceptability of ion channel models is re-lated to parameter identifiability. A model is said to have structurallyunidentifiable parameters when multiple parameters are equally pow-erful in explaining observed data in noise-free perfect experiments[5, 36, 2]. Since it is well-known that a large segment of ion-channelmodels in the literature lack parameter identifiability in noise experi-ments [12], fundamental studies of structural identifiability are scarcedue to the difficulty of solving the associated symbolic set of equations[8, 9].It should be also noted that in usual experimental protocols optimisedto improve the signal-to-noise ratio [4, 37], provide two sets of sep-arated data: one used to characterise the steady-state, and the otherwith time constants, when dynamics are not sufficiently fast to be disre-garded. With this in mind, Hodgkin-Huxley models employ empiricalexpressions for both steps that are not always realistic [37]. On the otherhand, estimation in models based on Markov chains often make use ofother types of experimental protocols, and do not exploit the large rangeof data available using the standard protocols.In this article, we describe a voltage-gated ion-channel model of sta-tionary conductance with three main characteristics. To begin with, themodel has a clear mechanistic motivation based on the underlying ther-modynamics. In addition, the model is relatively simple to implement,with a small number of easily identifiable parameters. Finally, we testthe model using experimental data from patch clamp and similar studiesreported in the literature. Our model can make use of classical exper-imental protocols and provide a mechanistic formulation for the cal-ibration of steady-state characteristics in the Hodgkin-Huxley settingand, by identifying the form of the stationary distribution, constrainthe estimation process in Markov processes, where the study of globalstructural identifiability is still an open problem. In addition, the simplemodel proposed here can be applied to study the gating of fast ion chan-nels, such as fast persistent and fast activated Na + currents or transientactivated K + current [18].
2. Background
The probability of a protein molecule adapting a particular conforma-tion in space is largely influenced by the electrical force field surround-ing it. When subject to an electric field, charged groups within a proteinwill experience a force and may attain a new electrostatic equilibriumby incorporating angular changes in the dipoles associated with the pep-tide bond [24]. If the thermodynamic kinetic energy is large enough toovercome the energy barrier, the protein takes up a new conformationalstate. Proteins can hence take up a large number of conformationalstates separated by small energy barriers. Despite the continuum ofintermediate states, this dynamical system would typically have only afew stable equilibria [19] and hence a limited number of experimentallyobservable states.
The rate at which the protein switches between such stable states ismainly determined by the driving forces that help in overcoming thethermodynamic energy barrier. However, there are some conforma-tional changes that are not regulated by either electric field or ligands.They may be considered to have a constant energy barrier in the givenenvironment, and hence may be thought to have a uniform rate. Classi-cal thermodynamics identifies the rate of transition between two reac-tion states based on the free energy barrier between them as k = k e − ( (cid:52) G ) / RT , (1)where k is the rate of transition between the two states, k is a constant, ∆ G is the free energy barrier between the two states, R is the universalgas constant and T the absolute temperature.More generally, the free energy barrier may be dependent on the elec-tric field, that is, the membrane potential V (figure 1). In this case wehave k ( V ) = k e − ∆ G ( V ) / RT (2)Here ∆ G(V) is the free energy barrier between the two states defined bythe voltage [26]. V ε c ε c ε c ε o ε o ε o V V C O k o k c Figure 1:
Energy profile for the transitions between open and closed states in anion channel, modulated by a change in membrane potential: a changein membrane potential brings in protein conformational changes andthereby free energies of protein conformation. The rates of channelclosing and opening, k c and k o , are a function of the free energy bar-riers ( ε c and ε o ) defined by structural configurations of the respectivestates. The figure describes how the open probability of a channel isenhanced by a change in the ambient voltage from V to V by virtueof its conformational energy. In the case of voltage based transition, the activation energy can beexpressed in a general form by using a Taylor series expansion [10, 29],as follows: ∆ G ( V ) = a + bV + cV + ··· and the rate of transition k(V) may be written as, k ( V ) = k e − ( a + bV + cV + ··· ) / RT (3)Here a corresponds to the free energy independent of the electric fieldand bV corresponds to interactions between the electric field and iso-lated charges and rigid dipoles on the protein. The higher order termscorrespond to the influence of polarization and deformation within theprotein structure as well as mechanical constraints. These effects areusually negligible as the trans-membrane voltage variations are gener-ally small [10].In what follows, we are mainly interested in models where the freeenergy barrier is linear in the membrane potential. The effect of temper-ature variations may be neglected for the system under consideration,as long as the physiological environment remains unaffected. In this2ase, the following trivial lemma, helps us derive a simple expressionfor voltage regulated conformational transitions. Lemma 1.
Consider a set { ∆ G ( V ) ,..., ∆ G M ( V ) } of activation en-ergies, where each ∆ G i ( V ) is affine in V. Denote the correspondingtransition rates by { k ,..., k M } . Then any expression of the form ∏ Pp = (cid:16) k ip ( V ) k jp ( V ) (cid:17) may be expressed as e [( V − V h ) s ] The series of molecular changes associated with the opening of an ionchannel is often described using Markov models with discrete states[20, 32, 19]. To begin, consider the simple case of a transition betweenan open and closed state; let O and C represent the probabilities thatthe molecule is in the corresponding open and closed state at a giventime. Since the equation governing changes in probabilities for a singlemolecule has a form similar to the rate equation for a large number ofmolecules, the transition may be represented by a kinetic scheme, asfollows: C k OC (cid:47) (cid:47) O k CO (cid:111) (cid:111) , where k i j represents the rate of transition from state j to state i . Theabove kinetic scheme has a transition intensity matrix K (cid:20) − k CO k OC k CO − k OC (cid:21) and the steady state probability for the channel to be in an open stateleads O = + k OC k CO . Using Lemma (1) we obtain, O = + e [( V − V h ) s ] (4)The expression is the modified Boltzmann’s expression used in model-ing ion-channel gating [6, 17, 39]. This sigmoidal function of voltage issymmetric about the half-activation voltage, V h . However, experimentaldata frequently show one or both of: (i) activation and inactivation be-haviour (ii) asymmetric behaviour, and hence it is hard to fit data witha single Boltzmann function. This implies that the two-state Markovchain with linear energy barriers does not ably model experimental ob-servations.Three possible ways to resolve this are (i) model the system as havingmore than two macro-states (stable conformational states), (ii) modelthe system as having transition rates that change with the dwell times,leading to a fractal model; or, (iii) use non-linear terms in the energyexpression.We focus on the first option here for the following reasons. The sec-ond option (ii) would give rise to time-inhomogeneous Markov chains,which are typically very difficult to analyse. Moreover, the identifi-cation task becomes less tractable in this case. As mentioned previ-ously, the last option (iii) includes higher order terms that are usuallyneglected. Also, as we shall see later in section (6) the use of higher or-der non-linear energy models does not appear to result in simpler mod-els, nor does it give more accurate fits to experimental data.In fact, usually models developed on the basis of Markovian dynam-ics have more than three states [option (i)] for a good agreement withexperimental observations [20, 26]. In such models the rate constant isoften defined with an exponential or a Boltzmann function in voltage,both obtained from a linear approximation of equation (3). We will now analyse how best we can extend the two-state model byoption (i) to have a network with limited and identifiable macro-states.
3. A Multiple Conformation Extension of the‘Modified Boltzmann Function’ for stationaryconductance
Here we consider the transition network of the ion-channel as a systemwith a set of N stable states, marked i = , ,..., n ; with n = N − S i denotes the probability for the protein to bein state i at any time t , the system obeys the master equation ˙ S = K S , (5)where S = ( O , S , S ,..., S n ) and K ∈ R N × N is a transition matrix with k i j ≥ j to state i . The diagonalelements satisfy, k ii = − ∑ i (cid:54) = j k i j , to ensure the system evolves on theprobability simplex. The entry k i j (cid:54) = j to state i . The stationary probability distribution S satisfies K S = T N S =
1, where N is the column vector of size N with allentries equal to one.As is standard, we associated a directed graph with the Markov pro-cess described by (5) consisting of the nodes { , ,..., n } with an edgefrom state j to state i ( i (cid:54) = j ) if and only if k i j (cid:54) = O takes a particularly simple form. Wealso describe a general condition on the structure of the graph associatedwith the Markov process that is sufficient for this simple form to hold. A three state transition diagram for channel opening is given below C k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:47) (cid:47) O k (cid:111) (cid:111) . Here O represents an open state and C and C represent closed or in-active conformations. The steady-state probability for the channel to bein the open state can be calculated as O = (cid:16) + k k + k k k k (cid:17) . By using lemma 1, there exist V h , s , V h , s such that O = + e [ ( V − V h ) s ] + e [ ( V − V h ) s ] . (6)In the case of a slightly different topology, C k (cid:47) (cid:47) O k (cid:111) (cid:111) k (cid:47) (cid:47) C k (cid:111) (cid:111) the probability of the channel to be in an open conformation would be3 = (cid:16) + k k + k k (cid:17) and once again, using lemma 1, the stationary probability of being opencan be represented in the form (6). This scheme can be generalized for n macrostates depending on theposition of the open state in the entire topology. A topology involvingthe open state on the network extremum, O k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:47) (cid:47) C k (cid:111) (cid:111) C i C n − k n , n − (cid:47) (cid:47) C nk n − , n (cid:111) (cid:111) would have an open state probability O = (cid:32) + n ∑ j = j ∏ i = k i , i − k i − , i (cid:33) and a topology C k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:47) (cid:47) ··· C p − k (cid:111) (cid:111) k , p − (cid:47) (cid:47) O k p − , (cid:111) (cid:111) k p + , (cid:47) (cid:47) ··· C n − k , p + (cid:111) (cid:111) k n , n − (cid:47) (cid:47) C nk n − , n (cid:111) (cid:111) would yield O = (cid:32) + p − ∑ j = j ∏ i = k i , i + k i + , i + n ∑ j = p j ∏ i = k i + , i k i , i + (cid:33) . In either case, with respect to the argument in lemma (1), the open stateprobability can be reduced and in general, a system with N macro-statesrelated in order of their conformational transitions would have an openstate probability O = + n ∑ i = e ( V − V h , i ) s i , where n = N −
1, is the number of transitions in the linear network.
Remark 1.
For the networks considered so far, rendering one of thetransitions irreversible, would make the network absorbing in nature.In other words, the system would reach and never leave a fixed confor-mation. Since such a possibility is not physically reasonable for ion-channels under consideration, this case is not considered further.
More generally, consider a transition network in which a unique simplepath exists from every stable state to the open state; so in the directedgraph associated with the matrix K , there is a unique path from every j (cid:54) = K is irreducible[16].The celebrated Markov Chain Tree Theorem [21] allows us to char-acterise the form of the steady state probability of the open state in thiscase. This result is usually stated for column stochastic matrices or dis-crete Markov chains; however, it is trivial to see that an exact analoguealso holds for continuous chains with matrices of the form K . We state arestricted version of this result below but first introduce some notation.For the directed graph G associated with K , a rooted spanning tree T i at i ∈ { ,..., n } consists of the vertices { , ,..., n } and has the follow-ing properties: (i) T i is acyclic; (ii) for every j (cid:54) = i , there exists exactly one outgoing edge from j ; (iii) there exist no edge outgoing from i .We denote by w ( T i ) the weight of T i , which is given by the product ofthe entries of K corresponding to the edges in T i . For 0 ≤ i ≤ n , let T i denote the set of all directed spanning trees rooted at i , and define w i = ∑ T i w ( T i ) . Note that each w i will be a sum of terms of the form k i j k i j ··· k i n j n (7)The following result is now a simple re-wording of the Markov ChainTree Theorem as presented in [21] and elsewhere. Theorem 1.
Assume the matrix K is irreducible. The unique stationaryprobability vector associated with K , π is given by π i = w i ∑ j w j where w i is defined as above for ≤ i ≤ n. If there is a unique path from every node j (cid:54) = T rooted at 0 (which represents the open state). It then follows thatthe steady state probability of the channel being in the open state is ofthe form O = w ( T ) ∑ j w j . (8)As there is only a single term of the form (7) in the numerator, it followsreadily by combining (8) with Lemma 1 that O can will take the form O = + N ∑ i = e ( V − V h , i ) s i . (9)A few examples of network topologies for which this form is guar-anteed by this analysis are illustrated in table 1.
4. Numerical analysis and low order approximation
We have observed so far that open-state probabilities of ion-channelsmay be expressed in a similar form to the modified Boltzmann equation,but with sum of exponentials replacing the single exponential term asin equation (9). The value of N is generally not less than the numberof transition macro-states. The ambiguity in the value of N (which is atlarge dependent on the network structure) together with computationaland identifiability issues for large N motivate us to consider modelswith small N . Of course, when considering such approximations, wepotentially open up inaccuracies in the model [23], and it is thereforeimportant to check that any such reduction still accurately captures theobservable behaviour of the system.A good way of approximating the model would be to search for theminimum number of exponential terms that yields a good fit for ex-perimentally observed ion channel current-voltage characteristic data.Let us denote the M -vector of ionic currents obtained from patch clampexperiments on a certain channel by I m ∈ R M . Let V m ∈ R M be thecorresponding voltage vector. In order to fit these data to equation (9),the distance between the experimental data and the model needs to beminimized. The distance may be measured as square of the L norm: J ( I m , I ( V m )) = (cid:107) I m − I ( V m ) (cid:107) = m ∑ j = (cid:16) I mj − I ( V mj ) (cid:17) where I ( V m ) ∈ R M is the vector of ionic current predicted by the modelgiven membrane potentials in the vector V m . I mj and V mj are respec-tively the j th position in vectors I m and V m . Let g max be the maximalconductance of the channel and V ∗ , the equilibrium potential of the ion4 k (cid:32) (cid:32) C k (cid:126) (cid:126) O k (cid:96) (cid:96) k (cid:62) (cid:62) k (cid:15) (cid:15) C k (cid:79) (cid:79) C k (cid:15) (cid:15) C k (cid:15) (cid:15) C k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:79) (cid:79) k (cid:32) (cid:32) C k (cid:79) (cid:79) k (cid:47) (cid:47) k (cid:126) (cid:126) C k (cid:111) (cid:111) O k (cid:96) (cid:96) k (cid:62) (cid:62) k (cid:15) (cid:15) C k (cid:126) (cid:126) k (cid:32) (cid:32) k (cid:79) (cid:79) C k (cid:62) (cid:62) C k (cid:96) (cid:96) O k (cid:126) (cid:126) k (cid:32) (cid:32) C k (cid:62) (cid:62) C k (cid:96) (cid:96) k (cid:126) (cid:126) k (cid:32) (cid:32) C k (cid:62) (cid:62) k (cid:126) (cid:126) k (cid:32) (cid:32) C k (cid:96) (cid:96) C k (cid:62) (cid:62) C k (cid:96) (cid:96) S t a r t opo l ogy , O = + ∑ i = e ( V − V h , i ) s i E x t e nd e d s t a r t opo l ogy , O = + ∑ i = e ( V − V h , i ) s i T r ee t opo l ogy , O = + ∑ i = e ( V − V h , i ) s i O k (cid:127) (cid:127) k (cid:31) (cid:31) C k (cid:63) (cid:63) k (cid:47) (cid:47) C k (cid:111) (cid:111) O k (cid:15) (cid:15) k (cid:47) (cid:47) C k (cid:15) (cid:15) C k (cid:79) (cid:79) k (cid:47) (cid:47) C k (cid:79) (cid:79) k (cid:111) (cid:111) O k (cid:15) (cid:15) k (cid:47) (cid:47) C k (cid:15) (cid:15) C k (cid:79) (cid:79) k (cid:47) (cid:47) C k (cid:111) (cid:111) O = + ∑ i = e ( V − V h , i ) s i O = + ∑ i = e ( V − V h , i ) s i O = + ∑ i = e ( V − V h , i ) s i R i ng t opo l og i e s w it h a un i qu e s i m p l e p a t h t o t h e op e n s t a t e O k (cid:127) (cid:127) k (cid:31) (cid:31) C k (cid:126) (cid:126) C k (cid:63) (cid:63) k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:62) (cid:62) O k (cid:127) (cid:127) k (cid:31) (cid:31) C k (cid:126) (cid:126) C k (cid:47) (cid:47) C k (cid:95) (cid:95) k (cid:111) (cid:111) k (cid:62) (cid:62) O k (cid:15) (cid:15) k (cid:47) (cid:47) C k (cid:47) (cid:47) C k (cid:15) (cid:15) k (cid:111) (cid:111) C k (cid:79) (cid:79) k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:47) (cid:47) C k (cid:111) (cid:111) k (cid:79) (cid:79) O = + ∑ i = e ( V − V h , i ) s i O = + ∑ i = e ( V − V h , i ) s i O = + ∑ i = e ( V − V h , i ) s i M e s h t opo l og i e s w it h a un i qu e s i m p l e p a t h t o t h e op e n s t a t e T a b l e : E x a m p l e s o fi on - c hann e lt opo l og i e s t ha t s a ti s f y t h ec r it e r i ao f aun i qu e s i m p l e pa t h t o t h e op e n s t a t e { V h , i } Ni = , { s i } Ni = , g max J ( I m , I ( V m )) (10)subject to the conductance expression: I ( V j ) = g max ( V j − V ∗ ) + N ∑ i = e ( V − V h , i ) s i ∀ j = ,.., m (11)and with bounds that are physiologically justified.This optimization problem brings forth two significant difficulties:(i) the cost function is not convex and (ii) the number of dominant con-formations that the protein adopts is usually unknown and hence theorder of the exponentials, N is unknown. To guarantee that the globaloptimum is achieved, an initial guess should be carefully selected. It ispossible to have a calculated guess for the V h values from the current-voltage characteristics, positioned centrally over the range of values atwhich the curve shows steady activation or inactivation. The steepnessof the tangent drawn at this estimated voltage may be used as an ini-tial estimate of the slope value. Alternatively, the recently developedglobal optimization Scatter Search based methodology, SSm GO [11],can be used. This algorithm combines a population-based metaheuristicmethod with a local optimization.Regarding the number of exponentials, the aim is to curtail the num-ber of conformational transitions so that it leads to minimum parametersto fit the data. For this purpose an algorithm was implemented whichcalculates the global optimum for several order of exponentials untilthe confidence in the fit is accomplished. The algorithm starts by N = ε times the original data. The measureconsidered is again the square of the L norm and mathematically, thiscriterion may be formulated as: J ( I m , I ( V m )) < (cid:107) ε I m (cid:107) For all the cases we have considered, a good fit is obtained by setting ε = . N = N =
1) for simplesymmetric data-plots, and for all other cases, a three dimensional statespace ( N =
2) could accurately explain the data. This would imply thatto a large extent stationary distributions of voltage gated ion-channelscan be adequately modeled as dwelling predominantly among three dif-ferent macro-states even if they are able to change between several con-formations. The time that the protein dwells in some of these states maybe considerably smaller than that of the three dominant conformations.The three macrostates may incorporate some of the effects of the mi-nor conformations rather than completely negating their influences; asreflected from the data-fits.With the above argument we propose the following equation for theopen state probability corresponding to a three state system, O = + e [ ( V − V h ) s ] + e [ ( V − V h ) s ] (12)as a good approximation to represent the stationary probability of volt-age gated ion-channels to remain open. The channel current may hencebe calculated with the following equation: I i = g max , i ( V − V ∗ i ) + e { V − V h , ( i ) } s ( i ) + e { V − V h , ( i ) } s ( i ) (13)where, I i is the channel current; V ∗ i , the reversal potential of ion i ; V , the membrane potential; g max , i , the maximal conductance of the channelfor the ion i and V h , ( i ) , s ( i ) , V h , ( i ) , s ( i ) being the parameters for theactivation function as defined by equation (12).
5. Fitting of experimental data to the model
Most papers in the literature follows the voltage step protocol wheresteady-state voltage-current characteristics are obtained independentlyof the ion channel dynamics. Several were selected for utilizable dataon single channel studies. We collect data from both, stationary conduc-tance obtained from peak currents (with less signal-to-noise ratio [3])or from raw steady-state measurements. Data were extracted from thepublished curves using the Enguage Digitizer 4.1. The digitized datasets were fit by the equation (13). The global scatter search algorithm,SSm GO [11] implemented in MATLAB was used for making the fits.The resulting data fits are presented in figures 2,3,4 and summarized intable 2. In all cases, with N ≤
2, we are able to represent the observeddata well using the model structure proposed.
6. Comparison with non-linear thermodynamicmodels
The estimation of the open-state probability of voltage-gated ion chan-nels to an equation of form as (9) would rather be equivalent to a steady-state approximation of the Hodgkin-Huxley formalism. However, as weseek to look beyond the Hodgkin-Huxley formalism for a good mech-anistic description, the given model need to be compared with effortsmade to represent the system on a realistic perspective. Thermody-namic formalism comes closest to this effort. Although the fundamen-tals are in place, such models [10] show very poor performance in evensmall extrapolations from the given data (Figure 5). We would alsoargue that the biophysical basis of the multiple conformation model ismuch clearer than that of a Taylor’s series dual conformation model.Ozer [29, 30] modified the non-linear model to a functional form bylumping the different transitions in the protein to a single event. Themodel which uses a sum of Gaussion distribution, seems to give thebest fits for experimental data. The steady-state open probability ofion-channels written with respect to this model, may be arrived usingequation (8) of [30] as, O = + n ∑ i = α , i e [( V − V α , i ) s α , i ] n ∑ i = β , i e [( V − V β , i ) s β , i ] (14)where, n was defined as the number of distinct transitions defined bydifferent energy barriers. It should be noted that in work [30] was alsoable to get acceptable fits with two transitions. The problem with thismodel however seems to be in the existence of unidentifiable parame-ters, as is the case with majority of the existing Markov Models [12].Further, if no manipulation is carried out to collect group of parametersthat are no identifiable only with stationary data, the model asks for fit-ting a minimum of twelve parameters (provided two macro-transitions( n =
2) gives acceptable fits).A mathematical model for a biological phenomenon is assured tohave uniquely estimated parameters if the model structure ensures iden-tifiability [36]. Conventionally, ion-channel data has limited interpre-tations by producing indistinguishable models [20]. In-silico represen-tation of ion-channel dynamics is yet to come up with a model for thereason that the model structure needs to be tweaked each time to incor-porate experimental observation. This indeed is a question of parameteridentifiability. The model that we have presented in this article, is found6
80 −60 −40 −20 0 20 40 60−1−0.8−0.6−0.4−0.20 Membrane Potential (mV) N o r m a li ze d C u rr e n t FitData (a) Normalized current-voltage relationship of the T-type calcium channel from thedigitized data of Talavera and Nilius (2006) [35] −60 −40 −20 0 20 40 60−1−0.8−0.6−0.4−0.20 Membrane Potential (mV) N o r m a li ze d C u rr e n t FitData (b) Normalized current-voltage relationship of the wild-type calcium channel [6] −60 −40 −20 0 20 40−0.2−0.15−0.1−0.050 Membrane Potential (mV) C u rr e n t ( µ A ) FitData (c) Averaged, peak current–voltage plots of Ca v . α L-type calcium channels [39] −80 −60 −40 −20 0 20 40−0.25−0.2−0.15−0.1−0.050 Membrane Potential (mV) C u rr e n t ( µ A ) FitData (d) Averaged, peak current–voltage plots of Ca v . α L-type calcium channels [39]
Figure 2:
Current-Volatge relationship of Voltage-gated Ca + channels fitted with equation (12) . Voltage-gated Ca + channels (VGCC) can be classified into high-voltage-activated (HVA) and low-voltage-activated (LVA) channels, which implies that LVA channels activate at 20 to 30 mV more negative potentialsthan HVA channels. T-type calcium channels are LVA and show fast macroscopic inactivation where as L-Type calcium channels are HVAs. −80 −60 −40 −20 0 20 40 60−800−700−600−500−400−300−200−1000 Membrane Potential (mV) C u rr e n t ( p A ) FitData (a) Peak inward current vs. voltage relation for the current records of sodium channel[34] fitted with equation (12) −60 −40 −20 0 20 40 60−1−0.8−0.6−0.4−0.20 Membrane Potential (mV) N o r m a li ze d P ea k C u rr e n t FitData (b) Normalized current-voltage relationship of sodium channel from the macropatchcurrent [31] fitted with equation (12)
Figure 3:
Current-Volatge relationship of Voltage-gated Na + channels fitted with equation 12.
100 −50 0 50−2024681012 Membrane Potential (mV) C u rr e n t ( µ A ) FitData (a) KCNH5 −100 −80 −60 −40 −20 0 20 40 60 80−0.500.511.522.53 Membrane Potential (mV) C u rr e n t ( µ A ) FitData (b) KCNH7 −60 −40 −20 0 20 4000.511.522.533.5 Membrane Potential (mV) C u rr e n t ( µ A ) FitData (c) KCNB1
Figure 4:
Current-voltage relations for KCNH5, KCNH7 and KCNB1 class of potassium channels expressed in oocytes [40] fitted with equation 12. −100 −50 0 50 100−1−0.9−0.8−0.7−0.6−0.5−0.4−0.3−0.2−0.100.1 Membrane Potential (mV) N o r m a li ze d C u rr e n t ExperimentalNon−linear thermodynamic modelOur model (a) Calcium channel[35] −100 −50 0 50 100−0.500.511.522.5 Membrane Potential (mV) C u rr e n t ( µ A ) ExperimentalNon−linear thermodynamic modelOur model (b) Potassium channel[40]
Figure 5:
Current-voltage relationship of ion-channels obtained from the studies [35, 40] fitted with a third order non-linear model and modified model presentedin this article according to equation (12) to have structurally output locally identifiable (s.o.l.i) parameters (seeA for formal definition), according to the following lemma.
Lemma 2.
The model for ion channel open state probability describedby means of a Markov chain with three macro-states [equation (12)] iss.o.l.i by the Taylor approximation of fourth order.Proof: see A
7. Conclusion
In this article, we bring forth an important criterion in modelling ionchannel activity using Markov jump schemes. We propose that the be-haviour of the system may be developed by an analytical design thatis simple and adaptable; and towards this, a minimum of three Markov8igure Reference Channel type V h ( mV ) s ( mV ) V h ( mV ) s ( mV ) g max Ca v . mV − )2(b) Beyl et al. (2007) [6] Ca v . mV − )2(c) Xu and Lipscombe (2001) [39] Ca v . µ A / mV )2(d) Xu and Lipscombe (2001) [39] Ca v . µ A / mV )3(a) StÃijhmer et al. (1987) [34] Na v . pA / mV )3(b) Ruben et al. (1997) [31] Na v . a -25.6539 -0.14232 38.2236 0.088449 0.0138 ( mV − )4(a) Zou et al. (2003) [40] K v . µ A / mV )4(b) Zou et al. (2003) [40] K v . µ A / mV )4(c) Zou et al. (2003) [40] K v . µ A / mV ) Table 2:
Fitted parameters of Channel current - membrane potential data of a few ion-channels using equation 12 states are necessary to model the dynamics.To illustrate the criterion, we describe the physical transitions ofa particular class of voltage gated ion channel proteins by means offree energy changes associated with perturbations in its environmen-tal electric field from a thermodynamic perspective. For this class ofion-channels we assume the existence of a unique simple path of tran-sition, from every stable state to the open state. The resulting model isa simple kinetic expression (equation 12) which approximates the sta-tionary open probability of the channel for a given membrane potentialand is a generalized form of the single exponent modified Boltzmann’sfunction.The model proposed in this article, is identifiable and requires smallnumber of parameters compared to existing models. This would in turnreduce the computational cost associated with regular Markov Mod-els, without compromising much from the mechanistic viewpoint. Itsmain limitation is that it has been developed for steady-state data andmay not represent ion-channel dynamics in its entirety but only fastion channels. Nonetheless, when dynamics are relevant, the proposedscheme may be extended and the resulting model calibrated by usingvoltage-current data from classical experimental protocols.
Acknowledgment
The authors are grateful to the Science Foundation of Ireland forfunding the research (Science Foundation of Ireland Research Grant07/PI/I1838). We are also indebted to Professor Peter Wellstead andWilhelm Huisinga for their valuable support and helpful discussions.
References [1] David J. Aidley and Peter R. Stanfield.
Ion channels: moleculesin action . Cambridge University Press, 1996.[2] E. Balsa-Canto, J. R. Banga, and M. R. García. Dynamic modelbuilding using optimal identification strategies, with applicationsin bioprocess engineering.
Process Systems Engineering , pages441–467, 2010.[3] J. Beaumont, FA Roberge, and DR Lemieux. Estimation ofthe steady-state characteristics of the hodgkin-huxley model fromvoltage-clamp data.
Mathematical biosciences , 115(2):145–186,1993.[4] J. Beaumont, FA Roberge, and LJ Leon. On the interpretation ofvoltage-clamp data using the hodgkin-huxley model.
Mathemati-cal biosciences , 115(1):65–101, 1993.[5] R. Bellman and KJ Åström. On structural identifiability.
Mathe-matical Biosciences , 7(3):329–339, 1970. [6] Stanislav Beyl, Eugen N. Timin, Annette Hohaus, Anna Stary,Michaela Kudrnac, Robert H. Guy, and Steffen Hering. Probingthe Architecture of an L-type Calcium Channel with a ChargedPhenylalkylamine.
J. Biol. Chem. , 282(6):3864–3870, 2007.[7] D. Colquhoun and A. G. Hawkes. On the stochastic properties ofsingle ion channels.
Proc. R. Soc. B-Biol. Sci. , 211(1183):205–235, 1981.[8] D. Csercsik, Szederkényi G., and K.M. Hangos. Identifiabilityof a Hodgkin-Huxley Type Ion Channel under Voltage Step Mea-surement Conditions.
Proceedings of the 9th International Sym-posium on Dynamics and Control of Process Systems (DYCOPS2010) , MoAT3.4:318–323, 2010.[9] D. Csercsik, K.M. Hangos, and G. Szederkényi. Identifiabilityanalysis and parameter estimation of a single hodgkin-huxley typevoltage dependent ion channel under voltage step measurementconditions.
Neurocomputing , 77:178–188, 2012.[10] A. Destexhe and J.R. Huguenard. Nonlinear thermodynamic mod-els of voltage-dependent currents.
Journal of Computational Neu-roscience , 9:259–270(12), 2000.[11] Jose A Egea, Rafael Martí, and Julio R Banga. An evolutionarymethod for complex-process optimization.
Computers & Opera-tions Research , 37(2):315–324, 2010.[12] Martin Fink and Denis Noble. Markov models for ion channels:versatility versus identifiability and speed.
Philosophical Trans-actions of the Royal Society A: Mathematical, Physical and Engi-neering Sciences , 367(1896):2161–2179, 2009.[13] Febe Francis, Míriam R García, and Richard H Middleton. A sin-gle compartment model of pacemaking in dissasociated substantianigra neurons.
Journal of computational neuroscience , pages 1–22, 2013.[14] Bertil Hille.
Ion Channels of Excitable Membranes . Sinauer As-sociates, Inc., 2001.[15] A. L. Hodgkin and A. F. Huxley. A quantitative description ofmembrane current and its application to conduction and excitationin nerve.
J Physiol , 117(4):500–544, August 1952.[16] Roger A. Horn and Charles R. Johnson.
Matrix Analysis . Cam-bridge University Press, 1985.[17] Shuyun Huang, Qing Cai, Weitian Liu, Xiaoling Wang, and TaoWang. Whole-cell recordings of voltage-gated calcium, potassiumand sodium currents in acutely isolated hippocampal pyramidalneurons.
Journal of Nanjing Medical University , 23(2):122 – 126,2009.918] E. M. Izhikevich.
Dynamical systems in neuroscience: the geome-try of excitability and bursting . Computational neuroscience. MITPress, 2007.[19] James Keener. Invariant manifold reductions for markovian ionchannel dynamics.
Journal of Mathematical Biology , 58:447–457,2009.[20] P. Kienker. Equivalence of Aggregated Markov Models of Ion-Channel Gating.
Proceedings of the Royal Society of London. B.Biological Sciences , 236(1284):269–309, 1989.[21] F. Leighton and R. Rivest. Estimating a probability using finitememory.
IEEE Transactions on Information Theory , 32(6), 1986.[22] L S Liebovitch, J Fischbarg, J P Koniarek, I Todorova, andM Wang. Fractal model of ion-channel kinetics.
Biochim Bio-phys Acta. , 896(2):173–180, Jan 1987.[23] L.S. Liebovitch. Testing fractal and markov models of ion channelkinetics.
Biophysical Journal , 55(2):373 – 377, 1989.[24] Smarajit Manna, Jyotirmoy Banerjee, and Subhendu Ghosh.Breathing of voltage dependent anion channel as revealed by thefractal property of its gating.
Physica A: Statistical Mechanicsand its Applications , 386(1):573 – 580, 2007.[25] Satoshi Matsuoka, Hikari Jo, Masanori Kuzumoto, AyakoTakeuchi, Ryuta Saito, and Akinori Noma. Modeling energet-ics of ion transport, membrane sensing and systems biology of theheart.
Molecular System Bioenergetics: Energy for Life , pages435–455, 2007.[26] LS Milescu, G Akk, and F Sachs. Maximum likelihood estimationof ion channel kinetics from macroscopic currents.
BiophysicalJournal , 88(4):2494–2515, APR 2005.[27] G L Millhauser, E E Salpeter, and R E Oswald. Diffusion mod-els of ion-channel gating and the origin of power-law distribu-tions from single-channel recording.
Proceedings of the NationalAcademy of Sciences of the United States of America , 85(5):1503–1507, 1988.[28] Ali Nekouzadeh and Yoram Rudy. Statistical properties of ionchannel records. part ii: Estimation from the macroscopic current.
Mathematical Biosciences , 210(1):315 – 334, 2007.[29] Mahmut Ozer. An improved non-linear thermodynamic model ofvoltage-dependent ionic currents.
Neuroreport , 15(12):1953–7,2004.[30] Mahmut Ozer. A comparative analysis of linear, nonlinear andimproved nonlinear thermodynamic models of voltage-dependention channel kinetics.
Physica A: Statistical Mechanics and itsApplications , 379(2):579 – 586, 2007.[31] Peter C. Ruben, Andrea Fleig, David Featherstone, John G.Starkus, and Martin D. Rayner. Effects of clamp rise-time on ratbrain iia sodium channels in xenopus oocytes.
Journal of Neuro-science Methods , 73(2):113 – 122, 1997.[32] M.S. Sansom, F.G. Ball, C.J. Kerry, R. McGee, R.L. Ramsey, andP.N. Usherwood. Markov, fractal, diffusion, and related modelsof ion channel gating. a comparison with experimental data fromtwo ion channels.
Biophysical Journal , 56(6):1229 – 1243, 1989. [33] Jonathan R Silva, Hua Pan, Dick Wu, Ali Nekouzadeh, Keith FDecker, Jianmin Cui, Nathan A Baker, David Sept, and YoramRudy. A multiscale model linking ion-channel molecular dynam-ics and electrostatics to the cardiac action potential.
Proceed-ings of the National Academy of Sciences , 106(27):11102–11106,2009.[34] W. StÃijhmer, C. Methfessel, B. Sakmann, M. Noda, andS. Numa. Patch clamp characterization of sodium channelsexpressed from rat brain cdna.
European Biophysics Journal ,14(3):131–138, 1987.[35] Karel Talavera and Bernd Nilius. Biophysics and structure-function relationship of t-type ca2+ channels.
Cell Calcium ,40(2):97 – 114, 2006. T-type calcium channels: from old physiol-ogy to novel functions.[36] E. Walter and L. Pronzato.
Identification of Parametric Modelsfrom Experimental Data . Springer, 1997.[37] A. R. Willms, D. J. Baro, R. M. Harris-Warrick, and J. Guck-enheimer. An improved parameter estimation method forhodgkin-huxley models.
Journal of Computational Neuroscience ,6(2):145–168, 1999.[38] Inc. Wolfram Research.
Mathematica . Wolfram Research, Inc.,version 7.0 edition, 2008.[39] Weifeng Xu and Diane Lipscombe. Neuronal CaV1.3alpha1 L-Type Channels Activate at Relatively Hyperpolarized MembranePotentials and Are Incompletely Inhibited by Dihydropyridines.
J. Neurosci. , 21(16):5944–5951, 2001.[40] Anruo Zou, Zhixin Lin, Margaret Humble, Christopher D. Creech,P. Kay Wagoner, Douglas Krafte, Timothy J. Jegla, and Alan D.Wickenden. Distribution and functional properties of humanKCNH8 (Elk1) potassium channels.
Am J Physiol Cell Physiol ,285(6):C1356–1366, 2003.
A. Model Identifiability, Proof of Lemma 2
Simplicity and identifiability of ion channel model would be of ad-vantage especially, when such models are employed to develop largerframeworks describing metabolic activity of a cellular system, such apacemaking neurons [13] or cardiac cells [25]. To study the parame-ter identifiability of model given by equation (12), let us consider thefollowing general mathematical structure: M ( θ ) : y ( θ ) = f ( θ , x ) where θ ∈ R n θ is the set of parameters to be estimated and y ( θ ) ∈ R and x ∈ R denote the output and input of the system, respectively. Inaddition, let us denote the search space of the parameter by Θ ⊆ R n θ . Definition 1.
The model M ( θ ) is said to be structurally output globallyidentifiable (s.o.g.i) , if for any (cid:101) θ ∈ Θ , except for the points of a subsetof measure zero, and for all x ∈ R :y ( θ , x ) ≡ y ( (cid:101) θ , x ) = ⇒ θ ≡ (cid:101) θ If the same conditions is fulfilled only in the neighborhood of θ then theparameters will be structurally output locally identifiable (s.o.l.i.) This definition is quite general and is often difficult to be put intopractice. One of the usual approach in such instance is the so-called"Taylor approach" [36] where both outputs y ( θ , x ) and y ( (cid:98) θ , x ) are ap-proximated by a Taylor expansion around a point x ∗ ∈ R . If we denote10y U nx ∗ the ball around x ∗ , where the Taylor approximation of order n holds; the general identifiability definition is therefore set as: Definition 2.
The model M ( θ ) is said to be s.o.g.i in x ∈ U nx ∗ if thesystem of equations: (cid:20) ∂ i y ( θ , x ) ∂ x i (cid:21) x = x = (cid:34) ∂ i y ( (cid:101) θ , x ) ∂ x i (cid:35) x = x ∀ i = ,... n (15) has an unique solution. In the same way, the model will be s.o.l.i inx ∈ U nx ∗ if a finite number of solution is obtained.Proof. It has to be proved that, the ion channel open state probabilityrepresented by the equation, G ( θ ) = + e ( V − V h ) s + e ( V − V h ) s (16)is s.o.l.i for any V ∈ U nV ∗ , by using the Taylor approximation of fourthorder. It may be noted that studying the identifiability property for G ( θ ) with parameters θ = [ V h , , s , V h , , s ] ∈ R (cid:101) θ = [ (cid:101) V h , , (cid:101) s , (cid:101) V h , , (cid:101) s ] ∈ R is equivalent to studying this property for the model M ( θ ) = e ( aV + b ) + e ( cV + d ) , θ = [ a , b , c , d ] where a = s , b = − V h s , c = s , d = − V h s (cid:101) a = (cid:101) s , (cid:101) b = − (cid:101) V h (cid:101) s , (cid:101) c = (cid:101) s , (cid:101) d = − (cid:101) V h (cid:101) s Therefore, the system of equations built by using (15) leads to: e aV + b + e cV + d = e (cid:101) aV + (cid:101) b + e (cid:101) cV + (cid:101) d (17a) ae aV + b + ce cV + d = (cid:101) ae (cid:101) aV + (cid:101) b + (cid:101) ce (cid:101) cV + (cid:101) d (17b) a e aV + b + c e cV + d = (cid:101) a e (cid:101) aV + (cid:101) b + (cid:101) c e (cid:101) cV + (cid:101) d (17c) a e aV + b + c e cV + d = (cid:101) a e (cid:101) aV + (cid:101) b + (cid:101) c e (cid:101) cV + (cid:101) d (17d)with the following two solutions obtained with Mathematica software[38]: a = (cid:101) c , b = (cid:101) d , c = (cid:101) a , d = (cid:101) ba = (cid:101) a , b = (cid:101) b , c = (cid:101) c , d = (cid:101) d Therefore, the system is s.o.l.i in the neighbourhood of V ∗ where theTaylor approximation holds. Corollary 1.
If the model (16) is completed with the following condi-tion: s > s it is trivial to see that the system (17)now has only one solution and themodel is s.o.g.i. for any V ∈ U nV ∗ ..