We are searching data for your request:

**Forums and discussions:**

**Manuals and reference books:**

**Data from registers:**

**Wait the end of the search in all databases.**

Upon completion, a link will appear to access the found materials.

Upon completion, a link will appear to access the found materials.

Let's take a channel that does not get inactivated, but is only in open and closed states. I understand what m (activation variable) represents. It's the fraction of the gates open.

But what does $ au$ (tau) , the time constant represent? Does it mean the time taken to open the gate? If tau is higher then it mean more time taken to open the gate?

I'm slightly confused.

Although the membrane time constant is very important in neuroscience, the time constant described in the original question is about *ion channel gating* rather than the rate of change of membrane voltage.

The "tau" when referring to ion channel gating is indeed referring to how long it takes to open/close.

When the membrane voltage changes, there is a new equilibrium value for the probability of the channel to be open versus closed (note: there can be other states for channels too, but for this explanation I will stick to open/closed). The *rate at which this equilibrium is approached* from the previous open probability is described by an exponential function, with "tau" as the time constant of that approach. A small value of tau means that the new equilibrium is reached very quickly. This site has a good explanation of the terminology and gating process, as do many others.

Time constants are involved in any system described by an exponential (these are called "linear time invariant" processes), whether decay or increase, and it is equal to the time it takes to progress approximately 63.2% of the way to the equilibrium. The cool thing about time constants (and linear time invariant processes in general) is that it doesn't matter what your starting point is or how far into the process you are, the time constant *stays constant* and every time "tau" passes you are now 63.2% closer to equilibrium than you were at time t-tau.

According to Wikipedia

In physics and engineering, the

time constant, usually denoted by the Greek letter τ (tau),is the parameter characterizing the response to a step input of a first-order, linear time-invariant (LTI) system.

**In Neuroscience:**

The time constant, $ au$, characterizes how rapidly current flow changes the membrane potential. $ au$ is calculated as:

$$ au = r_mc_m$$

where r$_m$ and c$_m$ are the resistance and capacitance, respectively, of the plasma membrane.

- The
**resistance across the membrane is a function of the number of open ion channels**and the capacitance is a function of the properties of the lipid bilayer

Note: the time constant is a "passive property" of membranes because they are intrinsic properties of all biological membranes. [Source].

My guess is that you are referring to some equation reference from Hodgkin-Huxley Model for the Generation of Action Potentials, but without further information about your specific equation, I cannot comment about Tau's specific role in your equation.

**Example Application:**

This constant *can* be applied to one of a few equations describing changes in membrane voltage. For example, a fall in membrane voltage is modeled as:

$$V(t) = V_{max}e^{-t/ au}$$

where voltage is in millivolts, time is in seconds, and $ au$ is in seconds.

- You can read more about the mathematics behind this and passive membrane properties in general here and here.

Again, From Wikipedia:

The larger a time constant is, the slower the rise or fall of the potential of a neuron. A long time constant can result in temporal summation, or the algebraic summation of repeated potentials. A short time constant rather produces a coincidence detector through spatial summation.

I Understand that you are referring to the 2-state model of ion channels. In this case, the transition from one steady-state to another is modeled as an exponential decay. $$ au * dm/dt = (m_inf - m) $$ The so called "fast" channels have $ au$ of a few milliseconds, "slow" ion channels have time constants of hundreds of milliseconds or more.

Of course, the 2-state model is a simplified mathematical model of ion channels, but for many phenomenological purposes, it is good enough.

## Stochastic models of ion channel gating

**Stochastic models of ion channel gating** date back as far as the pioneering work of Hodgkin and Huxley (1952), whose gating variables are often interpreted as probabilities. Sakmann and Neher directed the efforts that led to the first single-channel recordings, for which they won the Nobel Prize in Physiology and Medicine in 1991. Single-channel recordings (Sakmann and Neher, 1995) demonstrate that there are typically only two conductance states of individual ion channels: the closed and open states. Fluctuations between closed and open states appear to occur randomly ( Figure 1, top trace), presumably driven by thermal noise. Under voltage clamp, the *coefficient of variation* (ratio of standard deviation to mean) of the current passing through a set of independent and identically distributed stochastic channels is proportional to (1/sqrt

## Theoretical Modeling

The complex structure of the multidimensional conformational space of proteins implies an intricate kinetics despite an apparently simple bistability that is observed (22). Two popular theoretical approaches have been developed to cope with this complexity. A first one uses a simple bistable dynamics as a basis. To model the complexity of the observed kinetics this dynamics is amended by using an additional stochastic time dependence of the energy profile, or kinetic constants. Such an approach is nowadays commonly known under the label of “fluctuating barriers” (23–27). Alternatively, one can attempt to model the complexity of the energy profile itself in the simplest possible way. Our strategy is to find such a minimal model of the second kind that does allow for a rigorous analysis and does reproduce some nontrivial features of the gating dynamics.

Let us assume that the conformational stochastic dynamics between the open and closed states can be described in terms of a one-dimensional reaction coordinate dynamics *x*(*t*) in a conformational potential *U*(*x*) (Figs. 2 and 3). Because the distribution of open residence time intervals assumes typically a single exponential (1), in the following we rather shall focus on the behavior of the closed residence time intervals. To evaluate the distribution of closed residence time intervals it suffices to restrict our analysis to the subspace of closed states by putting an absorbing boundary at the interface, *x* = *b*, between the closed and open conformations (see Fig. 3). We next assume that the gating dynamics is governed by two gates: an inactivation gate and an activation gate. The inactivation gate corresponds to the manifold of *voltage-independent* closed substates. It is associated with the flat part, −*L* < *x* < 0, of the potential *U*(*x*) in Fig. 3. In this respect, our modeling resembles that in ref. 28. The mechanism of inactivation in potassium channels is quite sophisticated and presently not totally established (1). It is well known that inactivation can occur on quite different time scales (1). The role of a fast inactivation gate in *Shaker* K + channels is taken over by the channel's extended N terminus, which is capable of plugging the channel's pore from the cytosol part while diffusing towards the pore center (29). The slow inactivation apparently is due to a conformational narrowing of the channel pore in the region of selectivity filter (1). In both cases, no net gating charge translocation occurs and the inactivation process does not depend on voltage. When the inactivating plug is outside of the pore, or the selectivity filter is open (*x* > 0 in Fig. 3) the channel can open only *if* the activation barrier is overcome.

Studied model and its diffusion counterpart.

The dynamics of the activation gate occurs on the linear part of the ramp of the potential *U*(*x*)—i.e., on 0 < *x* < *b* in Fig. 3, as in refs. 18 and 19. Note that for 0 < *x* < *b*, the inactivating plug diffuses outside of the channel's pore and the selectivity filter is open. During the activation step a gating charge *q* moves across the membrane this feature renders the overall gating dynamics voltage dependent. The channel opens when the reaction coordinate reaches the location *x* = *b* in Fig. 3. This fact is accounted for by putting an absorbing boundary condition at *x* = *b*. Moreover, the channel closes immediately when the inactivation gate closes (*x* ≤ 0), or when the activation gate closes. To account for this behavior in extracting the closed residence time distribution we assume that the channel is reset into the state *x* = 0 after each closure (see below).

The diffusional motion of the inactivated gate is restricted in conformational space. We characterize this fact by the introduction of a conformational diffusion length *L* (Fig. 3) and the diffusion constant *D* ∼ *k*_{B}*T* that are combined into a single parameter—the conformational diffusion time 2 This quantity constitutes an essential parameter for the theory. We assume that the activation barrier height *U*_{0} is linearly proportional to the voltage bias *V* (18, 19)—i.e., in terms of the gating charge *q* we have 3 Moreover, *U*_{0} is positive for negative voltages—i.e., for *V* < *V*_{c}—vanishes at *V* = *V*_{c}, and becomes negative for *V* > *V*_{c}. Thus, for *V* > *V*_{c} the channel “slips” in its open state, rather than overcomes a barrier. In addition, the fraction ξ of the voltage-dependent substates in the whole manifold of the closed states should be very small, implying that 4

### Analytical Solution.

The corresponding Fokker–Planck equation for the probability density of closed states *P*(*x*, *t*) reads 5 where β = 1/(*k*_{B}*T*). To find the distribution of closed residence times *f*_{c}(*t*), we solve Eq. 5 with the initial condition *P*(*x*,0) = δ(*x*), in combination with a reflecting boundary condition *dP*(*x*, *t*)/*dx* |_{x=−L} = 0, and an absorbing boundary condition, *P*(*x*, *t*)|_{x=b} = 0 (4). The closed residence time distribution then follows as 6 where Φ_{c}(*t*) = ∫ *P*(*x*,*t*)*dx* is the survival probability of the closed state.

By use of the standard Laplace transform method we arrive at the following *exact* solution: 7 where 8 9 The explicit result in **7–9** allows one to find all moments of the closed residence time distribution. In particular, the mean closed residence time 〈*T*_{c}〉 = lim_{s→0}[1 − f̃_{c}(*s*)]/*s* reads 10 This very same result **10** can be obtained alternatively if we invoke the well-known relation for the mean first-passage time 〈*T*_{c}〉 = 1/*D* ∫ *dxe* β*U*(*x*) ∫ *dye* −β*U*(*y*) (4). This alternative scheme provides a successful validity check for our analytical solution in **7–9**.

### Elucidation of the Voltage Dependence in Eq. 1.

Upon observing the condition **4**, Eq. 10 by use of **3** reads in leading order of ξ 11 With the parameter identifications 12 and 13 the result in **11** precisely coincides with Eq. 1. The fact that our approach yields the puzzling voltage dependence in Eq. 1 constitutes a first important result of this work.

Let us next estimate the model parameters for a *Shaker IR* K + channel from ref. 6. In ref. 6, the voltage dependence of *k*_{o}(*V*) at *T* = 18°C has been parameterized by Eq. 1 with the parameters given in the caption of Fig. 1. Then, from Eq. 12 the gating charge can be estimated as *q* ≈ 20*e* (*e* is the positive-valued elementary charge). As to the diffusion time τ_{D}, we speculate that it corresponds to the time scale of inactivation the latter is in the range of seconds and larger (6). Therefore, we use τ_{D} = 1 sec as a lower bound for our estimate. The fraction of voltage-dependent states ξ is then extracted from Eq. 13 to yield, ξ ≈ 0.0267. This value, indeed, is rather small and thus proves our finding in Eq. 11 to be consistent.

### Analysis for the Closed Residence Time Distribution.

The exact results in Eqs. **7–9** appear rather entangled. To extract the behavior in real time one needs to invert the Laplace transform numerically. With ξ ≪ 1, however, Eqs. **7–9** are formally reduced to 14 This prominent leading order result can be inverted *analytically* in terms of an infinite sum of exponentials, yielding: 15 where the rate constants 0 < λ_{1} < λ_{2} < … are solutions of the transcendental equation 16 and the expansion coefficients *c*_{n}, respectively, are given by 17 Note from Eq. 6 that the set *c*_{n} is normalized to unity, i.e. Σ *c*_{n} = 1.

The analytical approximation, Eqs. **15–17**, is compared in Fig. 4 with the precise numerical inversion of the exact Laplace transform in Eqs. **7–9**. The numerical inversion has been performed with the Stehfest algorithm (30). As can be deduced from Fig. 4, for *t* > 10 msec the agreement is very good indeed. A serious discrepancy occurs only in the range 0.01 msec < *t* < 0.1 msec, which lies outside the range of the patch clamp experiments (*t* ≳ 0.1 msec). Moreover, the agreement improves with increasing τ_{D} (not shown).

Closed residence time distribution for a diffusion-limited case. The numerical precise result (solid line) is compared with the analytical approximation in Eqs. **15–17** (broken line). The latter one coincides with the exact solution of the diffusion model by Millhauser *et al.* (12) in the scaling limit.

### Origin of the Power Law Distribution.

The features displayed by the closed residence time distribution *f*_{c}(*t*) depend sensitively on the applied voltage *V*. When *V* > *V*_{c}—e.g., *V* = −45 mV, as in Fig. 4—the activation barrier towards the channel opening disappears and the opening dynamics becomes diffusion limited. In this case, the diffusion time τ_{D} = 1 sec largely exceeds the mean closed residence time 〈*T*_{c}〉 ≈ 18.4 msec. Put differently, τ_{D} ≫ 〈*T*_{c}〉 and the closed residence time distribution exhibits an intricate behavior with three distinct regions (see Fig. 4). Most importantly, for the intermediate time scale 18 we find from Eq. 14 (by considering the limit τ_{D} → ∞) that the closed residence time distribution obeys a power law, reading 19 This type of behavior is clearly detectable in Fig. 4, where it covers about two decades of time. As follows from Eq. 18, an increase of τ_{D} by *one* order of magnitude (while keeping 〈*T*_{c}〉 fixed) extends the power law region by *two* orders of magnitude. This conclusion is fully confirmed by our numerics (not shown). This power law dependence, which extends over four orders of magnitude, has been seen experimentally for a K + channel in NG 108-15 cells (11). On the contrary, for channels, where τ_{D} is smaller, the power law region **18** shrinks and eventually disappears, whereas the mean opening rate defined by Eq. 10 still exhibits a steep dependence on the voltage. Thus, our model is capable of describing for different channels both the emergence of the power law and its absence.

On the time scale *t* ≥ τ_{D} the discussed power law distribution crosses over into the exponential tail the latter is fully described by the smallest exponent λ_{1} in Eq. 15, i.e., by 20 This feature is clearly manifest in Fig. 4. The transition towards the exponential tail in the closed residence time-interval distribution can be used to estimate the diffusion time τ_{D} on pure *experimental* grounds!

Finally, let us consider the opposite limit, τ_{D} ≪ 〈*T*_{c}〉, *for* *V* ≪ *V*_{c}. For the considered set of parameters this occurs—e.g., for *V* = −55 mV—when the channel is predominantly closed. Then, the diffusion step in the opening becomes negligible and in the experimentally relevant range of closed residence times, defined by 〈*T*_{c}〉, the corresponding distribution can be approximated by a single exponential, **20**. A perturbation theory in Eq. 16 yields λ_{1} ≈ *k*_{o} (1 − (*k*_{o}τ_{D})/3). For the parameters used we have λ_{1} ≈ 0.96*k*_{o} and, from Eq. 17, *c*_{1} ≈ 0.95. This is in a perfect agreement with the precise numerical results obtained from Eqs. **7–9**. Thus, the distribution of closed residence times is single exponential to a very good degree. Consequently, one and the same channel can exhibit both an exponential and a power-law distribution of closed residence times, as a function of the applied transmembrane voltage. With an increase of τ_{D} the voltage range of the exponential behavior shifts towards more negative voltages, *V* < *V*_{c}, and vice versa.

### Reduction to a Diffusion Model.

Let us relate our model to that introduced previously by Millhauser *et al.* (12). The latter one is depicted with the lower part in Fig. 3. It assumes a discrete number *N* of closed substates with the same energy. The gating particle jumps with the equal forward and backward rates *k* between the adjacent states, which are occupied with probabilities *p*_{n}(*t*). At the right edge of the chain of closed states the ion channel undergoes transition into the open state with the voltage-dependent rate constant γ. To calculate the closed residence time distribution *f*_{c}(*t*) one assumes *p*_{0}(0) = 1, *p*_{n≠0}(0) = 0, and *d*Φ_{c}(*t*)/*dt* = −γ*p*_{0}(*t*), where Φ_{c}(*t*) = Σ *p*_{n}(*t*) is the survival probability (12, 17).

We consider the continuous diffusion variant of this model (31) in a scaling limit: we put Δ*x* → 0, *k* → ∞, γ → ∞, *N* → ∞, keeping the diffusion length *L* = *N* Δ*x*, the diffusion constant *D* = *k*(Δ*x*) 2 , and the constant *k*_{o} = γ/*N* all finite. The latter one has the meaning of mean opening rate (see below). Note that in clear contrast with our approach, the rate parameter *k*_{o} in the diffusion model is of pure phenomenological origin. The problem of finding the closed residence time distribution is reduced to solving the diffusion equation 21 with the initial condition *P*(*x*,0) = δ(*x* − 0_{−}), the reflecting boundary condition ∂*P*(*x*, *t*)/∂*x* |_{x=−L} = 0 and the radiation boundary condition (32). 22 We emphasize that the radiation boundary condition **22** is not postulated, but is rather *derived* from the original discrete model in the considered scaling limit. Using the Laplace transform method, we solved this problem exactly and obtained the result in Eq. 14. In conclusion, our approximate result in Eqs. **14–17** provides the *exact* solution of the diffusion model (12, 17) in the scaling limit! Note, however, that this diffusion model so obtained is not able to resolve the puzzling voltage dependence in Eq. 1.

## Materials and Methods

Methods have previously been described 18 and are presented in detail in the online supplement. The mAno1 **a, c** splice variant (Accession: <"type":"entrez-protein","attrs":<"text":"Q8BHY3","term_id":"148887069","term_text":"Q8BHY3">> Q8BHY3) was used. Mutations were made using PCR-based mutagenesis. mAno1 was transfected into HEK293 cells using Fugene-6 (Roche Molecular Biochemicals, Indianapolis, IN). Transfected HEK293 cells were patch clamped using conventional whole-cell and excised inside-out patches. The zero Ca 2+ intracellular solution contained (mmol/L): 146 CsCl, 2 MgCl_{2}, 5 EGTA, 10 HEPES, 10 sucrose, pH 7.3, adjusted with NMDG. High Ca 2+ pipette solution contained 5 mmol/L Ca 2+ -EGTA, instead of EGTA (free Ca 2+ 20 μM). The 126 μM and 2 mmol/L Ca 2+ were made by adding 0.2 mmol/L and 2 mmol/L CaCl_{2} to high-Ca 2+ solution. The standard extracellular solution contained (mmol/L): 140 NaCl, 5 KCl, 2 CaCl_{2}, 1 MgCl_{2}, 15 glucose, 10 HEPES, pH 7.4 with NaOH. Relative anion permeability was determined by measuring the shift in zero-current *E*_{rev} after changing the bath solution from 151 mmol/L Cl − to 140 mmol/L substitute anion (X − ) plus 11 mmol/L Cl − . The permeability (P) ratio was calculated by using the Goldman-Hodgkin-Katz equation. Permeability of Na + or Cs + relative to Cl − was determined by measuring changes in zero-current *E _{rev}* when the concentration of extracellular NaCl or CsCl was changed (𠇍ilution potential” method). The fast application of Ca 2+ to excised inside-out patches was performed using a double-barreled theta tubing with a tip diameter of

50 μm attached to a piezobimorph on a micromanipulator. One barrel was filled with standard zero-[Ca 2+ ]_{i} solution, and the other barrel was filled with intracellular solution containing the indicated free Ca 2+ . The time course of solution exchange across the laminar flow interface was estimated by liquid junction potential measurements to be 0.5 ms (10�% rise time) for a 10-fold difference in ionic strength. mAno1 with inserted human influenza hemaglutinin (HA) epitopes at various positions were used to determine the topology by assessing the accessibility of the HA epitope to extracellularly applied antibody. Cells were fixed for 15 min at room temperature in 1% paraformaldehyde in 0.1mol/L phosphate buffer pH 7. Non-permeabilized cells were washed in 3 times in blocking buffer. Permeabilized cells were incubated in blocking buffer containing 0.15% - 0.3% Triton-X100. Cells were incubated with anti-HA antibody diluted 1:750 in blocking buffer for 2 h at room temperature, washed 3 times in blocking buffer, and then incubated in goat-anti-rabbit IgG conjugated with Dylight-549 (Jackson Immunochemicals).

## Abstract

Mathematical models of cardiac action potentials have become increasingly important in the study of heart disease and pharmacology, but concerns linger over their robustness during long periods of simulation, in particular due to issues such as model drift and non-unique steady states. Previous studies have linked these to violation of conservation laws, but only explored those issues with respect to charge conservation in specific models. Here, we propose a general and systematic method of identifying conservation laws hidden in models of cardiac electrophysiology by using bond graphs, and develop a bond graph model of the cardiac action potential to study long-term behaviour. Bond graphs provide an explicit energy-based framework for modelling physical systems, which makes them well suited for examining conservation within electrophysiological models. We find that the charge conservation laws derived in previous studies are examples of the more general concept of a ‘conserved moiety’. Conserved moieties explain model drift and non-unique steady states, generalizing the results from previous studies. The bond graph approach provides a rigorous method to check for drift and non-unique steady states in a wide range of cardiac action potential models, and can be extended to examine behaviours of other excitable systems.

### 1. Introduction

Models of the cardiac action potential have been developed to study cardiac diseases, such as arrhythmia [1–3], ischaemia [4] and acidosis [5]. Increasing model complexity has led to concerns over the occurrence of drift and non-unique steady states [6–8], particularly for extensions of the DiFrancesco & Noble [9] and Luo–Rudy [1,3] models. While solutions to these issues have been proposed using conservation principles [10,11], they have not been universally applied for more recent models, many of which still use non-conservative stimulus currents that predispose them to drift [12–14]. More recently, the Food and Drug Administration (FDA) has initiated plans to use cardiac action potential models to assess potential drug side-effects on cardiac instability through the human *ether-à-go-go*-related gene (hERG) K + channel. Thus, with an increasing emphasis on model robustness and accuracy, there is a renewed incentive to resolve the issues of drift and non-unique steady states [15,16].

Drift is the failure of a model to reach a consistent limit cycle when simulated over long periods, and is often caused by a non-conservative stimulus containing current with no charge carrier [6,10]. Hund *et al.* [10] derived a charge conservation law, and found that non-conservative stimulus currents violate this conservation law, hence they proposed K + ions as the current charge carrier to resolve this. A related issue in many models where drift has been resolved is that steady-state limit cycles under constant pacing depend upon the initial conditions and are therefore non-unique [7,8,10]. Thus, depending on the initial conditions, the same model may lead to different conclusions. Like drift, authors have suggested that charge conservation can constrain initial conditions such that they lead to the same steady state [8,10,11].

While the studies by Hund *et al.* [10] and Livshitz & Rudy [11] suggest measures to eliminate drift and attain a unique steady state by using conservation laws, their analyses are limited in their scope and not a comprehensive solution for all models. Because existing studies [10,11] explore charge conservation only in specific models, and the conservation laws were derived from physical intuition rather than a principled mathematical approach, it is difficult to generalize their findings to other models where charge conservation is routinely neglected. Furthermore, because these studies focus only on conservation of charge, they may miss other conservation laws relevant for long-term behaviour, such as those corresponding to ions, ion channels and buffers. A general approach is, therefore, desirable to deal with the issues of drift and steady states in a more systematic manner and for a broader range of models.

To facilitate a general approach, we propose the use of bond graphs which explicitly model energy transfer across physical systems to ensure compliance with conservation principles. Bond graphs were initially invented to model hydroelectric systems [17] and they have subsequently been extended to model chemical [18], biochemical [19,20] and electrochemical systems [21]. As with all physical systems, biological processes must obey the fundamental principles of physics and thermodynamics [22] therefore, bond graphs are well suited for constraining models of biological systems to physically plausible solutions [23], and also for inferring the energetic cost of biological processes [21,23–25]. Because the bond graph representation emphasizes analogies between different physical domains, electrophysiological systems can be analysed as an analogous biochemical system with a stoichiometric matrix that describes the stoichiometry of each reaction within its columns [20,26–28]. In this context, the ‘conservation principle’ described in earlier studies is an example of the more general principle of a conserved moiety in metabolic and bond graph analysis [23,29].

In this study, we develop an abridged bond graph model of the cardiac action potential and outline a general approach to study the effects of conserved moieties on drift and steady-state behaviour. Our bond graph model simulates the essential features of the cardiac action potential, and because bond graphs are energy-based this easily provides an estimate of the energetic cost (in Joules) of the cardiac action potential. Our analysis reveals conservation of charge as one of the conserved moieties of our model, along with other conserved moieties corresponding to ions, channels, transporters and buffers. We observed that our model solution was subject to drift when the stimulus current violated any conservation laws corresponding to the conserved moieties, and that changes to the initial conditions led to different steady states if the value of any conserved moiety was changed. To demonstrate that our approach is general, we analyse variants of our bond graph model where different ions have been fixed at a constant concentration (corresponding to ‘chemostats’). It should be noted that fixing an ion concentration can change the conserved moieties of a system, therefore influencing a model's susceptibility to drift and non-unique steady states. The bond graph approach is a useful and general method to identify and interpret conservation principles, and it can link conserved moieties to individual steady states. Our approach can be used to automatically derive charge conservation laws that are frequently neglected in existing models of cardiac electrophysiology, and we build upon existing reports [10,11] to propose solutions for drift and non-unique steady states which work for all cardiac action potential models that can be represented using bond graphs.

### 2. Methods

#### (a) Model components

To study the issues of drift and non-unique steady states, we built an abridged bond graph model of the cardiac action potential, with the minimal number of channels and pumps required to simulate the essential features of a cardiac action potential, and maintain ionic concentrations over long periods of simulation. We note that a series of mathematical equations can only be described using bond graphs if they describe a thermodynamically consistent system, and that no existing model of cardiac electrophysiology is entirely thermodynamically consistent. As seen in Gawthrop *et al.* [21], an exact translation of existing ion channel models into the bond graph framework is generally not possible. Therefore, for some transport processes and particularly for ion channels, it was necessary to build new models under the bond graph framework, although their parameters and structure can be chosen based on equations in existing models. Accordingly, many of the essential components of our model were based primarily on the Luo–Rudy 1994 dynamic model [1], although it is possible to use other models and/or model more sub-cellular processes. We argue in §2c that some aspects of the equations for the ion channels described in the Luo and Rudy model are not thermodynamically consistent, and unable to be represented as bond graphs. We describe our approach to modelling ion channels in §2d. As we use simple and elementary bond graph structures based on physical principles, our bond graph representation of ion channels is a physically constrained approximation of the equations described in Luo and Rudy, although it is still able to simulate the essential features of the cardiac action potential.

Model components are shown in figure 1*a*, together with the overall bond graph structure (figure 1*b*). Ion channels and Ca 2+ buffering components were based upon their representations in Luo & Rudy [1]. The L-type Ca + channel in the Luo–Rudy model is permeable to calcium, sodium and potassium, but we neglected its sodium conductance as this has a relatively small contribution to the action potential. The Na + /K + ATPase model was based on the model by Terkildsen *et al.* [4], with modifications suggested by Pan *et al.* [30] to allow conversion into a bond graph model. The equation for the Na + –Ca 2+ exchanger (NCX) current in Luo and Rudy did not have an obvious correspondence to a bond graph structure, thus we modelled this component using a simple bond graph module that was fitted to experimental data [31,32]. The general approach for modelling ion channels is described later in the methods, with details on the other components, as well as further detail on ion channel modelling given in the electronic supplementary material.

Figure 1. Action potential model. (*a*) Cell schematic, (*b*) overall bond graph structure. The bond graph modules Na_channel, NaK, K_channels, LCC, NCX and Ca_buffer contain more detailed aspects of the bond graph structure which are described further in the Methods and electronic supplementary material. Coloured bonds link bond graph modules to the appropriate chemical species. Definitions: *I*_{Na}, sodium current *I*_{K1}, time-independent K + current *I*_{K}, time-dependent K + current *I*_{Kp}, plateau K + current *I*_{LCC}, L-type Ca 2+ current NCX, Na + –Ca 2+ exchanger Na/K, Na + /K + ATPase TRPN, troponin CMDN, calmodulin. (Online version in colour.)

#### (b) Bond graph modelling

Here, we briefly outline bond graph components as used in electrophysiological modelling. For a more comprehensive introduction, the texts by Gawthrop & Smith [33] and Borutzky [34] provide detailed descriptions of bond graph theory, and Gawthrop & Bevan [35] provide a short tutorial for engineers. Theory for bond graph modelling of biochemical systems can be found in [19,20,23,36].

Bond graphs consist of components (representing physical objects and processes), bonds (representing the transfer of energy) and junctions (representing network structure). Each bond carries two variables: an effort *e* and a flow *f* , such that their product determines the power of the bond (i.e. *p* = *ef*). Thus, bond graphs explicitly account for energy transfer, and are thermodynamically consistent. Because effort and flow are generalized variables, they can represent quantities from a variety of physical systems, including mechanical (*e* = force (N), *f* = velocity (m s −1 )), electrical (*e* = voltage (V), *f* = current (A)) and hydraulic systems (*e* = pressure (Pa), *f* = volumetric flow rate (m 3 s −1 )) [34].

The network structure of a bond graph is specified by 0 and 1 junctions. The 0 (or effort) junctions specify that efforts of all connected bonds are equal, and thus to ensure conservation of energy through this junction, the flows of the bonds must sum to zero. In the electrical and hydraulic domains, 0 junctions represent parallel connections, whereas they represent series connections in the mechanical domain. By a similar principle, 1 (or flow) junctions specify that the flows of all connected bonds are equal, ensuring that their efforts sum to zero. Thus, 1 junctions correspond to series connections in the electrical and hydraulic domains, and parallel connections in the mechanical domain.

To illustrate the use of a bond graph for electric circuit analysis, we consider the electric circuit where two capacitors are connected to a resistor in series (figure 2*a*). All components are linear, described by the following equations:

Figure 2. Conceptual representations of physical systems. (*a*) A bond graph for the illustrated simple electric circuit with two capacitors and a resistor in series. (*b*) A bond graph analogous to the electric circuit in (*a*) can also represent the chemical reaction A ⇌ B . (*c*) Bond graphs can also model interaction of components in both the chemical and physical domains, such as the transport of an ion across a membrane. (*d*) Transport of an ion across a membrane through an ion channel involves gating which modulates the rate of reaction. Thus the ion channel is analogous to a potentiometer. (Online version in color.)

More recently, bond graphs have been extended to model biochemical systems [19,20] where the chemical potential *μ* (J mol −1 ) is the effort variable, and molar flow rate *v* (mol s −1 ) is the flow variable. As temperature and pressure are assumed to be constant in biochemical systems, the measure of thermodynamic potential *μ* corresponds to Gibbs free energy [37]. For bond graph modelling of more general systems where temperature and pressure change, the reader is directed to the text by Thoma & Bouamama [38].

Because mass action equations are nonlinear in general, nonlinear components are required to model biochemical systems. Each chemical species is represented as a capacitor. However, in contrast to the electrical domain, the constitutive equation for the capacitor representing each species is logarithmic:

The bond graph framework for biochemistry can be extended to electrochemical systems [21] as demonstrated in figure 2*c*, which models the transport of a positively charged species *X* across a membrane. It should be noted that chemical species are described with C components that have a logarithmic association, whereas the C component corresponding to the (electric) membrane potential has a linear constitutive relationship. A transformer (TF) is used to convert the membrane voltage into an equivalent chemical potential through Faraday's constant *F* = 96485 C mol −1 , such that:

We chose to represent ion channels such that conductance was modulated by membrane voltage, both directly and indirectly through gating processes. A bond graph representation for this relationship is given in figure 2*d*. As shown, this model has the same electrical representation as figure 2*c* however it uses a variable resistor. The bond graph representation contains the same states, with C : xi, C : xe, and C : mem (with a transformer) connected through 0 junctions. In this case however, the Re components that describe the constitutive relation have been changed, such that Re_GHK : r1 is connected to an additional effort that modulates its velocity, and the gating affinity *A* g is added to both the forward and reverse affinities to describe changes in permeability due to gating. Further detail on modelling ion channels using bond graphs is given in §2d.

#### (c) Modelling approach

Because bond graphs constrain the equations of a model to ensure thermodynamic consistency, they only allow the representation of the subset of general mathematical models that are thermodynamically consistent. In such models, dissipative processes such as reactions and ionic currents can only proceed in the direction of decreasing chemical or electrochemical potential, and when there is no potential gradient (i.e. the process is at equilibrium), the process must stop. However, many existing models of cardiac electrophysiology are not constrained to ensure this behaviour, and in some cases, they describe physically infeasible systems that create energy out of nowhere [20]. Therefore, many existing models do not have a direct bond graph representation [23]. An advantage of the bond graph approach is that the discipline required to convert an existing model into a bond graph helps to highlight thermodynamic issues and inconsistencies in existing models that would have otherwise been missed or ignored (see Gawthrop *et al.* [23] for an example of this process applied to a model of glycolysis). For the currents in this study, the Luo–Rudy equations for the time-dependent K + and L-type Ca 2+ channels are not thermodynamically consistent, and therefore cannot be directly converted into a bond graph model. As shown in equation (2.11), a consequence of the thermodynamic consistency of the bond graph framework is that there is a constraint on the equilibrium point for each ion channel that is determined by the Nernst equation. As the equilibrium points of the time-dependent K + and L-type Ca 2+ channels do not correspond to the Nernst potential (see electronic supplementary material, A.1.3 and A.1.5), they are not thermodynamically consistent and therefore their equations cannot be translated into a dissipative bond graph component. These issues were addressed by finding a thermodynamically consistent approximation that satisfies the Nernst equation. A second issue is that it is difficult to simultaneously model open-channel currents and channel gating using the elementary components, described in §2b and figure 2, thus we approximated these relations using the components described in that section rather than using more complicated components. Therefore, rather than attempting to reproduce the Luo–Rudy equations exactly, we built a bond graph structure as implied by the equations in the Luo and Rudy model, and chose parameters of our bond graph model to fit aspects of the Luo–Rudy model as closely as possible, specifically the current-voltage (I-V) curves and gating parameters. For all other components conversion into a bond graph model was more straightforward, and we used the methods of Gawthrop *et al.* [23]. Further information on the bond graph model, and parameter identification is given in the electronic supplementary material.

#### (d) Ion channel modelling

##### (i) Bond graph structure.

In this section, we discuss decisions made in developing models of ion channels. The bond graph structure for the Kp channel is shown in figure 3. The other channels have similar structures that follow from the discussion in this section.

Figure 3. The bond graph model of a plateau K + channel. (*a*) The **channel_Kp** module describes the current through the ion channel. (*b*) The **gate_en_Kp** contains the states required for gating. (*c*) The **vRe** module contains a voltage-dependent reaction used to describe channel state transitions. (*d*) The channel current and gating modules are combined into an ion channel model ( **Kp_channel** ). (Online version in colour.)

##### (ii) Current–voltage relations.

While thermodynamic properties can be used to determine how membrane voltage and ionic concentrations relate at equilibrium, they do not specify behaviour away from equilibrium. For this purpose, the current-voltage (I-V) relationship defines how the membrane voltage relates to the current through a specific channel. In many electrophysiology models [1], currents through ion channels are described through the use of linear I-V relationships based on empirical fits to data rather than fundamental physical principles. Using bond graphs, it is difficult to incorporate the effects of gating using a linear I-V equation. While linear I-V equations can be thermodynamically consistent provided they satisfy the Nernst equation in equation (2.11), and such equations can be represented by using a linear R component, the modulation of its conductance would require the use of signal bonds that do not necessarily represent physical processes [34]. Therefore, we use the Goldman–Hodgkin–Katz (GHK) equation to model ion channels, as it enables relatively simple incorporation of ion channel gating as a physics-based biochemical module [21]. The GHK equation defines a nonlinear relationship between current *I* and membrane voltage *V* :

As many ion channels in the Luo–Rudy model are described using a linear I-V relationship, the use of GHK equations requires approximation. A comparison of the resulting I-V curves is given in figure 4. Details of the fitting process, as well as the fitted bond graph parameters, are given in the electronic supplementary material. The Na + channel I-V curves appeared to match reasonably well (figure 4*a*), with some discrepancies at positive membrane potentials. For K + channels (figure 4*b*–*d*), we attempted to optimize the fit across voltages that correspond to their physiological function, so that their currents would be most similar to those of Luo and Rudy when the ion channels are open. Accordingly, for *I*_{K1} (−90 mV≤*V* < − 30 mV), *I*_{K} (−20 mV≤*V* ≤30 mV) and *I*_{Kp} (*V* > 0 mV) the I-V curves matched reasonably well in these regions. Discrepancies occurred outside these ranges of voltages, but appeared to only cause minor differences to the currents. In their implementation of *I*_{K}, Luo & Rudy [1] use a thermodynamically inconsistent I-V equation where the current is non-zero at the Nernst potential for K + . Despite this, bond graph parameters could still be chosen to give a reasonable fit to this I-V equation (figure 4*c*). Because the Luo–Rudy model based their L-type Ca 2+ I-V curves on the GHK equation, there was a far closer match between the bond graph and Luo–Rudy models for these currents, (figure 4*e*,*f* ) and the K + curve was matched exactly (figure 4*f* ).

Figure 4. Comparison of I-V curves between the Luo–Rudy (LRd) and bond graph (BG) models. (*a*) *I*_{Na}, (*b*) *I*_{K1}, (*c*) *I*_{K}, (*d*) *I*_{Kp}, (*e*) *I*_{Ca,L}, (*f*) *I*_{K,L}. The standard concentrations in Luo & Rudy [1] ([Na + _{i}] = 10 mM, [Na + _{e}] = 140 mM, [ K + _{i} ] = 145 mM , [K + _{e}] = 5.4 mM, [Ca + _{i}] = 0.12 μM, [Ca + _{e}] = 1.8 mM) were used to match I-V curves.

##### (iii) Modulation.

While the I-V curves describe currents through open ion channels, a formulation for gating is required to describe the number of open ion channels at any given time. In the Hodgkin–Huxley framework, gating is modelled as differential equations that give the proportion of open gates at any given time. We incorporated the effects of gating through a gating affinity *A* g , which is added to both the forward and reverse affinities of a reaction (figure 3*a*) to modulate its rate without changing the equilibrium [21].

##### (iv) State models.

Ion channel models must account for gating and bond graphs require the use of physical components to achieve this. We model gating as transitions between channel states, known in the literature as Markov models [41,42]. To illustrate, we use the example of a typical Na + channel in which the current *I* is described by the equation

*S*

_{31}represents the open channel. Because individual channel states are modelled, the current depends only on the amount of

*S*

_{31}and not any of the other closed states. Thus, incorporation into the gating framework described above is intuitive each state represents a structural conformation of the ion channel and the number of channels in each state is explicitly tracked, facilitating a simple approach to account for the energetics of gating under varying ion channel densities.

Figure 5. Channel states of a Na + channel.

##### (v) Voltage dependence of state transitions.

The transition rates between open and closed states are voltage dependent for ion channels. Hodgkin–Huxley models describe state transitions using ODEs of the form

To assess the quality of fit, we compare steady-state open probabilities *g*_{ss} = *α*(*V* )/(*α*(*V* ) + *β*(*V* )) and time constants *τ* = 1/(*α*(*V* ) + *β*(*V* )) (figure 6). The curves for *g*_{ss} and *τ* were generally in agreement however, there were some exceptions. In particular, time constants for the Na + channel gates have lower peaks in the bond graph model when compared to the Luo–Rudy model (figure 6*a*–*c*), but this did not appear to significantly affect Na + channel function as the peaks were all decreased by a similar proportion, facilitating coordination between opening and closing. Similarly, the time constant *τ*_{d} (figure 6*h*) was lower in the bond graph model for some voltages, but given that discrepancies occur at time constants much smaller than the time course of a cardiac action potential we expect that the effects would be negligible. Finally, for the time-dependent K + current *X*_{ss} is substantially higher at negative voltages so that the bond graph model can provide a better match at positive voltages (figure 6*e*). The effects of this difference are partially offset by the lower GHK current at negative voltages which are still above the Nernst potential of K + (figure 4*c*).

Figure 6. Fits for bond graph (BG) parameters against corresponding gating equations from the Luo–Rudy (LRd) model. Steady-state open probabilities are shown on the left panels, and time constants are shown on the right. The bond graph equations are plotted with solid lines, and the Luo and Rudy equations in dashed lines. Gates include (*a*) *m*, sodium activation (*b*) *h*, sodium inactivation (*c*) *j*, slow sodium inactivation (*d*) K1, time-independent K + activation (*e*) X, time-dependent K + activation (*f*) X i , time-dependent K + inactivation (*g*) Kp, plateau K + activation (*h*) *d*, L-type Ca 2+ channel activation. Note that the X i and Kp gates were originally formulated as steady-state equations, thus time constants are shown only for matched bond graph parameters.

#### (e) Finding conserved moieties

Within a biochemical model, conserved moieties are chemical structures that are neither created, removed nor broken down. A common example in energy-dependent metabolic networks is the adenosine moiety found in AMP, ADP and ATP [23,29]. Mass balance specifies that the total amount of each conserved moiety remains constant, and if information on the molecular structure of each species of a reaction network is available, these conservation laws can be derived by counting the number of moieties across all species [29]. In practice, many models do not contain this structural information and this approach cannot be used however, the conservation laws still hold. Here, we outline a method to find conserved moieties using stoichiometric information rather than chemical structures.

Models of cardiac electrophysiology can be represented by the differential equation

#### (f) Stimulus currents

The cardiac action potential model was stimulated using a constant current stimulus that contained enough charge to raise the membrane potential by 30 mV over 0.1 ms. The non-conservative stimulus currents consist of non-specific charge that is not carried by any specific ion. As recommended by Kneller *et al.* [8], conservative stimulus currents contained K + ions as the charge carrier.

### 3. Results

#### (a) Simulation of a single action potential

To verify that our bond graph model reproduced the features of a typical action potential, we simulated the model over a single beat (figure 7*a*–*c*). The membrane potential (figure 7*a*, with stimulation indicated by the arrow) resembled a typical cardiac action potential, with a distinct peak and plateau phase. The contributions of ion channel currents reproduce some common features of cardiac action potentials (figure 7*b*). Once the action potential is initiated by a stimulus current, the sodium current *I*_{Na} briefly activates to give rise to a voltage spike. Following this, the plateau phase occurs where depolarizing L-type Ca 2+ currents oppose the repolarizing K + currents *I*_{K} and *I*_{Kp}. Towards the end of the action potential, *I*_{K1} activates to restore the resting potential [50]. Our model also simulates the reversal of NCX current across the action potential, and the consistent outward current of the Na + /K + ATPase to maintain ionic gradients (figure 7*c*). As a consequence of incorporating the voltage-dependence of gating transitions in a physical framework, transitions between channel states are associated with a gating current resulting from the movement of such charged residues. Our model reveals that the total gating current across all channels *I*_{gate} has minimal contribution to total current (figure 7*c*).

Figure 7. A simulation of the cardiac action potential using a bond graph model. (*a*) Membrane voltage, following stimulation with a conservative stimulus current (arrow) (*b*) ion channel currents (*c*) transporter and gating currents (*d*) membrane voltage over three cycles, for comparison with (*e*) and (*f*) (*e*) power consumption (*f*) energy dissipated, with the variable *E* representing the energy consumption over the duration of the action potential. The model was run initially for 300 ms to allow the membrane potential and channel gates to stabilize. The intracellular ion concentrations were dynamic variables with initial concentrations [Na + _{i}] = 10 mM, [K + _{i} ] = 145 mM and [Ca + _{i}] = 0.12 μM. Constant concentrations were [Na + _{e}] = 140 mM, [K + _{e}] = 5.4 mM, [Ca + _{e} ] = 1.8 mM , [MgATP] = 6.95 mM, [MgADP] = 0.035 mM, [P _{i} ] = 0.3971 mM and pH = 7.095. *T* = 310 K . (Online version in colour.)

Figure 7*e* shows the power consumption of the membrane model over three cardiac cycles which was integrated to estimate the energetic cost of the cardiac action potential (figure 7*f* ). Note that energy continues to be consumed even during the resting state due to the presence of currents associated with ion transporters. Thus, while energy is predominantly consumed during the action potential, there is a positive gradient between action potentials (figure 7*f* ). By setting the energy consumption at the start of the second action potential to zero (figure 7*f* , dotted blue line), we calculated the energetic cost over the duration of the action potential to be 46.8 pJ. As the capacitive area of membrane for this model is 1.534 × 10 −4 cm 2 , the energy consumed per unit membrane area is 305 nJ cm −2 . When compared to Gawthrop *et al.*'s [21] estimate of 173 nJ cm −2 for the energetic cost of an action potential in the giant axon of a squid, the cardiac action potential uses 76% more energy. The main reason for this difference is that in contrast to a neuron, the cardiac action potential contains a plateau phase with opposing currents. Despite the relatively slow rate of change in voltage, the Ca 2+ and K + currents remain relatively high therefore, a large amount of energy is dissipated during the plateau phase.

#### (b) Chemostats influence the conserved moieties of cardiac action potential models

Because the earliest models of the cardiac action potential did not include active transporters, they used constant intracellular concentrations to maintain ionic gradients across multiple cardiac cycles [9,51]. Later models incorporated ion transporters, allowing them to represent physiological conditions with dynamic intracellular ion concentrations, and constant extracellular ion concentrations to model washout from the circulatory system [1,3]. Under ischaemic conditions, washout is greatly inhibited, thus models of ischaemia use dynamic extracellular ion concentrations [4]. We investigated the issue of drift in three classes of model: those with (A) dynamic ion concentrations on both sides of the membrane, representing models of myocytes under ischaemic conditions (B) dynamic intracellular ion concentrations but constant extracellular ion concentrations, representing models of myocytes under physiological conditions and (C) constant ion concentrations, representing models without transporters.

We used our bond graph model to represent these classes of models, selecting ions to fix at constant concentrations that resulted in three variants representative of the classes listed above. Thus, variant A represents models of cardiomyocytes under ischaemic conditions [4], variant B represents models of cardiomyocytes under physiological conditions [1,3] and variant C represents models without ion transporters [9,51]. Conserved moieties of each variant were found using the left nullspace matrix of the stoichiometric matrix (table 1), and these include for example, the total amount of K1 channel (moiety 1). Because the channel is neither synthesized nor degraded in our model, the total amount of channel, i.e. the sum of its closed (C_{K1}) and open (O_{K1}) states, remains constant over the course of a simulation.

Table 1. Conserved moieties associated with chemostat selection. Across some biochemical subgroups (moiety), there are constraints (conserved quantity) on a corresponding sum of species representing the total of the moiety. The conserved quantities remain constant over the course of a simulation. Q represents contributions of other species to charge imbalance across the membrane. The symbol *Σ* represents charge contributions from Markov states of channels and transporters. The definition of *Σ*, and all species can be found in the electronic supplementary material and code.

Similarly, moiety 10 for variant (A) represents the total amount of K + ions, which includes intracellular K + , extracellular K + and the K + ions bound to Na + /K + ATPase. The total amount of K + is constant when ion concentrations are dynamic. However, because fixing the concentration of K + requires an additional external flux, the conservation law is broken in variants (B) and (C). Because the membrane capacitance is included in the stoichiometry of the system, our method automatically identifies a charge conservation law (moiety 13 for variant (A), and moiety 10 for variant (B)).

While it is reassuring that our approach reveals the obvious conserved moieties described above, it also reveals the nontrivial charge conserved moiety that is missed by many existing cardiac electrophysiology models. The overall amount of intracellular charge can be described as a sum of contributions from intracellular K + , Na + , Ca 2+ (and its buffers) and Markov states from ion channels and transporters (*Σ*), similar to forms found in previous studies [10,52]. It should be noted, however, that when all ion concentrations were held constant, charge conservation was broken, as indicated by the absence of a conserved charge moiety in the bottom partition of table 1. In general, holding the concentration of a species constant breaks conservation laws [39] and the number of conserved moieties progressively decreases as more ion concentrations are modelled as chemostats. We discuss the consequences of this in latter sections.

#### (c) Non-conservative stimulus currents cause drift in models with a chargeconservation law

An important feature of cardiac electrophysiology models is that they must be simulated for extended periods to examine physiologically relevant changes in behaviour, thus we tested how the type of stimulus current affected each variant of the cardiac action potential model by pacing at 1 Hz for 30 min. As illustrated (figure 8*a*,*b*), a non-conservative stimulus resulted in drift when the model had dynamic ion concentrations either for all compartments, or only within the intracellular compartment. The drift was particularly pronounced when all ion concentrations were dynamic (figure 8*a*), as extracellular concentrations changed faster than intracellular concentrations. In contrast, the model was resistant to drift from a non-conservative stimulus when all ion concentrations were held constant (figure 8*c*).

Figure 8. Effect of stimulus type and variable ion concentrations on drift in the bond graph model. (*a*) Dynamic ion concentrations (*b*) Dynamic intracellular ion concentrations (*c*) Constant ion concentrations. Results are shown for stimuli that conserve overall charge (blue) and those that do not conserve charge (red). Non-conservative stimulus currents contain non-specific charge, and conservative stimulus currents contain charge and potassium ions. Charge values are given as differences from the initial value of − 5882.2 fmol . T = 310 K . Definitions: V dia , diastolic membrane potential APD, action potential duration at 90% repolarization. (Online version in colour.)

These results suggested that drift arose due to violations of the conserved charge moiety. Charge is a conserved moiety (table 1) in model variants where drift occurred with a non-conservative stimulus. In this situation, non-conservative stimulus currents cause drift because every stimulus causes a stepwise increase in the value of the conserved charge moiety (figure 8*a*,*b* bottom panels). However, because conservation laws are broken as more species are represented as chemostats [39], charge is no longer a conserved moiety when all ion concentrations are constant (table 1). Thus, an observation which may not be obvious to intuition is that under these conditions charge is no longer constant between stimuli, and therefore free to return to its original value after each stimulus (figure 8*c*, bottom panel), allowing such models to achieve a steady-state limit cycle.

Using the observations from the bottom row of figure 8, it is possible to develop a systematic and automated check for drift. Let *v*_{s} be a row matrix representing the stoichiometry of the stimulus current (with chemostats removed), *N* cd be the stoichiometric matrix after removing rows corresponding to chemostats, and *G* be the left nullspace matrix of *N* cd . As seen in this section, a stimulus current will cause drift if it results in any change to the conserved moieties *GX*. Therefore to avoid altering any of the conserved moieties, the stimulus current must have zero contribution to them, i.e. *Gv*_{s} = 0 (or equivalently, *v*_{s} needs to lie in the image of *N* cd ). Thus, the model drifts if *Gv*_{s}≠0.

#### (d) Initial conditions influence steady states through conserved moieties and chemostats

Next, for different sets of conserved moieties (as determined by constrained/dynamic ionic concentrations) we tested how the steady-state behaviour of the cardiac action potential was altered under three different initial conditions (figure 9). The first set of initial conditions (IC1) are common values for comparison (figure 5 at bottom). IC2 is the same as IC1 but with 1mM intracellular K + exchanged for 1 mM of intracellular Na + , such that charge is conserved but K + and Na + are not conserved. Similarly, IC3 is the same as IC1, but with some K + extruded and an equal amount of Na + moved into the cell such that charge, Na + , and K + are all conserved. When all ion concentrations are dynamic, IC1 and IC3 lead to the same steady state, but IC2 results in a different steady state (figure 9*a*). If only intracellular ion concentrations are dynamic, however, IC1 and IC2 result in identical steady states, but IC3 leads to a different steady state (figure 9*b*). Finally, keeping all ion concentrations constant leads to different steady states for all initial conditions (figure 9*c*).

Figure 9. Effect of initial conditions on steady-state behaviour. (*a*) Dynamic ion concentrations (*b*) Dynamic intracellular ion concentrations (*c*) Constant ion concentrations. The models were paced at 1 Hz for 30 min using a conservative stimulus current. [MgATP] = 6.95 mM, [MgADP] = 0.035 mM, [P_{i}] = 0.3971 mM, pH = 7.095, *T* = 310 K. Definitions: *V*_{dia}, diastolic membrane potential APD, action potential duration at 90% repolarization. (Online version in colour.)

These results demonstrate that the summed amount for each conserved moiety and/or chemostat value determines the steady-state behaviour of cardiac action potential models. To investigate this further, we calculated the values for conserved moieties and chemostats that resulted from each initial condition (table 2 *differences from IC1 indicated in italics*). For two sets of initial conditions to achieve identical steady states, all conserved moieties and chemostats must have the same value. Thus, under dynamic ion concentrations (figure 9*a*), IC3 results in the same steady state as IC1 because all conserved moieties have been preserved (table 2), whereas IC2 causes a different steady state because the *K* + and *Na* + conserved moieties take on different values. Similarly, when only intracellular ion concentrations are dynamic, IC2 preserves the value of all conserved moieties and chemostats, but IC3 changes the values of the chemostats corresponding to extracellular Na + and K + concentrations (table 2), hence the different steady state. When all ion concentrations were held constant, changes in the chemostat values (table 2) were associated with different steady states for all three initial conditions (figure 9*c*).

Table 2. The values of conserved moieties and chemostats under different initial conditions. All values are in fmol. Chemostats are indicated with (cs). Values different from IC1 are shown in italics.

### 4. Discussion

In this study, we developed a bond graph model of the cardiac action potential with the aim of resolving the issues of drift and non-unique steady states. Analysis using conserved moieties enabled the discovery of all conservation laws within the model, including the charge conservation law neglected by many existing studies. In addition to the conservation of charge law from previous studies [10,52,53], we found conservation laws corresponding to ions, states of Markov models of channels and transporters, and buffers, demonstrating the comprehensiveness of our approach. Because the bond graph approach requires species and processes to be resolved in biophysical detail, calculation of conservation laws is more straightforward than working purely with the mathematical equations of a model. Two key advantages of our approach over existing analyses are that it reveals all conservation laws in a comprehensive and systematic manner, and that it is general for all models of the cardiac action potential that can be represented as bond graphs. As a result, our approach can be scaled to more complex cardiac models without manual examination of the equations of a model to identify each individual conservation law, and the nontrivial charge conservation law will be accounted for when it appears as a conserved moiety. When simulated over long periods with a non-conservative stimulus, our bond graph model displayed solution drift, but it became resistant to drift when ion concentrations were held constant, demonstrating that changes in the value of a conserved charge moiety drive model drift. We also found that two sets of initial conditions can lead to different steady states if the values of their corresponding conserved moieties and chemostats are different, suggesting a strong link between conserved moieties and the steady-state limit cycles of cardiac action potential models. To demonstrate that our approach is general, we tested how the selection of chemostats (i.e. fixed concentrations) influenced drift and steady states by using variants of our model that were representative of existing models in the literature. Our approach highlights the subtle but relevant observation that holding ion concentrations constant changes the conserved moieties of the model, which in turn has an effect on the susceptibility of a model to drift and non-unique steady states. Because chemostats represent connections between a system and its external environment, they are essential to coupling together biological processes [23]. The coupling of biological processes generally causes changes in the conservation laws of a system which may be difficult to capture through observation. Our approach using bond graphs provides a systematic method of dealing with changes to conservation laws as a result of coupling models together.

#### (a) Drift

When paced with a non-conservative stimulus, variants of the model with a charge conservation law (A and B) underwent drift (figure 8*a*,*b*) consistent with previous studies on the Luo–Rudy model [10,11]. Our observations also explain why models with constant ion concentrations (similar to variant C) are more likely to be resistant to drift [10]. By observing changes in the charge conserved moiety, the bond graph approach attributes drift to regular perturbations in charge that cannot be restored due to the presence of a conservation law. Whereas previous analyses relied solely on intuition to derive a conservation law corresponding to charge [10,11], we note that our approach automatically derives this conservation law, and can also detect other conservation laws that may be relevant for drift.

We note that in order to avoid drift, all conserved moieties (and not only the conserved moiety corresponding to charge) must be preserved by the stimulus current. However, in the examples explored in this study, we found that the stimulus currents preserved the value of all conserved moieties apart from charge, and therefore were not plotted. While intuition may suggest that modelling a stimulus current with K + would cause an accumulation of intracellular K + , the K + is passively transported back through K + channels. Furthermore, the total amount of K + (the relevant conserved moiety) remains constant because the loss of K + from the extracellular side is offset by accumulation on the intracellular side. By contrast, for charge conservation in variants A and B, the model contains no mechanism to reverse the change in charge difference caused by the non-conservative current, thus drift results.

As demonstrated, the bond graph method requires construction of a stoichiometric matrix, providing a simple approach to check whether a stimulus current will cause drift. While it is common practice to use K + as the charge carrier for stimulus currents, it is likely that multiple species contribute to the current [8,10]. Thus, the automated approach suggested here is well suited for checking whether more complex stimulus currents satisfy conservation of charge, as well as other conservation laws within the model. It should be noted, however, that while a model satisfying *Gv*_{s} = 0 will not drift due to violating conservation laws, drift may still occur due to an imbalance of currents throughout the action potential, for instance, in the absence of Na + /K + ATPase, the ionic gradients would gradually disappear in a model with dynamic ion concentrations.

Finally, we believe that this analysis provides a link between the issues of drift and steady states. Our models show that drift due to a non-conservative stimulus current can be attributed to changes in the value of the charge conserved moiety with every stimulus, and accordingly the steady state of the model changes. Model drift then occurs as the solution continually chases a moving steady state.

#### (b) Effects of initial conditions on steady states

We also found that initial conditions of cardiac action potential models change their steady states through the values of chemostats and conserved moieties (figure 9 and table 2). Accordingly, the same perturbation to initial conditions can have different effects on the steady state depending on which species are held constant. Therefore, in addition to ensuring that the concentration of ions is physiological, care should be taken to correctly initialize each state of buffers and Markov models of ion channels and ion transporters, as they may contain a significant fraction of total ion abundance. For example, Ca 2+ buffers and SERCA can sequester a significant amount of Ca 2+ and they should be initialized with the correct amount of bound Ca 2+ when multi-state models are used [54]. We note that the difficulty of manually deriving conservation laws increases exponentially as models of cardiac electrophysiology become more complex, and we believe that our approach extends on existing analyses [10,11] to provide a general method for assessing steady-state behaviour by comparing the values of conserved moieties and chemostats that result from each initial condition.

In the field of biochemical network analysis, there is a well-established dependence of quiescent steady states on conserved moieties, and numerous mathematical techniques for assessing the uniqueness and stability of these steady states have been developed [55,56]. However, the influence of conserved moieties on limit cycles in an oscillating system that is regularly stimulated has yet to be investigated. Our results hint at similarities between these two fields, and while we only tested the uniqueness of steady states using relatively small perturbations to the initial conditions, it is possible that a set of conserved moieties may have multiple steady states, and greater perturbations may lead to other limit cycles.

#### (c) The ‘differential’ and ‘algebraic’ methods

The discovery of conservation principles in cardiac electrophysiology has led to a debate over whether to use the differential or algebraic methods of simulation [7,10,11,52,53]. The differential method is the calculation of membrane voltage by integrating total current, and the algebraic method is the calculation of membrane voltage using an algebraic relationship derived from charge conservation. We chose the differential method over the algebraic method because it better supports model reuse and modularity - in particular it is easier to modify the equations to select different species as chemostats, and to combine equations when two models are coupled. We note, however, that the algebraic method may reduce computational complexity [10,20]. In bond graph modelling, the algebraic method can be implemented by using conserved moieties to turn the system of ODEs into an index-0 differential algebraic equation (DAE) (eqn 3.48 of [20]). This method generalizes existing algebraic methods to reduce the system of differential equations by using all conserved moieties and not just the conserved charge moiety. While we did not use the algebraic approach, we emphasize that the choice of method relates to numerical approaches for model simulation rather than the underlying physics of the system [10]. Therefore, the differential and algebraic methods are equivalent in conservative systems provided that the initial conditions and values of conserved moieties are consistent.

#### (d) Integration into whole-cell models

Our bond graph model of the cardiac action potential is the first step towards a fully integrated whole-cell bond graph model of a cardiomyocyte that couples electrophysiology, signalling, metabolism and mechanics. As energy drives all biological processes, and energy supply can be limited under certain pathophysiological states, it is of interest to examine how cardiomyocytes allocate their energy, and also to estimate the efficiency of processes that are essential for cardiomyocyte function [57,58]. Modelling studies for the energetic regulation of a cardiac cell exist across the literature [58], but while some components used in these models are thermodynamically consistent [59,60], existing whole-cell models are neither energy-based nor thermodynamically consistent throughout the entire model. Furthermore, because existing experimental and modelling studies use ATP consumption as a proxy for energy consumption, they can only estimate the energy consumption of major energy sinks: the Na + /K + ATPase, SERCA and crossbridge cycling [58,61]. A bond graph approach may thus provide more detailed insights into how a cardiac cell uses energy downstream of ATP hydrolysis processes, and help to identify energy-consuming processes. Because the bond graph approach is energy-based, it not only provides the necessary constraints to develop a thermodynamically consistent model, but also allows us to directly assess energy consumption of the model (in Joules). We found that when normalized against membrane area, the cardiac action potential consumes approximately 76% more energy than an action potential in the axon of a giant squid. To the authors' knowledge, this is the first account of energy consumed by electrochemical processes during the cardiac action potential.

#### (e) Limitations

A limitation of our approach is that the physical constraints imposed by the bond graph approach prevent the direct translation of existing models of cardiac electrophysiology into bond graphs. This points at thermodynamic issues that exist within existing models of cardiac electrophysiology while many models contain components that are thermodynamically consistent, the authors are not aware of any models that are entirely thermodynamically consistent. Nonetheless, the inability to perfectly replicate properties of existing models as bond graphs impedes the creation of bond graph models of cardiac electrophysiology, as it requires new models to be built from the bottom up with existing equations only used as guides for parameter fitting. As a result of this limitation, our bond graph model represents only a subset of the ion transport processes within a cardiomyocyte, and only approximates the behaviour described in existing models. While our model is able to reproduce many essential features of the action potential, it is unlikely to be physiologically realistic under conditions different to those used for our fitting process. Thus in future studies, it would be interesting to explore the use of more complex bond graph components to generate more physiologically realistic models. We chose to use a restricted set of components and constitutive equations to keep the model simple, although the bond graph framework is flexible enough to account for a wider range of equations provided that they are thermodynamically consistent. While many ion channels (including those for *I*_{Na}, *I*_{K1} and *I*_{Kp} in the Luo and Rudy model) are modelled using linear I-V relationships, the choice of I-V relationship is generally chosen on the basis of providing an empirical fit to data rather than as a result of physical principles. Thus, a more ideal approach would be to use nonlinear bond graph components derived from physical principles, and fit their parameters directly to experimental data. Additionally, empirical equations for gating transition parameters are generally not expressed in a thermodynamic framework, and are therefore impossible to replicate exactly with a bond graph model. We used a simple gating mechanism to reduce equation complexity while maintaining thermodynamic consistency however, the quality of fit could be improved by using more complicated gating mechanisms, or other thermodynamically consistent constitutive equations.

Because of physical restrictions imposed by the bond graph framework, we were forced to model ion channels and transporters using Markov states to faithfully represent their underlying physics. However, this produced a model that had numerous states compared to the number of biological processes. While it is reassuring to find that our method of identifying conserved moieties remained robust despite this complexity, simulation of the model was computationally expensive. For the purpose of integrating this action potential model into a larger whole-cell model, it would be useful to have simple model components that reduce computational cost. While current methods for reducing biochemical models in the bond graph framework are not advanced enough to apply to the biological components in this study, we note that bond graphs provide a useful foundation for applying model simplification while ensuring that thermodynamic consistency is maintained [20].

We also decided to limit the transport processes included in our model to those considered essential for producing a cardiac action potential, while maintaining a limit cycle using dynamic ion concentrations. Our bond graph model omitted many ionic currents due to their small amplitudes however, these channels may have greater contributions under conditions which vary from those tested here. Thus, an obvious extension of this work would be the integration of other electrogenic processes within the cardiac membrane. It would be interesting to investigate whether coupling other models requires further tuning of parameters [62], and whether the presence of physical bond graph parameters changes this process. A related limitation is that our model does not account for ion concentrations such as *H* + and *Cl* − , as well as pH buffers [5]. While including these ions and their transporters would lead to a more accurate model, the omission of these ions did not cause inconsistencies in the conservation laws described in this study. These ions were assumed to be membrane impermeable in our model, and thus their constant contributions to the membrane potential were accounted for in the initial value of the charge conserved moiety. Similar to how calcium and its buffers were accounted for in our current approach (table 1), our bond graph approach is sufficiently general to account for other ions and their buffers.

When formulating the structure and parameters for a bond graph model of the cardiac action potential (or most other biological processes), it is possible to either fit against existing mathematical models or the underlying experimental measurements. For all processes in this study excluding the NCX, we developed our bond graph model to reproduce the behaviour of an existing model, in an attempt to re-use existing knowledge about these processes. This approach poses constraints on the bond graph structure used, especially for gating structure. Therefore, it would be interesting to develop an approach that assesses bond graph structures as well as bond graph parameters, based on their fits to data [62]. Such an approach may provide a better fit to the data, and uncover insights into the physical mechanisms of ion channels.

### 5. Conclusion

In this study, we have developed a bond graph model of the cardiac action potential and used this to explore the issues of drift and non-unique steady states. We demonstrate that the analysis of conserved moieties generalizes the concept of charge conservation used in earlier studies, and found that changes in conserved moieties can explain drift as well as changes in steady-state behaviour. Importantly, holding ion concentrations constant can have significant consequences on both drift and steady states as they change the conserved moieties in the model. Our approach to resolving drift and non-unique steady states is sufficiently general that it can be applied to any bond graph model of the cardiac action potential model. We hope that the bond graph approach outlined here will prove useful for the development of future cardiac electrophysiology models, and eventually whole-cell models of the cardiomyocyte.

### Data accessibility

The code associated with this study is available from GitHub (https://github.com/uomsystemsbiology/bond_graph_cardiac_AP), and archived on Zenodo (https://doi.org/10.5281/zenodo.1172205) [63]. The repository contains Matlab (The MathWorks, Natick, MA) code that generates the figures, CellML code containing parameters, initial conditions and equations of the model, and full details of the bond graph structure.

### Author's contributions

M.P., P.J.G., J.C. and E.J.C. developed the theory. M.P. performed the research. K.T. provided conceptual advice and helped interpret the results. All authors contributed to the text of the manuscript and gave final approval for publication.

### Competing interests

We have no competing interests.

### Funding

M.P. acknowledges the financial support provided by an Australian Government Research Training Program Scholarship. P.J.G. thanks the Melbourne School of Engineering for its support via a Professorial Fellowship. K.T. is supported by the Heart Foundation of New Zealand (Research Fellowship 1692) and the Marsden Fund Council from Government funding, managed by Royal Society Te Ap a ¯ rangi (Marsden Fast-Start UOA1703. This research was in part conducted and funded by the Australian Research Council Discovery Projects funding scheme (project DP170101358).

## Introduction

Neurons communicate via chemical and electrical signals, and an integral part of this communication is the transmission of action potentials along their axons. The velocity of action potentials is crucial for the right timing in information processing and depends on the dynamics of ion channels studding the axon, but also on its geometrical properties. For instance, the velocity increases approximately linearly with the diameter of myelinated axons [1]. Myelin sheaths around axons are an evolutionary trait in most vertebrates and some invertebrates, which developed independently in several taxa [2]. The presence of a myelin sheath increases the velocity of action potentials by enabling saltatory conduction [3]. Long-term, activity-dependent changes in the myelination status of axons are related to learning [4]. The functional role of differentiated myelination is to regulate and synchronise signal transmission across different axonal fibres to enable cognitive function, sensory integration and motor skills [5]. White-matter architecture has also been found to affect the peak frequency of the alpha rhythm [6]. Axons and their supporting cells make up the white matter, which has, for a long time, only been accessible to histological studies [7, 8]. With the advent of advanced MRI techniques, some of the geometric parameters of axonal fibre bundles have become accessible to non-invasive methods. Techniques have been proposed to determine the orientation of fibre bundles in the white matter [9] as well as to estimate the distribution of axonal diameters [10], the packing density of axons in a fibre bundle [11, 12], and the ratio of the diameters of the axon and the myelin sheath (g-ratio) [13].

First quantitative studies were done by Hursh [14] who established the (approximately) linear relationship between action potential velocity and axonal radius in myelinated axons, and Tasaki [3] who first described saltatory conduction in myelinated axons. Seminal work on ion channel dynamics was later done by Hodgkin and Huxley, establishing the voltage-dependence of ion channel currents [15]. The general result of voltage-dependent gating has been confirmed in vertebrates [16], yet a recent result for mammals suggests that the gating dynamics of sodium channels is faster than described by the original Hodgkin-Huxley model, thereby enabling faster generation and transmission of action potentials [17]. In general, parameters determining channel dynamics differ widely across neuron types [18].

Seminal studies into signal propagation in myelinated axons using computational techniques were done by FitzHugh [19] and Goldman and Albus [20]. Goldman and Albus gave the first computational evidence for the linear increase of the conduction velocity with the radius of the axon, provided that the length of myelinated segments also increases linearly with the axonal radius. The linear relationship is supported by experimental evidence [21], although other studies suggest a slightly nonlinear relationship [22]. More recently, computational studies have investigated the role of the myelin sheath and the relationship between models of different complexity with experimental results [23]. One of the key findings here was that only a myelin sheath with finite capacitance and resistance reproduced experimental results for axonal conduction velocity. Other studies investigated the role of the width of the nodes of Ranvier on signal propagation [24, 25], or the effect of ephaptic coupling on signal propagation [26–33].

Most computational studies employ numerical schemes, i.e. they discretise the mathematical problem in space and time and use numerical integration methods to investigate the propagation of action potentials. One problem that arises here is that the spatial discretisation must be relatively coarse to ensure numerical stability, which can be remedied to some extent by advanced numerical methods and computational effort [34]. The other problem, however, cannot be remedied that easily: it is the lack of insight into how the model parameters influence the results, since there is a large number of parameters involved. A way to illustrate parameter dependencies in an efficient manner is to use analytical techniques all the while simplifying the model equations and extracting essential features. Studies that use analytical methods are few and far between [35–38] yet it is also worth noting that from a mathematical perspective, myelinated axons are similar to spine-studded dendrites, in the sense that active units are coupled by passive leaky cables. An idea that we pick up from the latter is to simplify the ionic currents crossing the membrane [39, 40], there at dendritic spines mediated by neurotransmitters, here at nodes of Ranvier mediated by voltage-gated dynamics.

The goal of this article is to use analytical methods to study the influence of parameters controlling action potential generation, and geometric and electrophysiological parameters of the myelinated axon, on the speed of action potentials. The main focus here is on parameters determining the axonal structure. This will be achieved by replacing the Hodgkin-Huxley dynamics with a spike-diffuse-spike model for action potential generation, i.e. ion currents are released at nodes of Ranvier when the membrane potential reaches a certain threshold. These ion channel currents are considered voltage-independent, but we investigate different forms of currents, ranging from instantaneous currents to currents that incorporate time delays. We also investigate ion currents that closely resemble sodium currents measured experimentally. Our aim is to derive closed-form solutions for the membrane potential along an axon, which yields the relationship of action potential velocity with model parameters.

The specific questions we seek to answer here are the following. First, we query how physiological parameters can be incorporated into our mathematical framework, especially parameters that control the dynamics of the ionic currents. We test if parameters from the literature yield physiologically plausible results for the shape and amplitude of action potentials, and test how the ionic currents from multiple nearby nodes of Ranvier contribute to the formation of action potentials. Secondly, we ask how geometric parameters of an axon affect the transmission speed in a single axon, and how sensitive the transmission speed is to changes in these parameters. We seek to reproduce known results from the literature, such as the dependence of the velocity on axon diameter. We also explore other dependencies, such as on the g-ratio, and other microscopic structural parameters resulting from myelination. We compare the results of our spike-diffuse-spike model with the results from a detailed biophysical model recently used to study the effect of node and internode length on action potential velocity [24]. Thirdly, we investigate how ephaptic coupling affects the transmission speed of action potentials, and what the conditions are for action potentials to synchronise. In particular, we examine how restricted extra-axonal space leads to coupling between two identical axons, and how action potentials travelling through the coupled axons interact.

## Proton-Dependent Inhibition, Inverted Voltage Activation, and Slow Gating of CLC-0 Chloride Channel

CLC-0, a prototype Cl − channel in the CLC family, employs two gating mechanisms that control its ion-permeation pore: fast gating and slow gating. The negatively-charged sidechain of a pore glutamate residue, E166, is known to be the fast gate, and the swinging of this sidechain opens or closes the pore of CLC-0 on the millisecond time scale. The other gating mechanism, slow gating, operates with much slower kinetics in the range of seconds to tens or even hundreds of seconds, and it is thought to involve still-unknown conformational rearrangements. Here, we find that low intracellular pH (pH_{i}) facilitates the closure of the CLC-0’s slow gate, thus generating current inhibition. The rate of low pH_{i}-induced current inhibition increases with intracellular H + concentration ([H + ]_{i})—the time constants of current inhibition by low pH_{i} = 4.5, 5.5 and 6 are roughly 0.1, 1 and 10 sec, respectively, at room temperature. In comparison, the time constant of the slow gate closure at pH_{i} = 7.4 at room temperature is hundreds of seconds. The inhibition by low pH_{i} is significantly less prominent in mutants favoring the slow-gate open state (such as C212S and Y512A), further supporting the fact that intracellular H + enhances the slow-gate closure in CLC-0. A fast inhibition by low pH_{i} causes an apparent inverted voltage-dependent activation in the wild-type CLC-0, a behavior similar to those in some channel mutants such as V490W in which only membrane hyperpolarization can open the channel. Interestingly, when V490W mutation is constructed in the background of C212S or Y512A mutation, the inverted voltage-dependent activation disappears. We propose that the slow kinetics of CLC-0’s slow-gate closure may be due to low [H + ]_{i} rather than due to the proposed large conformational change of the channel protein. Our results also suggest that the inverted voltage-dependent opening observed in some mutant channels may result from fast closure of the slow gate by the mutations.

## Abstract

Mathematical models of cardiac action potentials have become increasingly important in the study of heart disease and pharmacology, but concerns linger over their robustness during long periods of simulation, in particular due to issues such as model drift and non-unique steady states. Previous studies have linked these to violation of conservation laws, but only explored those issues with respect to charge conservation in specific models. Here, we propose a general and systematic method of identifying conservation laws hidden in models of cardiac electrophysiology by using bond graphs, and develop a bond graph model of the cardiac action potential to study long-term behaviour. Bond graphs provide an explicit energy-based framework for modelling physical systems, which makes them well suited for examining conservation within electrophysiological models. We find that the charge conservation laws derived in previous studies are examples of the more general concept of a ‘conserved moiety’. Conserved moieties explain model drift and non-unique steady states, generalizing the results from previous studies. The bond graph approach provides a rigorous method to check for drift and non-unique steady states in a wide range of cardiac action potential models, and can be extended to examine behaviours of other excitable systems.

### 1. Introduction

Models of the cardiac action potential have been developed to study cardiac diseases, such as arrhythmia [1–3], ischaemia [4] and acidosis [5]. Increasing model complexity has led to concerns over the occurrence of drift and non-unique steady states [6–8], particularly for extensions of the DiFrancesco & Noble [9] and Luo–Rudy [1,3] models. While solutions to these issues have been proposed using conservation principles [10,11], they have not been universally applied for more recent models, many of which still use non-conservative stimulus currents that predispose them to drift [12–14]. More recently, the Food and Drug Administration (FDA) has initiated plans to use cardiac action potential models to assess potential drug side-effects on cardiac instability through the human *ether-à-go-go*-related gene (hERG) K + channel. Thus, with an increasing emphasis on model robustness and accuracy, there is a renewed incentive to resolve the issues of drift and non-unique steady states [15,16].

Drift is the failure of a model to reach a consistent limit cycle when simulated over long periods, and is often caused by a non-conservative stimulus containing current with no charge carrier [6,10]. Hund *et al.* [10] derived a charge conservation law, and found that non-conservative stimulus currents violate this conservation law, hence they proposed K + ions as the current charge carrier to resolve this. A related issue in many models where drift has been resolved is that steady-state limit cycles under constant pacing depend upon the initial conditions and are therefore non-unique [7,8,10]. Thus, depending on the initial conditions, the same model may lead to different conclusions. Like drift, authors have suggested that charge conservation can constrain initial conditions such that they lead to the same steady state [8,10,11].

While the studies by Hund *et al.* [10] and Livshitz & Rudy [11] suggest measures to eliminate drift and attain a unique steady state by using conservation laws, their analyses are limited in their scope and not a comprehensive solution for all models. Because existing studies [10,11] explore charge conservation only in specific models, and the conservation laws were derived from physical intuition rather than a principled mathematical approach, it is difficult to generalize their findings to other models where charge conservation is routinely neglected. Furthermore, because these studies focus only on conservation of charge, they may miss other conservation laws relevant for long-term behaviour, such as those corresponding to ions, ion channels and buffers. A general approach is, therefore, desirable to deal with the issues of drift and steady states in a more systematic manner and for a broader range of models.

To facilitate a general approach, we propose the use of bond graphs which explicitly model energy transfer across physical systems to ensure compliance with conservation principles. Bond graphs were initially invented to model hydroelectric systems [17] and they have subsequently been extended to model chemical [18], biochemical [19,20] and electrochemical systems [21]. As with all physical systems, biological processes must obey the fundamental principles of physics and thermodynamics [22] therefore, bond graphs are well suited for constraining models of biological systems to physically plausible solutions [23], and also for inferring the energetic cost of biological processes [21,23–25]. Because the bond graph representation emphasizes analogies between different physical domains, electrophysiological systems can be analysed as an analogous biochemical system with a stoichiometric matrix that describes the stoichiometry of each reaction within its columns [20,26–28]. In this context, the ‘conservation principle’ described in earlier studies is an example of the more general principle of a conserved moiety in metabolic and bond graph analysis [23,29].

In this study, we develop an abridged bond graph model of the cardiac action potential and outline a general approach to study the effects of conserved moieties on drift and steady-state behaviour. Our bond graph model simulates the essential features of the cardiac action potential, and because bond graphs are energy-based this easily provides an estimate of the energetic cost (in Joules) of the cardiac action potential. Our analysis reveals conservation of charge as one of the conserved moieties of our model, along with other conserved moieties corresponding to ions, channels, transporters and buffers. We observed that our model solution was subject to drift when the stimulus current violated any conservation laws corresponding to the conserved moieties, and that changes to the initial conditions led to different steady states if the value of any conserved moiety was changed. To demonstrate that our approach is general, we analyse variants of our bond graph model where different ions have been fixed at a constant concentration (corresponding to ‘chemostats’). It should be noted that fixing an ion concentration can change the conserved moieties of a system, therefore influencing a model's susceptibility to drift and non-unique steady states. The bond graph approach is a useful and general method to identify and interpret conservation principles, and it can link conserved moieties to individual steady states. Our approach can be used to automatically derive charge conservation laws that are frequently neglected in existing models of cardiac electrophysiology, and we build upon existing reports [10,11] to propose solutions for drift and non-unique steady states which work for all cardiac action potential models that can be represented using bond graphs.

### 2. Methods

#### (a) Model components

To study the issues of drift and non-unique steady states, we built an abridged bond graph model of the cardiac action potential, with the minimal number of channels and pumps required to simulate the essential features of a cardiac action potential, and maintain ionic concentrations over long periods of simulation. We note that a series of mathematical equations can only be described using bond graphs if they describe a thermodynamically consistent system, and that no existing model of cardiac electrophysiology is entirely thermodynamically consistent. As seen in Gawthrop *et al.* [21], an exact translation of existing ion channel models into the bond graph framework is generally not possible. Therefore, for some transport processes and particularly for ion channels, it was necessary to build new models under the bond graph framework, although their parameters and structure can be chosen based on equations in existing models. Accordingly, many of the essential components of our model were based primarily on the Luo–Rudy 1994 dynamic model [1], although it is possible to use other models and/or model more sub-cellular processes. We argue in §2c that some aspects of the equations for the ion channels described in the Luo and Rudy model are not thermodynamically consistent, and unable to be represented as bond graphs. We describe our approach to modelling ion channels in §2d. As we use simple and elementary bond graph structures based on physical principles, our bond graph representation of ion channels is a physically constrained approximation of the equations described in Luo and Rudy, although it is still able to simulate the essential features of the cardiac action potential.

Model components are shown in figure 1*a*, together with the overall bond graph structure (figure 1*b*). Ion channels and Ca 2+ buffering components were based upon their representations in Luo & Rudy [1]. The L-type Ca + channel in the Luo–Rudy model is permeable to calcium, sodium and potassium, but we neglected its sodium conductance as this has a relatively small contribution to the action potential. The Na + /K + ATPase model was based on the model by Terkildsen *et al.* [4], with modifications suggested by Pan *et al.* [30] to allow conversion into a bond graph model. The equation for the Na + –Ca 2+ exchanger (NCX) current in Luo and Rudy did not have an obvious correspondence to a bond graph structure, thus we modelled this component using a simple bond graph module that was fitted to experimental data [31,32]. The general approach for modelling ion channels is described later in the methods, with details on the other components, as well as further detail on ion channel modelling given in the electronic supplementary material.

Figure 1. Action potential model. (*a*) Cell schematic, (*b*) overall bond graph structure. The bond graph modules Na_channel, NaK, K_channels, LCC, NCX and Ca_buffer contain more detailed aspects of the bond graph structure which are described further in the Methods and electronic supplementary material. Coloured bonds link bond graph modules to the appropriate chemical species. Definitions: *I*_{Na}, sodium current *I*_{K1}, time-independent K + current *I*_{K}, time-dependent K + current *I*_{Kp}, plateau K + current *I*_{LCC}, L-type Ca 2+ current NCX, Na + –Ca 2+ exchanger Na/K, Na + /K + ATPase TRPN, troponin CMDN, calmodulin. (Online version in colour.)

#### (b) Bond graph modelling

Here, we briefly outline bond graph components as used in electrophysiological modelling. For a more comprehensive introduction, the texts by Gawthrop & Smith [33] and Borutzky [34] provide detailed descriptions of bond graph theory, and Gawthrop & Bevan [35] provide a short tutorial for engineers. Theory for bond graph modelling of biochemical systems can be found in [19,20,23,36].

Bond graphs consist of components (representing physical objects and processes), bonds (representing the transfer of energy) and junctions (representing network structure). Each bond carries two variables: an effort *e* and a flow *f* , such that their product determines the power of the bond (i.e. *p* = *ef*). Thus, bond graphs explicitly account for energy transfer, and are thermodynamically consistent. Because effort and flow are generalized variables, they can represent quantities from a variety of physical systems, including mechanical (*e* = force (N), *f* = velocity (m s −1 )), electrical (*e* = voltage (V), *f* = current (A)) and hydraulic systems (*e* = pressure (Pa), *f* = volumetric flow rate (m 3 s −1 )) [34].

The network structure of a bond graph is specified by 0 and 1 junctions. The 0 (or effort) junctions specify that efforts of all connected bonds are equal, and thus to ensure conservation of energy through this junction, the flows of the bonds must sum to zero. In the electrical and hydraulic domains, 0 junctions represent parallel connections, whereas they represent series connections in the mechanical domain. By a similar principle, 1 (or flow) junctions specify that the flows of all connected bonds are equal, ensuring that their efforts sum to zero. Thus, 1 junctions correspond to series connections in the electrical and hydraulic domains, and parallel connections in the mechanical domain.

To illustrate the use of a bond graph for electric circuit analysis, we consider the electric circuit where two capacitors are connected to a resistor in series (figure 2*a*). All components are linear, described by the following equations:

Figure 2. Conceptual representations of physical systems. (*a*) A bond graph for the illustrated simple electric circuit with two capacitors and a resistor in series. (*b*) A bond graph analogous to the electric circuit in (*a*) can also represent the chemical reaction A ⇌ B . (*c*) Bond graphs can also model interaction of components in both the chemical and physical domains, such as the transport of an ion across a membrane. (*d*) Transport of an ion across a membrane through an ion channel involves gating which modulates the rate of reaction. Thus the ion channel is analogous to a potentiometer. (Online version in color.)

More recently, bond graphs have been extended to model biochemical systems [19,20] where the chemical potential *μ* (J mol −1 ) is the effort variable, and molar flow rate *v* (mol s −1 ) is the flow variable. As temperature and pressure are assumed to be constant in biochemical systems, the measure of thermodynamic potential *μ* corresponds to Gibbs free energy [37]. For bond graph modelling of more general systems where temperature and pressure change, the reader is directed to the text by Thoma & Bouamama [38].

Because mass action equations are nonlinear in general, nonlinear components are required to model biochemical systems. Each chemical species is represented as a capacitor. However, in contrast to the electrical domain, the constitutive equation for the capacitor representing each species is logarithmic:

The bond graph framework for biochemistry can be extended to electrochemical systems [21] as demonstrated in figure 2*c*, which models the transport of a positively charged species *X* across a membrane. It should be noted that chemical species are described with C components that have a logarithmic association, whereas the C component corresponding to the (electric) membrane potential has a linear constitutive relationship. A transformer (TF) is used to convert the membrane voltage into an equivalent chemical potential through Faraday's constant *F* = 96485 C mol −1 , such that:

We chose to represent ion channels such that conductance was modulated by membrane voltage, both directly and indirectly through gating processes. A bond graph representation for this relationship is given in figure 2*d*. As shown, this model has the same electrical representation as figure 2*c* however it uses a variable resistor. The bond graph representation contains the same states, with C : xi, C : xe, and C : mem (with a transformer) connected through 0 junctions. In this case however, the Re components that describe the constitutive relation have been changed, such that Re_GHK : r1 is connected to an additional effort that modulates its velocity, and the gating affinity *A* g is added to both the forward and reverse affinities to describe changes in permeability due to gating. Further detail on modelling ion channels using bond graphs is given in §2d.

#### (c) Modelling approach

Because bond graphs constrain the equations of a model to ensure thermodynamic consistency, they only allow the representation of the subset of general mathematical models that are thermodynamically consistent. In such models, dissipative processes such as reactions and ionic currents can only proceed in the direction of decreasing chemical or electrochemical potential, and when there is no potential gradient (i.e. the process is at equilibrium), the process must stop. However, many existing models of cardiac electrophysiology are not constrained to ensure this behaviour, and in some cases, they describe physically infeasible systems that create energy out of nowhere [20]. Therefore, many existing models do not have a direct bond graph representation [23]. An advantage of the bond graph approach is that the discipline required to convert an existing model into a bond graph helps to highlight thermodynamic issues and inconsistencies in existing models that would have otherwise been missed or ignored (see Gawthrop *et al.* [23] for an example of this process applied to a model of glycolysis). For the currents in this study, the Luo–Rudy equations for the time-dependent K + and L-type Ca 2+ channels are not thermodynamically consistent, and therefore cannot be directly converted into a bond graph model. As shown in equation (2.11), a consequence of the thermodynamic consistency of the bond graph framework is that there is a constraint on the equilibrium point for each ion channel that is determined by the Nernst equation. As the equilibrium points of the time-dependent K + and L-type Ca 2+ channels do not correspond to the Nernst potential (see electronic supplementary material, A.1.3 and A.1.5), they are not thermodynamically consistent and therefore their equations cannot be translated into a dissipative bond graph component. These issues were addressed by finding a thermodynamically consistent approximation that satisfies the Nernst equation. A second issue is that it is difficult to simultaneously model open-channel currents and channel gating using the elementary components, described in §2b and figure 2, thus we approximated these relations using the components described in that section rather than using more complicated components. Therefore, rather than attempting to reproduce the Luo–Rudy equations exactly, we built a bond graph structure as implied by the equations in the Luo and Rudy model, and chose parameters of our bond graph model to fit aspects of the Luo–Rudy model as closely as possible, specifically the current-voltage (I-V) curves and gating parameters. For all other components conversion into a bond graph model was more straightforward, and we used the methods of Gawthrop *et al.* [23]. Further information on the bond graph model, and parameter identification is given in the electronic supplementary material.

#### (d) Ion channel modelling

##### (i) Bond graph structure.

In this section, we discuss decisions made in developing models of ion channels. The bond graph structure for the Kp channel is shown in figure 3. The other channels have similar structures that follow from the discussion in this section.

Figure 3. The bond graph model of a plateau K + channel. (*a*) The **channel_Kp** module describes the current through the ion channel. (*b*) The **gate_en_Kp** contains the states required for gating. (*c*) The **vRe** module contains a voltage-dependent reaction used to describe channel state transitions. (*d*) The channel current and gating modules are combined into an ion channel model ( **Kp_channel** ). (Online version in colour.)

##### (ii) Current–voltage relations.

While thermodynamic properties can be used to determine how membrane voltage and ionic concentrations relate at equilibrium, they do not specify behaviour away from equilibrium. For this purpose, the current-voltage (I-V) relationship defines how the membrane voltage relates to the current through a specific channel. In many electrophysiology models [1], currents through ion channels are described through the use of linear I-V relationships based on empirical fits to data rather than fundamental physical principles. Using bond graphs, it is difficult to incorporate the effects of gating using a linear I-V equation. While linear I-V equations can be thermodynamically consistent provided they satisfy the Nernst equation in equation (2.11), and such equations can be represented by using a linear R component, the modulation of its conductance would require the use of signal bonds that do not necessarily represent physical processes [34]. Therefore, we use the Goldman–Hodgkin–Katz (GHK) equation to model ion channels, as it enables relatively simple incorporation of ion channel gating as a physics-based biochemical module [21]. The GHK equation defines a nonlinear relationship between current *I* and membrane voltage *V* :

As many ion channels in the Luo–Rudy model are described using a linear I-V relationship, the use of GHK equations requires approximation. A comparison of the resulting I-V curves is given in figure 4. Details of the fitting process, as well as the fitted bond graph parameters, are given in the electronic supplementary material. The Na + channel I-V curves appeared to match reasonably well (figure 4*a*), with some discrepancies at positive membrane potentials. For K + channels (figure 4*b*–*d*), we attempted to optimize the fit across voltages that correspond to their physiological function, so that their currents would be most similar to those of Luo and Rudy when the ion channels are open. Accordingly, for *I*_{K1} (−90 mV≤*V* < − 30 mV), *I*_{K} (−20 mV≤*V* ≤30 mV) and *I*_{Kp} (*V* > 0 mV) the I-V curves matched reasonably well in these regions. Discrepancies occurred outside these ranges of voltages, but appeared to only cause minor differences to the currents. In their implementation of *I*_{K}, Luo & Rudy [1] use a thermodynamically inconsistent I-V equation where the current is non-zero at the Nernst potential for K + . Despite this, bond graph parameters could still be chosen to give a reasonable fit to this I-V equation (figure 4*c*). Because the Luo–Rudy model based their L-type Ca 2+ I-V curves on the GHK equation, there was a far closer match between the bond graph and Luo–Rudy models for these currents, (figure 4*e*,*f* ) and the K + curve was matched exactly (figure 4*f* ).

Figure 4. Comparison of I-V curves between the Luo–Rudy (LRd) and bond graph (BG) models. (*a*) *I*_{Na}, (*b*) *I*_{K1}, (*c*) *I*_{K}, (*d*) *I*_{Kp}, (*e*) *I*_{Ca,L}, (*f*) *I*_{K,L}. The standard concentrations in Luo & Rudy [1] ([Na + _{i}] = 10 mM, [Na + _{e}] = 140 mM, [ K + _{i} ] = 145 mM , [K + _{e}] = 5.4 mM, [Ca + _{i}] = 0.12 μM, [Ca + _{e}] = 1.8 mM) were used to match I-V curves.

##### (iii) Modulation.

While the I-V curves describe currents through open ion channels, a formulation for gating is required to describe the number of open ion channels at any given time. In the Hodgkin–Huxley framework, gating is modelled as differential equations that give the proportion of open gates at any given time. We incorporated the effects of gating through a gating affinity *A* g , which is added to both the forward and reverse affinities of a reaction (figure 3*a*) to modulate its rate without changing the equilibrium [21].

##### (iv) State models.

Ion channel models must account for gating and bond graphs require the use of physical components to achieve this. We model gating as transitions between channel states, known in the literature as Markov models [41,42]. To illustrate, we use the example of a typical Na + channel in which the current *I* is described by the equation

*S*

_{31}represents the open channel. Because individual channel states are modelled, the current depends only on the amount of

*S*

_{31}and not any of the other closed states. Thus, incorporation into the gating framework described above is intuitive each state represents a structural conformation of the ion channel and the number of channels in each state is explicitly tracked, facilitating a simple approach to account for the energetics of gating under varying ion channel densities.

Figure 5. Channel states of a Na + channel.

##### (v) Voltage dependence of state transitions.

The transition rates between open and closed states are voltage dependent for ion channels. Hodgkin–Huxley models describe state transitions using ODEs of the form

To assess the quality of fit, we compare steady-state open probabilities *g*_{ss} = *α*(*V* )/(*α*(*V* ) + *β*(*V* )) and time constants *τ* = 1/(*α*(*V* ) + *β*(*V* )) (figure 6). The curves for *g*_{ss} and *τ* were generally in agreement however, there were some exceptions. In particular, time constants for the Na + channel gates have lower peaks in the bond graph model when compared to the Luo–Rudy model (figure 6*a*–*c*), but this did not appear to significantly affect Na + channel function as the peaks were all decreased by a similar proportion, facilitating coordination between opening and closing. Similarly, the time constant *τ*_{d} (figure 6*h*) was lower in the bond graph model for some voltages, but given that discrepancies occur at time constants much smaller than the time course of a cardiac action potential we expect that the effects would be negligible. Finally, for the time-dependent K + current *X*_{ss} is substantially higher at negative voltages so that the bond graph model can provide a better match at positive voltages (figure 6*e*). The effects of this difference are partially offset by the lower GHK current at negative voltages which are still above the Nernst potential of K + (figure 4*c*).

Figure 6. Fits for bond graph (BG) parameters against corresponding gating equations from the Luo–Rudy (LRd) model. Steady-state open probabilities are shown on the left panels, and time constants are shown on the right. The bond graph equations are plotted with solid lines, and the Luo and Rudy equations in dashed lines. Gates include (*a*) *m*, sodium activation (*b*) *h*, sodium inactivation (*c*) *j*, slow sodium inactivation (*d*) K1, time-independent K + activation (*e*) X, time-dependent K + activation (*f*) X i , time-dependent K + inactivation (*g*) Kp, plateau K + activation (*h*) *d*, L-type Ca 2+ channel activation. Note that the X i and Kp gates were originally formulated as steady-state equations, thus time constants are shown only for matched bond graph parameters.

#### (e) Finding conserved moieties

Within a biochemical model, conserved moieties are chemical structures that are neither created, removed nor broken down. A common example in energy-dependent metabolic networks is the adenosine moiety found in AMP, ADP and ATP [23,29]. Mass balance specifies that the total amount of each conserved moiety remains constant, and if information on the molecular structure of each species of a reaction network is available, these conservation laws can be derived by counting the number of moieties across all species [29]. In practice, many models do not contain this structural information and this approach cannot be used however, the conservation laws still hold. Here, we outline a method to find conserved moieties using stoichiometric information rather than chemical structures.

Models of cardiac electrophysiology can be represented by the differential equation

#### (f) Stimulus currents

The cardiac action potential model was stimulated using a constant current stimulus that contained enough charge to raise the membrane potential by 30 mV over 0.1 ms. The non-conservative stimulus currents consist of non-specific charge that is not carried by any specific ion. As recommended by Kneller *et al.* [8], conservative stimulus currents contained K + ions as the charge carrier.

### 3. Results

#### (a) Simulation of a single action potential

To verify that our bond graph model reproduced the features of a typical action potential, we simulated the model over a single beat (figure 7*a*–*c*). The membrane potential (figure 7*a*, with stimulation indicated by the arrow) resembled a typical cardiac action potential, with a distinct peak and plateau phase. The contributions of ion channel currents reproduce some common features of cardiac action potentials (figure 7*b*). Once the action potential is initiated by a stimulus current, the sodium current *I*_{Na} briefly activates to give rise to a voltage spike. Following this, the plateau phase occurs where depolarizing L-type Ca 2+ currents oppose the repolarizing K + currents *I*_{K} and *I*_{Kp}. Towards the end of the action potential, *I*_{K1} activates to restore the resting potential [50]. Our model also simulates the reversal of NCX current across the action potential, and the consistent outward current of the Na + /K + ATPase to maintain ionic gradients (figure 7*c*). As a consequence of incorporating the voltage-dependence of gating transitions in a physical framework, transitions between channel states are associated with a gating current resulting from the movement of such charged residues. Our model reveals that the total gating current across all channels *I*_{gate} has minimal contribution to total current (figure 7*c*).

Figure 7. A simulation of the cardiac action potential using a bond graph model. (*a*) Membrane voltage, following stimulation with a conservative stimulus current (arrow) (*b*) ion channel currents (*c*) transporter and gating currents (*d*) membrane voltage over three cycles, for comparison with (*e*) and (*f*) (*e*) power consumption (*f*) energy dissipated, with the variable *E* representing the energy consumption over the duration of the action potential. The model was run initially for 300 ms to allow the membrane potential and channel gates to stabilize. The intracellular ion concentrations were dynamic variables with initial concentrations [Na + _{i}] = 10 mM, [K + _{i} ] = 145 mM and [Ca + _{i}] = 0.12 μM. Constant concentrations were [Na + _{e}] = 140 mM, [K + _{e}] = 5.4 mM, [Ca + _{e} ] = 1.8 mM , [MgATP] = 6.95 mM, [MgADP] = 0.035 mM, [P _{i} ] = 0.3971 mM and pH = 7.095. *T* = 310 K . (Online version in colour.)

Figure 7*e* shows the power consumption of the membrane model over three cardiac cycles which was integrated to estimate the energetic cost of the cardiac action potential (figure 7*f* ). Note that energy continues to be consumed even during the resting state due to the presence of currents associated with ion transporters. Thus, while energy is predominantly consumed during the action potential, there is a positive gradient between action potentials (figure 7*f* ). By setting the energy consumption at the start of the second action potential to zero (figure 7*f* , dotted blue line), we calculated the energetic cost over the duration of the action potential to be 46.8 pJ. As the capacitive area of membrane for this model is 1.534 × 10 −4 cm 2 , the energy consumed per unit membrane area is 305 nJ cm −2 . When compared to Gawthrop *et al.*'s [21] estimate of 173 nJ cm −2 for the energetic cost of an action potential in the giant axon of a squid, the cardiac action potential uses 76% more energy. The main reason for this difference is that in contrast to a neuron, the cardiac action potential contains a plateau phase with opposing currents. Despite the relatively slow rate of change in voltage, the Ca 2+ and K + currents remain relatively high therefore, a large amount of energy is dissipated during the plateau phase.

#### (b) Chemostats influence the conserved moieties of cardiac action potential models

Because the earliest models of the cardiac action potential did not include active transporters, they used constant intracellular concentrations to maintain ionic gradients across multiple cardiac cycles [9,51]. Later models incorporated ion transporters, allowing them to represent physiological conditions with dynamic intracellular ion concentrations, and constant extracellular ion concentrations to model washout from the circulatory system [1,3]. Under ischaemic conditions, washout is greatly inhibited, thus models of ischaemia use dynamic extracellular ion concentrations [4]. We investigated the issue of drift in three classes of model: those with (A) dynamic ion concentrations on both sides of the membrane, representing models of myocytes under ischaemic conditions (B) dynamic intracellular ion concentrations but constant extracellular ion concentrations, representing models of myocytes under physiological conditions and (C) constant ion concentrations, representing models without transporters.

We used our bond graph model to represent these classes of models, selecting ions to fix at constant concentrations that resulted in three variants representative of the classes listed above. Thus, variant A represents models of cardiomyocytes under ischaemic conditions [4], variant B represents models of cardiomyocytes under physiological conditions [1,3] and variant C represents models without ion transporters [9,51]. Conserved moieties of each variant were found using the left nullspace matrix of the stoichiometric matrix (table 1), and these include for example, the total amount of K1 channel (moiety 1). Because the channel is neither synthesized nor degraded in our model, the total amount of channel, i.e. the sum of its closed (C_{K1}) and open (O_{K1}) states, remains constant over the course of a simulation.

Table 1. Conserved moieties associated with chemostat selection. Across some biochemical subgroups (moiety), there are constraints (conserved quantity) on a corresponding sum of species representing the total of the moiety. The conserved quantities remain constant over the course of a simulation. Q represents contributions of other species to charge imbalance across the membrane. The symbol *Σ* represents charge contributions from Markov states of channels and transporters. The definition of *Σ*, and all species can be found in the electronic supplementary material and code.

Similarly, moiety 10 for variant (A) represents the total amount of K + ions, which includes intracellular K + , extracellular K + and the K + ions bound to Na + /K + ATPase. The total amount of K + is constant when ion concentrations are dynamic. However, because fixing the concentration of K + requires an additional external flux, the conservation law is broken in variants (B) and (C). Because the membrane capacitance is included in the stoichiometry of the system, our method automatically identifies a charge conservation law (moiety 13 for variant (A), and moiety 10 for variant (B)).

While it is reassuring that our approach reveals the obvious conserved moieties described above, it also reveals the nontrivial charge conserved moiety that is missed by many existing cardiac electrophysiology models. The overall amount of intracellular charge can be described as a sum of contributions from intracellular K + , Na + , Ca 2+ (and its buffers) and Markov states from ion channels and transporters (*Σ*), similar to forms found in previous studies [10,52]. It should be noted, however, that when all ion concentrations were held constant, charge conservation was broken, as indicated by the absence of a conserved charge moiety in the bottom partition of table 1. In general, holding the concentration of a species constant breaks conservation laws [39] and the number of conserved moieties progressively decreases as more ion concentrations are modelled as chemostats. We discuss the consequences of this in latter sections.

#### (c) Non-conservative stimulus currents cause drift in models with a chargeconservation law

An important feature of cardiac electrophysiology models is that they must be simulated for extended periods to examine physiologically relevant changes in behaviour, thus we tested how the type of stimulus current affected each variant of the cardiac action potential model by pacing at 1 Hz for 30 min. As illustrated (figure 8*a*,*b*), a non-conservative stimulus resulted in drift when the model had dynamic ion concentrations either for all compartments, or only within the intracellular compartment. The drift was particularly pronounced when all ion concentrations were dynamic (figure 8*a*), as extracellular concentrations changed faster than intracellular concentrations. In contrast, the model was resistant to drift from a non-conservative stimulus when all ion concentrations were held constant (figure 8*c*).

Figure 8. Effect of stimulus type and variable ion concentrations on drift in the bond graph model. (*a*) Dynamic ion concentrations (*b*) Dynamic intracellular ion concentrations (*c*) Constant ion concentrations. Results are shown for stimuli that conserve overall charge (blue) and those that do not conserve charge (red). Non-conservative stimulus currents contain non-specific charge, and conservative stimulus currents contain charge and potassium ions. Charge values are given as differences from the initial value of − 5882.2 fmol . T = 310 K . Definitions: V dia , diastolic membrane potential APD, action potential duration at 90% repolarization. (Online version in colour.)

These results suggested that drift arose due to violations of the conserved charge moiety. Charge is a conserved moiety (table 1) in model variants where drift occurred with a non-conservative stimulus. In this situation, non-conservative stimulus currents cause drift because every stimulus causes a stepwise increase in the value of the conserved charge moiety (figure 8*a*,*b* bottom panels). However, because conservation laws are broken as more species are represented as chemostats [39], charge is no longer a conserved moiety when all ion concentrations are constant (table 1). Thus, an observation which may not be obvious to intuition is that under these conditions charge is no longer constant between stimuli, and therefore free to return to its original value after each stimulus (figure 8*c*, bottom panel), allowing such models to achieve a steady-state limit cycle.

Using the observations from the bottom row of figure 8, it is possible to develop a systematic and automated check for drift. Let *v*_{s} be a row matrix representing the stoichiometry of the stimulus current (with chemostats removed), *N* cd be the stoichiometric matrix after removing rows corresponding to chemostats, and *G* be the left nullspace matrix of *N* cd . As seen in this section, a stimulus current will cause drift if it results in any change to the conserved moieties *GX*. Therefore to avoid altering any of the conserved moieties, the stimulus current must have zero contribution to them, i.e. *Gv*_{s} = 0 (or equivalently, *v*_{s} needs to lie in the image of *N* cd ). Thus, the model drifts if *Gv*_{s}≠0.

#### (d) Initial conditions influence steady states through conserved moieties and chemostats

Next, for different sets of conserved moieties (as determined by constrained/dynamic ionic concentrations) we tested how the steady-state behaviour of the cardiac action potential was altered under three different initial conditions (figure 9). The first set of initial conditions (IC1) are common values for comparison (figure 5 at bottom). IC2 is the same as IC1 but with 1mM intracellular K + exchanged for 1 mM of intracellular Na + , such that charge is conserved but K + and Na + are not conserved. Similarly, IC3 is the same as IC1, but with some K + extruded and an equal amount of Na + moved into the cell such that charge, Na + , and K + are all conserved. When all ion concentrations are dynamic, IC1 and IC3 lead to the same steady state, but IC2 results in a different steady state (figure 9*a*). If only intracellular ion concentrations are dynamic, however, IC1 and IC2 result in identical steady states, but IC3 leads to a different steady state (figure 9*b*). Finally, keeping all ion concentrations constant leads to different steady states for all initial conditions (figure 9*c*).

Figure 9. Effect of initial conditions on steady-state behaviour. (*a*) Dynamic ion concentrations (*b*) Dynamic intracellular ion concentrations (*c*) Constant ion concentrations. The models were paced at 1 Hz for 30 min using a conservative stimulus current. [MgATP] = 6.95 mM, [MgADP] = 0.035 mM, [P_{i}] = 0.3971 mM, pH = 7.095, *T* = 310 K. Definitions: *V*_{dia}, diastolic membrane potential APD, action potential duration at 90% repolarization. (Online version in colour.)

These results demonstrate that the summed amount for each conserved moiety and/or chemostat value determines the steady-state behaviour of cardiac action potential models. To investigate this further, we calculated the values for conserved moieties and chemostats that resulted from each initial condition (table 2 *differences from IC1 indicated in italics*). For two sets of initial conditions to achieve identical steady states, all conserved moieties and chemostats must have the same value. Thus, under dynamic ion concentrations (figure 9*a*), IC3 results in the same steady state as IC1 because all conserved moieties have been preserved (table 2), whereas IC2 causes a different steady state because the *K* + and *Na* + conserved moieties take on different values. Similarly, when only intracellular ion concentrations are dynamic, IC2 preserves the value of all conserved moieties and chemostats, but IC3 changes the values of the chemostats corresponding to extracellular Na + and K + concentrations (table 2), hence the different steady state. When all ion concentrations were held constant, changes in the chemostat values (table 2) were associated with different steady states for all three initial conditions (figure 9*c*).

Table 2. The values of conserved moieties and chemostats under different initial conditions. All values are in fmol. Chemostats are indicated with (cs). Values different from IC1 are shown in italics.

### 4. Discussion

In this study, we developed a bond graph model of the cardiac action potential with the aim of resolving the issues of drift and non-unique steady states. Analysis using conserved moieties enabled the discovery of all conservation laws within the model, including the charge conservation law neglected by many existing studies. In addition to the conservation of charge law from previous studies [10,52,53], we found conservation laws corresponding to ions, states of Markov models of channels and transporters, and buffers, demonstrating the comprehensiveness of our approach. Because the bond graph approach requires species and processes to be resolved in biophysical detail, calculation of conservation laws is more straightforward than working purely with the mathematical equations of a model. Two key advantages of our approach over existing analyses are that it reveals all conservation laws in a comprehensive and systematic manner, and that it is general for all models of the cardiac action potential that can be represented as bond graphs. As a result, our approach can be scaled to more complex cardiac models without manual examination of the equations of a model to identify each individual conservation law, and the nontrivial charge conservation law will be accounted for when it appears as a conserved moiety. When simulated over long periods with a non-conservative stimulus, our bond graph model displayed solution drift, but it became resistant to drift when ion concentrations were held constant, demonstrating that changes in the value of a conserved charge moiety drive model drift. We also found that two sets of initial conditions can lead to different steady states if the values of their corresponding conserved moieties and chemostats are different, suggesting a strong link between conserved moieties and the steady-state limit cycles of cardiac action potential models. To demonstrate that our approach is general, we tested how the selection of chemostats (i.e. fixed concentrations) influenced drift and steady states by using variants of our model that were representative of existing models in the literature. Our approach highlights the subtle but relevant observation that holding ion concentrations constant changes the conserved moieties of the model, which in turn has an effect on the susceptibility of a model to drift and non-unique steady states. Because chemostats represent connections between a system and its external environment, they are essential to coupling together biological processes [23]. The coupling of biological processes generally causes changes in the conservation laws of a system which may be difficult to capture through observation. Our approach using bond graphs provides a systematic method of dealing with changes to conservation laws as a result of coupling models together.

#### (a) Drift

When paced with a non-conservative stimulus, variants of the model with a charge conservation law (A and B) underwent drift (figure 8*a*,*b*) consistent with previous studies on the Luo–Rudy model [10,11]. Our observations also explain why models with constant ion concentrations (similar to variant C) are more likely to be resistant to drift [10]. By observing changes in the charge conserved moiety, the bond graph approach attributes drift to regular perturbations in charge that cannot be restored due to the presence of a conservation law. Whereas previous analyses relied solely on intuition to derive a conservation law corresponding to charge [10,11], we note that our approach automatically derives this conservation law, and can also detect other conservation laws that may be relevant for drift.

We note that in order to avoid drift, all conserved moieties (and not only the conserved moiety corresponding to charge) must be preserved by the stimulus current. However, in the examples explored in this study, we found that the stimulus currents preserved the value of all conserved moieties apart from charge, and therefore were not plotted. While intuition may suggest that modelling a stimulus current with K + would cause an accumulation of intracellular K + , the K + is passively transported back through K + channels. Furthermore, the total amount of K + (the relevant conserved moiety) remains constant because the loss of K + from the extracellular side is offset by accumulation on the intracellular side. By contrast, for charge conservation in variants A and B, the model contains no mechanism to reverse the change in charge difference caused by the non-conservative current, thus drift results.

As demonstrated, the bond graph method requires construction of a stoichiometric matrix, providing a simple approach to check whether a stimulus current will cause drift. While it is common practice to use K + as the charge carrier for stimulus currents, it is likely that multiple species contribute to the current [8,10]. Thus, the automated approach suggested here is well suited for checking whether more complex stimulus currents satisfy conservation of charge, as well as other conservation laws within the model. It should be noted, however, that while a model satisfying *Gv*_{s} = 0 will not drift due to violating conservation laws, drift may still occur due to an imbalance of currents throughout the action potential, for instance, in the absence of Na + /K + ATPase, the ionic gradients would gradually disappear in a model with dynamic ion concentrations.

Finally, we believe that this analysis provides a link between the issues of drift and steady states. Our models show that drift due to a non-conservative stimulus current can be attributed to changes in the value of the charge conserved moiety with every stimulus, and accordingly the steady state of the model changes. Model drift then occurs as the solution continually chases a moving steady state.

#### (b) Effects of initial conditions on steady states

We also found that initial conditions of cardiac action potential models change their steady states through the values of chemostats and conserved moieties (figure 9 and table 2). Accordingly, the same perturbation to initial conditions can have different effects on the steady state depending on which species are held constant. Therefore, in addition to ensuring that the concentration of ions is physiological, care should be taken to correctly initialize each state of buffers and Markov models of ion channels and ion transporters, as they may contain a significant fraction of total ion abundance. For example, Ca 2+ buffers and SERCA can sequester a significant amount of Ca 2+ and they should be initialized with the correct amount of bound Ca 2+ when multi-state models are used [54]. We note that the difficulty of manually deriving conservation laws increases exponentially as models of cardiac electrophysiology become more complex, and we believe that our approach extends on existing analyses [10,11] to provide a general method for assessing steady-state behaviour by comparing the values of conserved moieties and chemostats that result from each initial condition.

In the field of biochemical network analysis, there is a well-established dependence of quiescent steady states on conserved moieties, and numerous mathematical techniques for assessing the uniqueness and stability of these steady states have been developed [55,56]. However, the influence of conserved moieties on limit cycles in an oscillating system that is regularly stimulated has yet to be investigated. Our results hint at similarities between these two fields, and while we only tested the uniqueness of steady states using relatively small perturbations to the initial conditions, it is possible that a set of conserved moieties may have multiple steady states, and greater perturbations may lead to other limit cycles.

#### (c) The ‘differential’ and ‘algebraic’ methods

The discovery of conservation principles in cardiac electrophysiology has led to a debate over whether to use the differential or algebraic methods of simulation [7,10,11,52,53]. The differential method is the calculation of membrane voltage by integrating total current, and the algebraic method is the calculation of membrane voltage using an algebraic relationship derived from charge conservation. We chose the differential method over the algebraic method because it better supports model reuse and modularity - in particular it is easier to modify the equations to select different species as chemostats, and to combine equations when two models are coupled. We note, however, that the algebraic method may reduce computational complexity [10,20]. In bond graph modelling, the algebraic method can be implemented by using conserved moieties to turn the system of ODEs into an index-0 differential algebraic equation (DAE) (eqn 3.48 of [20]). This method generalizes existing algebraic methods to reduce the system of differential equations by using all conserved moieties and not just the conserved charge moiety. While we did not use the algebraic approach, we emphasize that the choice of method relates to numerical approaches for model simulation rather than the underlying physics of the system [10]. Therefore, the differential and algebraic methods are equivalent in conservative systems provided that the initial conditions and values of conserved moieties are consistent.

#### (d) Integration into whole-cell models

Our bond graph model of the cardiac action potential is the first step towards a fully integrated whole-cell bond graph model of a cardiomyocyte that couples electrophysiology, signalling, metabolism and mechanics. As energy drives all biological processes, and energy supply can be limited under certain pathophysiological states, it is of interest to examine how cardiomyocytes allocate their energy, and also to estimate the efficiency of processes that are essential for cardiomyocyte function [57,58]. Modelling studies for the energetic regulation of a cardiac cell exist across the literature [58], but while some components used in these models are thermodynamically consistent [59,60], existing whole-cell models are neither energy-based nor thermodynamically consistent throughout the entire model. Furthermore, because existing experimental and modelling studies use ATP consumption as a proxy for energy consumption, they can only estimate the energy consumption of major energy sinks: the Na + /K + ATPase, SERCA and crossbridge cycling [58,61]. A bond graph approach may thus provide more detailed insights into how a cardiac cell uses energy downstream of ATP hydrolysis processes, and help to identify energy-consuming processes. Because the bond graph approach is energy-based, it not only provides the necessary constraints to develop a thermodynamically consistent model, but also allows us to directly assess energy consumption of the model (in Joules). We found that when normalized against membrane area, the cardiac action potential consumes approximately 76% more energy than an action potential in the axon of a giant squid. To the authors' knowledge, this is the first account of energy consumed by electrochemical processes during the cardiac action potential.

#### (e) Limitations

A limitation of our approach is that the physical constraints imposed by the bond graph approach prevent the direct translation of existing models of cardiac electrophysiology into bond graphs. This points at thermodynamic issues that exist within existing models of cardiac electrophysiology while many models contain components that are thermodynamically consistent, the authors are not aware of any models that are entirely thermodynamically consistent. Nonetheless, the inability to perfectly replicate properties of existing models as bond graphs impedes the creation of bond graph models of cardiac electrophysiology, as it requires new models to be built from the bottom up with existing equations only used as guides for parameter fitting. As a result of this limitation, our bond graph model represents only a subset of the ion transport processes within a cardiomyocyte, and only approximates the behaviour described in existing models. While our model is able to reproduce many essential features of the action potential, it is unlikely to be physiologically realistic under conditions different to those used for our fitting process. Thus in future studies, it would be interesting to explore the use of more complex bond graph components to generate more physiologically realistic models. We chose to use a restricted set of components and constitutive equations to keep the model simple, although the bond graph framework is flexible enough to account for a wider range of equations provided that they are thermodynamically consistent. While many ion channels (including those for *I*_{Na}, *I*_{K1} and *I*_{Kp} in the Luo and Rudy model) are modelled using linear I-V relationships, the choice of I-V relationship is generally chosen on the basis of providing an empirical fit to data rather than as a result of physical principles. Thus, a more ideal approach would be to use nonlinear bond graph components derived from physical principles, and fit their parameters directly to experimental data. Additionally, empirical equations for gating transition parameters are generally not expressed in a thermodynamic framework, and are therefore impossible to replicate exactly with a bond graph model. We used a simple gating mechanism to reduce equation complexity while maintaining thermodynamic consistency however, the quality of fit could be improved by using more complicated gating mechanisms, or other thermodynamically consistent constitutive equations.

Because of physical restrictions imposed by the bond graph framework, we were forced to model ion channels and transporters using Markov states to faithfully represent their underlying physics. However, this produced a model that had numerous states compared to the number of biological processes. While it is reassuring to find that our method of identifying conserved moieties remained robust despite this complexity, simulation of the model was computationally expensive. For the purpose of integrating this action potential model into a larger whole-cell model, it would be useful to have simple model components that reduce computational cost. While current methods for reducing biochemical models in the bond graph framework are not advanced enough to apply to the biological components in this study, we note that bond graphs provide a useful foundation for applying model simplification while ensuring that thermodynamic consistency is maintained [20].

We also decided to limit the transport processes included in our model to those considered essential for producing a cardiac action potential, while maintaining a limit cycle using dynamic ion concentrations. Our bond graph model omitted many ionic currents due to their small amplitudes however, these channels may have greater contributions under conditions which vary from those tested here. Thus, an obvious extension of this work would be the integration of other electrogenic processes within the cardiac membrane. It would be interesting to investigate whether coupling other models requires further tuning of parameters [62], and whether the presence of physical bond graph parameters changes this process. A related limitation is that our model does not account for ion concentrations such as *H* + and *Cl* − , as well as pH buffers [5]. While including these ions and their transporters would lead to a more accurate model, the omission of these ions did not cause inconsistencies in the conservation laws described in this study. These ions were assumed to be membrane impermeable in our model, and thus their constant contributions to the membrane potential were accounted for in the initial value of the charge conserved moiety. Similar to how calcium and its buffers were accounted for in our current approach (table 1), our bond graph approach is sufficiently general to account for other ions and their buffers.

When formulating the structure and parameters for a bond graph model of the cardiac action potential (or most other biological processes), it is possible to either fit against existing mathematical models or the underlying experimental measurements. For all processes in this study excluding the NCX, we developed our bond graph model to reproduce the behaviour of an existing model, in an attempt to re-use existing knowledge about these processes. This approach poses constraints on the bond graph structure used, especially for gating structure. Therefore, it would be interesting to develop an approach that assesses bond graph structures as well as bond graph parameters, based on their fits to data [62]. Such an approach may provide a better fit to the data, and uncover insights into the physical mechanisms of ion channels.

### 5. Conclusion

In this study, we have developed a bond graph model of the cardiac action potential and used this to explore the issues of drift and non-unique steady states. We demonstrate that the analysis of conserved moieties generalizes the concept of charge conservation used in earlier studies, and found that changes in conserved moieties can explain drift as well as changes in steady-state behaviour. Importantly, holding ion concentrations constant can have significant consequences on both drift and steady states as they change the conserved moieties in the model. Our approach to resolving drift and non-unique steady states is sufficiently general that it can be applied to any bond graph model of the cardiac action potential model. We hope that the bond graph approach outlined here will prove useful for the development of future cardiac electrophysiology models, and eventually whole-cell models of the cardiomyocyte.

### Data accessibility

The code associated with this study is available from GitHub (https://github.com/uomsystemsbiology/bond_graph_cardiac_AP), and archived on Zenodo (https://doi.org/10.5281/zenodo.1172205) [63]. The repository contains Matlab (The MathWorks, Natick, MA) code that generates the figures, CellML code containing parameters, initial conditions and equations of the model, and full details of the bond graph structure.

### Author's contributions

M.P., P.J.G., J.C. and E.J.C. developed the theory. M.P. performed the research. K.T. provided conceptual advice and helped interpret the results. All authors contributed to the text of the manuscript and gave final approval for publication.

### Competing interests

We have no competing interests.

### Funding

M.P. acknowledges the financial support provided by an Australian Government Research Training Program Scholarship. P.J.G. thanks the Melbourne School of Engineering for its support via a Professorial Fellowship. K.T. is supported by the Heart Foundation of New Zealand (Research Fellowship 1692) and the Marsden Fund Council from Government funding, managed by Royal Society Te Ap a ¯ rangi (Marsden Fast-Start UOA1703. This research was in part conducted and funded by the Australian Research Council Discovery Projects funding scheme (project DP170101358).

Download and print this article for your personal scholarly, research, and educational use.

Buy a single issue of *Science* for just $15 USD.

### Science

Vol 336, Issue 6078

13 April 2012

### Article Tools

Please log in to add an alert for this article.

By Morten Ø. Jensen , Vishwanath Jogini , David W. Borhani , Abba E. Leffler , Ron O. Dror , David E. Shaw

*Science* 13 Apr 2012 : 229-233

Molecular dynamics simulations show how a voltage-gated channel closes.

## Award presentation by Denis Baylor

It’s a pleasure to introduce the winners of this year’s Albert Lasker Basic Medical Research Award: Clay Armstrong, Bertil Hille, and Roderick MacKinnon. Clay, Bertil and Rod have made outstanding discoveries about some of the most important molecules in living cells: ion channels. Ion channels are specialized protein molecules in cell membranes. They generate electrical signals in nerve cells and muscle fibers, produce the cardiac rhythm, and instruct glands to secrete. An individual channel molecule is an on/off switch that controls the flow of ions across the membrane. When switched on, it allows a specific type of ion—Na, K, or Cl—to cross the membrane. These ion movements produce small electrical currents that change the membrane voltage and transmit information from one place to another.

An ion channel is a remarkable molecular machine. Typically it sorts ions one at a time, yet allows the preferred species to pass through at a rate of a million per second, far higher than the speed at which the best enzymes can process their substrates. One of several types of signal can activate channels. Some switched on by the membrane voltage itself have a voltage sensitivity tenfold higher than that of the logic gates in a spanking new computer. Wonderful as they are, channels do not always work properly, and malfunctions produce an expanding list of diseases that includes epilepsy, cardiac arrythmias, and cystic fibrosis.

Before Clay, Bertil and Rod’s work, physiologists knew that ion channels are important, but what channels actually are and how they work were only understood at the level of the formalisms that Hodgkin and Huxley had enunciated two decades earlier. Hodgkin and Huxley’s work, which earned the Nobel Prize in 1963, had revealed that separate pathways allow ions of different types to cross the membrane, and that, for certain neuronal channels, the voltage across the membrane determines whether ions are allowed to pass or not. It was not clear how ions pass through channels, how channels achieve their remarkable selectivity for ions, nor how the voltage across the membrane instructs the channel to allow ion flow or to block it. The complete molecular structure of a channel? Forget it! This of course was before our awardees came on the scene.

#### More Info

Bertil Hille, from the University of Washington, provided a physical mechanism for ion selectivity. He did so by characterizing the gauntlet an ion must run as it traverses the channel’s pore. The approach was beautifully direct: find the size of the hole in a sodium channel by trying a variety of ions and determining the largest that can squeeze through. The bottleneck he measured would allow a sodium ion to pass only after it had shed all but a few of the water molecules that normally surround it. Loss of the waters was energetically compensated by negatively charged carboxyl oxygen atoms in the pore wall. This picture explained how a sodium channel effectively allows only sodium ions to pass, yet lets these “chosen” ions through at a very high rate. Bertil extended the analysis to potassium channels and showed that they achieve selectivity in a more subtle way, stripping potassium ions of water and allowing multiple ions to dwell in the pore at one time. Electrical repulsion among the ions allows them to traverse the pore at very high rates even though they interact strongly with the pore walls. More recently, Bertil discovered that some channels are regulated directly by interactions with G-proteins, ubiquitous signaling molecules of the cell interior. He continues distinguished investigations on fundamental mechanisms of cell signaling. Bertil’s book, *Ionic Channels of Excitable Membranes*, is prized by students and senior scientists alike for its depth and clarity.

Rod MacKinnon, at The Rockefeller University, did the unthinkable: He crystallized a potassium channel and determined its entire molecular structure by X-ray diffraction. Earlier he had used physiology and molecular biology to demonstrate that a potassium channel consists of four subunits and to identify the amino acids that line the channel’s outer vestibule and pore. The next logical step was to visualize the complete structure of a potassium channel. Doing so required him to test countless crystallization conditions, construct many mutants, and in his off hours learn X-ray crystallography. He did all these things in short order, finding a beautiful structure built from four subunits that come together to make an inverted teepee spanning the membrane. The form and makeup of the pore provide a firm structural basis for the high potassium selectivity and ion transport rates in these channels. Indeed, as Hille, Hodgkin and others had inferred, multiple potassium ions reside in the pore. Now the channel that Rod crystallized was from the bacterium *Streptomyces lividans*. In a beautiful sequel to the structural work, Rod showed that the bacterial channel can be blocked by a scorpion toxin that blocks potassium channels in higher organisms. This indicates that the structure of the potassium channel’s mouth was conserved as higher organisms evolved. Rod is extending this work to other channels, and we eagerly await his next structure.

In closing, it is clear that Clay, Bertil and Rod’s fundamental insights into these key molecules will become ever more important in understanding and conquering human disease.

## Discussion

Calcium-activated K + channels are essential players in cell physiology, as they link intracellular signaling with the electrical properties of the membrane. The mechanism of channel modulation by Ca 2+ has been studied in detail for eukaryotic Ca 2+ -activated K + channels (BK and SK channels) but an important piece of the puzzle is still missing: structural information. The structure of MthK was postulated to be a model for the structure of BK channels and a starting point for understanding the mechanism of channel opening upon Ca 2+ modulation. However, little functional information is actually available for MthK (Jiang et al., 2003). Here, we characterize in detail the gating of MthK channels and we describe the mechanism of Ca 2+ modulation in the framework of the available structural information.

### MthK and BK Channels: Functional Comparison

What are the similarities between MthK and BK channels? They are both potassium channels with a homologous pore region and a GYG signature sequence for potassium selectivity. The COOH-terminal cytoplasmic domain of both channels is composed of two RCK domains (identical in MthK, different in BK, but homologous between the two species Jiang et al., 2001), and calcium favors channel opening in both channel types. However, while BK channels are activated by micromolar Ca 2+ concentrations (Latorre, 1989), MthK channels do not respond until Ca 2+ is raised to millimolar concentrations.

One of the most intriguing features of MthK is the very steep activation by Ca 2+ . We examined MthK activation as Ca 2+ was raised from 0 mM (with added EDTA) to 10 mM. We found that the channel activity increased sharply between 3.5 and 5 mM Ca 2+ . This jump in channel activity was well described by a Hill coefficient of at least 8. The value was constrained by the number of Ca 2+ ions present in the high resolution structure of the channel eight Ca 2+ ions were found in the MthK structure, two per functional subunit and one per RCK domain. This is the highest Hill coefficient found in any ligand-gated channel for comparison, the largest Hill coefficient found for the Ca 2+ activation of BK channels is 4–5 (Golowasch et al., 1986 McManus and Magleby, 1991 Cox et al., 1997 Nimigean and Magleby, 1999), for cyclic nucleotide activation of CNG channels it is 2–3 (Zagotta and Siegelbaum, 1996 Ruiz et al., 1999), and for Ca 2+ activation in SK channels it is 2–4 (Kohler et al., 1996 Vergara et al., 1998).

What does this unusually high Hill coefficient mean for the mechanism of Ca 2+ modulation of MthK channels? The simplest answer is that MthK almost never opens unless all eight Ca 2+ binding sites are occupied. This mechanism is different from the Ca 2+ activation of BK channels in which the channel can open in partially Ca 2+ -occupied states and the coupling between Ca 2+ binding and channel opening, although strong, is not as strong as we observe with MthK (McManus and Magleby, 1991 Horrigan and Aldrich, 2002). An alternative to the above simple explanation is that the extreme steepness of the Ca 2+ dose–response is due to multiple overlapping processes that lead to channel opening from closed/dormant/blocked states (i.e., Ca 2+ activation plus one or more additional independent processes that modulate MthK channel gating). Inactivation was recently proposed to be the cause for the apparently low open probability of KcsA channels during steady-state gating (Cordero-Morales et al., 2006). However, it is not obvious that an inactivated state exists for MthK channels (which lack the glutamate, E71, proposed to be involved in KcsA inactivation), and, in the absence of supporting data for a more complex process, we favor the simpler mechanism described above for the origin of the steep Hill slope.

The apparent all-or-none Ca 2+ activation mechanism for MthK is also supported by our kinetic analysis of the open and closed states with Ca 2+ . First, we found that MthK activity has a “bursty” appearance at the single channel level that we catalogued as two gating modes: slow (bursts and gaps between bursts) and fast gating (flickers within a burst). We found that the durations of the open and closed intervals within a burst (the components of fast gating) are Ca 2+ independent and that the increase in Po with Ca 2+ (>100-fold) comes mainly from the marked decrease of the durations of the gaps between bursts (40-fold) and to a lesser extent from the increase in the duration of bursts (two to threefold). Ca 2+ affects only the slow gating of MthK by mainly altering the frequency of opening while the fast gating is Ca 2+ independent. For BK channels, while Ca 2+ similarly increases the frequency of openings, it also markedly increases the durations of the open intervals (McManus and Magleby, 1991).

Therefore, while BK and MthK appear similar on the surface, as they are both activated by Ca 2+ and possess RCK Ca 2+ binding domains, the details of the mechanism of Ca 2+ activation differ between the two channels. This is not surprising, as the structure of the eukaryotic BK channel is likely to be more complex than that of the more primitive MthK. BK channels have an additional five transmembrane domains upstream of the pore region, while MthK has only the two transmembrane regions that form the pore. In addition, the RCK-containing large COOH terminus of BK channels also includes other functional domains as well as other Ca 2+ binding structures/sites (Schreiber and Salkoff, 1997 Shi and Cui, 2001 Zhang et al., 2001 Bao et al., 2002 Xia et al., 2002).

### A Modified Monod-Wyman-Changeux Model Describes MthK Ca 2+ Gating

Using our kinetic data, we derived a qualitative model to describe the Ca 2+ activation of MthK. This model must fulfill the following conditions: (a) the channel opens in 0 Ca 2+ (Fig.1) (b) the fast gating in 0 Ca 2+ is identical to the fast gating in saturating Ca 2+ (Fig. 4, durations of opening and closings within a burst are Ca 2+ independent), and (c) the channel opens more frequently in the presence of Ca 2+ , which increases the open probability ∼100-fold (Fig. 3). The purpose of this modeling effort is not to provide an exact kinetic treatment with precisely determined rate constants, but to show that a minimal allosteric model derived from our experimental findings (with structural constraints) can reasonably reproduce our data.

We tested several allosteric schemes based on the Monod-Wyman-Changeux model (MWC Monod et al., 1965) that were shown to be appropriate for other ligand-gated channels (Cox et al., 1997 Li et al., 1997) and the proposed model in Scheme 1

is a minimal model that satisfies all the above conditions. There must be at least nine closed states and nine open states on the top two rows in order to satisfy the requirement of having a distinct kinetic state for each Ca 2+ occupancy state (open and closed states corresponding to 0, 1, 2, 3, 4, 5, 6, 7, and 8 Ca 2+ ions bound). These rows of states are connected by rate constants that are Ca 2+ dependent and are involved in the slow Ca 2+ -dependent gating of MthK. The Ca 2+ binding equilibrium constant is described by K_{Ca} (K_{Ca} = α[Ca 2+ ]/β, where the forward [α] and backward [β] rate constants represent the intrinsic on and off rates of Ca 2+ binding to both the closed and the open conformation of MthK. For simplicity, we assume that the Ca 2+ binding constants are identical for the open and closed channel conformations (similar results are obtained if we assume they are different). L represents the equilibrium constant (L = k_{o}/k_{c}, where k_{o} and k_{c} represent the opening and closing rates, respectively, which define the vertical transitions between the top two rows). The vertical transitions between the top two rows are identical for each Ca 2+ -occupied conformation, with the exception of the eight Ca 2+ -bound conformation, which opens with a higher probability conferred by a factor f (to satisfy microscopic reversibility, f also multiplies the equilibrium constant of the final transition O_{7Ca} to O_{8Ca}). This maneuver, which distinguishes our proposed model from a classical MWC model ensures that the channel is favored to open mainly when all Ca 2+ binding sites are occupied, as demanded by our results, and at the same time allows a nonzero open probability in the absence of Ca 2+ . The cooperativity coefficient (μ) allows each sequential Ca 2+ binding event to occur with higher probability. In the absence of this cooperativity coefficient, we were unable to simulate a Ca 2+ dose–response curve with Hill coefficients higher than 4. The closed states on the bottom row represent “flicker states.” This equilibrium is described by K_{f} (K_{f} = k_{fo}/k_{fc}, where k_{fo} and k_{fc} are rate constants that do not depend on Ca 2+ and define the fast Ca 2+ -independent gating within the burst).

We were able to constrain most of these rates by examining our results Table II lists the key parameters from our experiments that were used to determine the rate constants (“real” and “simu” denote kinetic parameters extracted from experimental and simulated data, respectively). Since the fast gating behavior is Ca 2+ independent, k_{f0} is determined by the inverse of the mean open time within the bursts (∼1.5 ms at −200 mV) and k_{fc} by the inverse of the mean closed time within the burst (∼0.4 ms at −200 mV). The rate of leaving the bursts (closing rate, k_{c}) was estimated as the inverse of the mean burst duration in 0 Ca 2+ (∼50 ms at −200 mV). The values of these experimentally determined rate constants are in Table I.

Thus, although the model in Scheme 1 may appear complex due to the total number of states presented (27), it is actually the simplest model that will satisfy all the structural and experimental constraints listed above. In addition, due to the fact that we were able to determine directly some of the microscopic rate constants from our kinetic analysis, the 27-state model in Scheme 1 has only five free parameters: the rates of Ca 2+ association/dissociation (α and β), the opening rate k_{o}, the cooperativity factor, μ, and the factor f.

Using Scheme 1 and the constrained parameters, we fit by eye our data to obtain the rest of the parameters (Table I). We then simulated single channel traces (Fig.7, A and B) and a dose–response curve (Fig. 7 C) using the QuB software as detailed in Materials and Methods. The simulated data illustrates visually that Scheme 1 recapitulates the salient Ca 2+ -dependent gating features of MthK (compare simulated single-channel traces in Fig. 7 with the real traces in Figs. 1, 3, and 5). The simulated data was then analyzed to obtain the kinetic parameters (shown together for comparison with the real data in Table II). This comparison shows that using Scheme 1 we were able to capture some detailed characteristics of MthK gating as well. (a) The simulated mean open and closed times within bursts are Ca 2+ independent and identical in values to real data. (b) Simulated burst durations are close in value to real data and change little with Ca 2+ . (c) There was a large change (>100-fold) in the durations of simulated gaps between bursts upon increasing Ca 2+ . Fig. 7 C illustrates that we also captured the steepness of the Ca 2+ dose–response (curve in Fig. 7 C was fit with a Hill coefficient of 7.5) achieved by introducing the cooperativity coefficient, μ. A maximal open probability value <1 (0.3 in the case illustrated in Fig. 7 C) was achieved by both adjusting the open–closed equilibrium L and by varying the factor f, which multiplies the open–closed equilibrium. In Fig. 7 C, with the parameters listed in Table I, the open probability increased ∼100-fold from 0 to 5 mM Ca 2+ , as in our experimental data.

Allosteric models similar to the modified MWC model proposed here for the Ca 2+ gating of MthK have been used to describe Ca 2+ gating in BK channels (McManus and Magleby, 1991 Cox et al., 1997 Rothberg and Magleby, 1998 Rothberg and Magleby, 2000 Horrigan and Aldrich, 2002). While our model is relatively simple, it captures the main features of MthK's Ca 2+ activation.

### Voltage-dependent Gating of MthK

For the classical voltage-gated channels (Shaker, BK channels), increasing membrane depolarization leads to a dramatic increase in the channel's open probability accompanied by an increase in the durations of the open intervals (Bezanilla and Stefani, 1994 Sigworth, 1994 Fedida and Hesketh, 2001). Since this large voltage dependence is caused by a voltage sensor domain contained within the first four to five transmembrane segments, which are absent in MthK, there is no surprise that voltage does not elicit a similar increase in open probability in MthK. At first glance, voltage appeared to only disrupt the flickering behavior of MthK channels, dramatically reducing the frequency of the short closing events within a burst. This has a strong visual impact by lengthening the amount of time the channel is open, while keeping the burst lengths constant. After further kinetic analysis, we found that the open probability also increases with voltage (to a lesser extent) but not because of the intraburst kinetics change. The open probability increased two to threefold for a 150-mV depolarization, due primarily to a twofold decrease in the durations of the gaps between bursts (Fig. 5). Thus it appears that two distinct mechanisms are responsible for the voltage effects observed here with MthK: a major one that decreases the flickering activity within the burst and a minor one that slightly increases the Po by decreasing the durations of the gaps between the burst. We will concentrate on the major voltage effect for the rest of the discussion.

Can the model described above in Scheme 1 be modified to also capture the major effect of voltage on fast gating? To do this, we assigned voltage dependence to the rate constants corresponding to the vertical flicker transitions (middle to bottom row) by assuming that only the closing rates depend on voltage (Eq. 5). We did not assign voltage dependence to the opening rates (k_{fo}), since the durations of the closing events within a burst do not vary with voltage (Fig. 5 E). We used z = 0.62, the value obtained from our analysis of the open duration dependence on voltage in Fig. 5 B. Using the parameters listed in Table I, we were able to simulate single-channel currents that capture the voltage effect on MthK gating (Fig. 7 B). As seen in Table II and Fig. 7 B, we succeeded in precisely capturing the increase in mean open time with depolarization by simply attaching voltage dependence (with a z value of 0.62) to the flicker closing rate. We did not attempt to capture the minor effect of voltage on the open probability using the above model.

### Fast Gate and Slow Gate of MthK

What does this kinetic analysis tell us in terms of mechanism? If we examine the kinetic model proposed here, we notice that there is a clear separation between the states and transitions involved in the Ca 2+ -dependent slow gating and those involved in the voltage-dependent fast gating. This effect is not model dependent, as we showed that Ca 2+ and voltage act on two distinct pathways because of the distinct open and closed interval populations they affect (Figs. 3 and 5).

We propose that there are two distinct gates underlying the gating of MthK. The first is the Ca 2+ -dependent gate proposed previously (Jiang et al., 2002b), located at the inverted teepee at the cytoplasmic end of the inner helices. It was suggested that when Ca 2+ binds, the gating ring formed by eight RCK domains switches conformations and pulls on the linkers that connect the gating ring to the inner helices, mechanically opening the channel. A similar gate was proposed for the BK channels in light of recent findings that different lengths of linkers, at a homologous position to the MthK linker, affect the open–closed equilibrium in a predictable fashion for a “spring”-like model (Niu et al., 2004). This mechanism appears plausible for the slow gate of MthK channels, and our kinetic analysis is indifferent to the actual structural changes that occur upon opening the Ca 2+ gate. However, our analysis puts some constraints on the intricacies of this mechanism. First, the Ca 2+ gate does not open when exposed to Mg 2+ it is a Ca 2+ specific phenomenon. Second, voltage actually affects the slow gate as well (see above, the minor effect of voltage on the open probability). This may occur through a different mechanism than Ca 2+ , and not by modulating the binding of Ca 2+ to open the gate if the voltage also increases the open probability and the duration of gaps between bursts in nominally 0 Ca 2+ (we do not have sufficient data at this time to refute this possibility). Third, our very steep dose–response curve and the low affinity of the binding sites for Ca 2+ suggest that this conformational change that “tugs” on the inner helices to open the channel is favored only when all eight Ca 2+ ions are bound. Due to calcium's low affinity for MthK, it is also likely that the Ca 2+ dissociation rate is very fast (as seen in Table I). If so, then in high Ca 2+ , entry into the all-eight-Ca 2+ -bound open states is favored, albeit predicted to be short lived due to the fast Ca 2+ dissociation rate coupled with a low probability of opening in partially liganded states. However, once an open state is reached, a second gate comes into play, allowing the channel to flicker very fast between open and closed flickery states (bursting, see Scheme 1) during the time the channel has the first gate open.

The second gate is the flicker gate modulated by voltage and responsible for the fast gating behavior of MthK. No indication of such gate is visible in the structure of MthK (Jiang et al., 2002a). There are several potential candidates for an additional gate. One of them is the actual selectivity filter that was proposed to be the only gate for several other K + channels (Flynn and Zagotta, 2001 Bruening-Wright et al., 2002 Kuo et al., 2003). Alternatively, these voltage-dependent flickers could represent open channel block by an external (molecules in the recording solution) or internal (protein structures/side chain residue of the channel proper) charged particle. The only external molecules present are K + , Cl − , Ca 2+ , Tris, and CTX. Tris, Cl − , Ca 2+ , and CTX can be disqualified from the race since the first one is uncharged, the second one is an anion that should not be able to get near the pore, and for the last two, recordings in the absence and presence of either (Figs. 1 and 2) show identical flickering behavior. That leaves K + as a viable cause for the generation of fast gating behavior. Precedents exist for rapid flickery block by permeant ions in CFTR channels (Linsdell and Hanrahan, 1996).

Other possibilities include block by a flexible charged amino acid side chain that is favored to protrude into the pore when the gating ring undergoes a conformational change. We cannot identify the potential culprit side chain due to the fact that the resolution of the channel pore was only high enough to allow backbone tracing of MthK while the amino acid side chains were not resolved (Jiang et al., 2002a).

### Conclusion

MthK is a low affinity Ca 2+ -gated and voltage-gated inwardly rectifying potassium channel. The MthK Ca 2+ dose–response curve has an impressively large Hill coefficient, which suggests that the channel opens mainly when all eight sites are bound by Ca 2+ . MthK has two distinct gating modes: the slow gate, affected by Ca 2+ , and a fast gate (flickers), affected by voltage. The increase in MthK activity with Ca 2+ is due mainly to an increase in the frequency of channel opening. While voltage does not lead to a significant increase in channel activity, it dramatically alters the frequency of flickering closures within a burst. We proposed a modified MWC kinetic model that approximates the main gating features of MthK.