In excitable cells, the molecular motions of ion channels are linked to the whole cell physiology via the transmembrane potential (Vm). To examine such events at the cellular level, we recently proposed a model that links the molecular dynamics (MD) and electrostatics of the slow delayed rectifier K+ channel (IKs) to the cell and tissue electrophysiology (Silva et al., 2009). Although the data for creating such a model are sparse and techniques to do a rigorous simulation are still being developed, the model provides a roadmap for future experiments and simulations that will provide a detailed understanding of the channel role in cardiac electrophysiology.
Here, we outline elements of the model that simulate processes at different time and spatial scales, and highlight recent and potential future advances that could improve the modeling at each level. These elements are:
Structure: Ion channel atomic structure, homology modeling, and MD. Finding energetically stable states of the channel based on atomic structure, and the energy imparted by the Vm.
Gating: Modeling the channel opening and closing in response to changing energetics.
Current: Incorporating the gating transitions of the molecular model into a kinetic scheme that simulates the channel current.
Action potential (AP) and higher levels: Studying the role of molecular interactions in the context of the whole cell AP. Simulating electrocardiographic waveforms that reflect the organ level electrophysiology.
During the past decade, publication of crystal structures for several ion channels has enabled insight into their assembly, permeation, and gating. Several physiologically relevant channels are structurally homologous to the recently published rat Kv1.2 channel (Long et al., 2005a,b). In the absence of a crystal structure, this similarity has facilitated the creation of homology-based models of the IKs channel (Smith et al., 2007) and the bacterial sodium channel NaChBac (DeCaen et al., 2008), among others. Homology modeling is accomplished by aligning the sequence of the channel of interest with a known structure. An initial model can then be created by minimizing an objective function that includes the alignment, van der Waal’s forces, dihedral angles, bond lengths and torsions, and coulombic interactions. This task is accomplished with computer software, for example Modeller (Sali and Blundell, 1993). The alignment we used for KCNQ1 (the α subunit that together with KCNE1 forms IKs) with Kv1.2 serving as a template is shown in Fig. 1 A; it shows good agreement with respect to the positions of the charged residues. Note that four KCNQ1 subunits form a homomeric functioning channel, with KCNE1 modulating its conductance and kinetic properties (Silva and Rudy, 2005).
Once a putative model is generated, it can be tested experimentally for accuracy using a variety of techniques, including metal bridges (Lainé et al., 2003; Webster et al., 2004; Campos et al., 2007; Lewis et al., 2008), disulphide bond formation (Liu et al., 2008; Chung et al., 2009), luminescence resonance energy transfer (Cha et al., 1999; Posson et al., 2005; Richardson et al., 2006), histidine scanning (Cymes et al., 2005), accessibility to methanethiosulphonate reagents (Liu et al., 1997; del Camino et al., 2000; del Camino and Yellen, 2001; Zhen et al., 2005) or fluorescent cysteine labeling (Zheng and Zagotta, 2000; Pathak et al., 2007), and paired amino acid substitutions (Tiwari-Woodruff et al., 1997; Chatelain et al., 2005; Kollewe et al., 2009). Because these experiments provide information regarding the spatial location of residues, they can be used effectively to define ion channel structures. Some examples that have been published recently include refinement of the Shaker channel open state (Lewis et al., 2008), sodium channel structure (DeCaen et al., 2008), the location of the KCNE1 β subunit in the IKs channel (Chung et al., 2009), and the structure of the K2P family of channels (Kollewe et al., 2009).
Additionally, MD simulations can be used to situate the model in a physiological environment, as we did for IKs with explicit lipid and water molecules surrounding the channel (Fig. 1 B). With this model in hand, we focused on the channel voltage-sensing regions, which include transmembrane segments S1–S4, and S5 of the neighboring subunit.
Transition rates between stable channel states are determined by the energy and availability of intermediate conformations. From a computational point of view, these transition rates can be estimated by performing many all-atom MD simulations in parallel for a given time duration. The proportion of these simulations that reach the target state, or a conformation with reasonable similarity, can then be used to estimate the transition rate (Voelz et al., 2010). This approach and others that accurately sample the conformational space pose a difficult challenge for ion channel simulations because of the long timescale of conformational transitions that is usually on the order of several milliseconds (or longer). Currently, long MD simulations typically represent no more than a microsecond in simulated time for a system of this size, and many of these simulations would be needed to probe a transition. It is also not guaranteed that any of the simulations would actually reach the target state. The orders of magnitude difference between possible MD-simulated time and protein transition time, even for fast channel transitions, makes proper sampling with current computational speeds and techniques extremely challenging.
On the horizon are several methodologies that could someday provide a means to alleviate this computational challenge. Most simply, computing speeds have been reliably improving over the last several decades. Also, specialized super computers for MD simulations, for example the Anton machine (Klepeis et al., 2009), have recently become available, allowing all-atom MD simulations on the order of milliseconds. New software approaches, with significant advancement, may also be of help. Some examples are accelerated MD (Hamelberg et al., 2004) and coarse-grained simulations (Scott et al., 2008) that are being developed to minimize the number of computing cycles necessary for long simulations. Finally, methodology has been recently published that uses the string method to find transition pathways in complex systems (Pan et al., 2008). It should be noted that even with substantial progress in these techniques, MD simulations will continue to be a significant limiting step in the modeling process.
In the absence of a full MD trajectory, the transition rates of the molecular motions that determine channel opening and closing can be inferred by calculating a multidimensional energy landscape (Sigg et al., 1999; Sigg and Bezanilla, 2003). This landscape is intended to reflect the stability of the protein as a function of reaction coordinates. Coordinates are chosen to track the progress of the channel between stable states. To describe gating of KCNQ1, we chose translation of S4 along its helical axis and rotation about it (Fig. 1 C). From a macroscopic perspective, possible degrees of freedom include tilting, horizontal translation, and bending of the helix in addition to the motion of other helices. Most simply, an estimate of the energy at a given point can be obtained from the electrostatic energy at each sampled point (Lecar et al., 2003), which is the approach we used by incorporating the voltage-sensing region into an implicit membrane (Fig. 1 D). The energy landscape is then constructed by assigning the electrostatic energy calculated from the Poisson–Boltzmann equation to each point along the reaction coordinates. Because we specified two degrees of freedom (translation and rotation; Fig. 1 C) a reduced (minimal) two-dimensional energy landscape resulted (Fig. 1 E). Future work will be necessary to determine how many degrees of freedom are required to accurately specify the voltage sensor movement.
The movement of positive charges on S4 across the membrane during channel opening and closing causes gating to depend on Vm. The change in the energy landscape, caused by the presence of Vm, can be computed using a modified Poisson–Boltzmann equation (Roux, 1997) (Fig. 1 F). During Vm application, the model correctly predicts favored occupancy in the open state at positive potentials and in the closed state at negative potentials (Fig. 1 G). The gating transitions associated with the movement of the S4 segment due to changes in Vm are then determined by the Smoluchowski equation, which specifies a transition rate reflecting the energy landscape in addition to a diffusive component (Sigg et al., 1999; Sigg and Bezanilla, 2003).
Once gating transitions have been successfully described, a model of a single KCNQ1 channel can be constructed by coupling the molecular transitions associated with the energy landscape to other known movements of the channel, such as inactivation (Silva et al., 2009). For simplicity, these transitions do not necessarily need to be described in as much detail. The main challenge at this level is representing the connectivity between transitions. For example, channel inactivation can be independent of or may occur preferentially from the open state of the channel. The final model (Fig. 2 A), which consists of the energy landscape and its linkage to additional channel states, can then be used to simulate the activity of a single channel via Monte Carlo simulation (Fig. 2 B).
Using the single-channel model as a foundation, the macroscopic current can be obtained by computing the probability of residency in each channel state, or as the sum of activity in many (∼1,000) single channels (Fig. 2 C). Because the current through the channel links it to the other channels in a cell (via the changes the current causes in Vm), the model at this stage is ideal for validation. Experimentally, current kinetics can be readily observed in native cells with drugs that block the channel selectively, or in mammalian and frog oocyte expression systems. Ideally, a molecule-based model will be able to predict the effect of mutations on the current, providing insight into the details of the mutation effects. In our model, we compared several mutations that alter charge at residue E160 on S2 (Fig. 3). The wild-type channel carries a negative charge, glutamate (E). The naturally occurring mutation E160K, that causes the arrhythmogenic long QT syndrome type 1 (LQT1), replaces this negative charge with a fully positive charge. Intermediate to these are neutral residues E160Q and E160A. The model was able to predict slowing of current activation for E160Q and E160A, as observed experimentally (Fig. 3). For E160K, no current was detected in oocytes, so the conductance was set to 0 in the simulation.
AP and beyond
As mentioned above, currents interact with each other at the cellular level through the changes they cause in Vm. This interrelationship is well characterized by the classical circuit model from Hodgkin and Huxley that represents the membrane as a capacitor in parallel with resistors that carry ionic current (Hodgkin and Huxley, 1952). To simulate cardiac ventricular cells, models based on the Luo–Rudy dynamic model (Luo and Rudy, 1994) are often adjusted to match the kinetics of a given species. Due to its similarity to the human AP, the canine ventricular AP has been extensively modeled. One of the most up to date models of the canine is the Hund–Rudy dynamic model (Hund and Rudy, 2004; Decker et al., 2009), which we used to compute the whole cell AP (Fig. 4 A).
The primary benefit of incorporating a molecular model into AP simulations is to understand the role of the modeled kinetic properties in shaping the AP (Clancy and Rudy, 1999). As shown previously for IKs, the most pertinent parameters are its conductance and the rate of its activation relative to KCNQ1 alone (Silva and Rudy, 2005). Fig. 4 B shows that even though the conductance of KCNQ1 is only one fifth that of IKs, KCNQ1 still displays more current initially, due to its faster rate of activation (Fig. 3). The consequence of this interplay between conductance and activation rate can be seen in the difference in AP prolongation when comparing the LQT1 mutation E160K and the LQT5 mutation in KCNE1, L51H. E160K eliminates IKs current completely, whereas L51H prevents the conductance increase from KCNE1 and reduces conductance further by 20% to a total conductance that is 16% that of IKs. Still, because L51H mutant current activates faster, the consequence is only a moderately prolonged AP compared with a much greater prolongation of AP with E160K.
To extrapolate these findings to the body surface, we simulated pseudo-electrocardiographic waveforms. These simulations show QT interval prolongation with a more moderate phenotype for the LQT5 mutant compared with LQT1 (Fig. 4 C). Although this result is to be expected, the ability to simulate the electrocardiogram (ECG) starting from the molecular movement of ion channels shows promise for enabling interpretation of more complex waveforms for clinical diagnosis (Gima and Rudy, 2002). Akin to the accurate modeling of molecular transitions, the simulation of an accurate body-surface ECG will require substantial effort. Some current challenges are modeling accurately the spatial heterogeneity of cardiac APs in different regions of the heart, representing the distribution of non-excitable cells and gap junctions and the connectivity of cells, and the computational challenge of simulating the electrical activity of 5 × 109 myocytes (Adler and Costabel, 1975).
By formulating a comprehensive multi-scale model, we put forward testable hypotheses. For example, we assume that the folding of the KCNQ1 channel protein and the motion of its voltage sensor are similar to those of KV1.2, but the amount of charge moved per channel is less due to a missing charged residue. We also propose that the interaction with the KCNE1 β subunit slows the gating transitions, and this slowing counters an increase in the macroscopic conductance conferred by KCNE1. At the level of electrocardiographic waveforms, we hypothesize that the removal of negative charge at residue 160 on S2 will prolong the electrocardiographic QT interval. The model simulations can predict the physiological consequences of these hypotheses and propose experiments to test them. It is also hoped that by linking specific molecular changes to clinical phenotypes, the modeling approach will assist with individualized patient treatment.
We thank A. Kollewe for helpful comments on this manuscript.
This work was funded by National Institutes of Health National Heart Lung and Blood Institute Merit Award R37-HL33343, grant R01-HL49054 (to Y. Rudy), and fellowship F31-HL68318 (to J. Silva). Y. Rudy is the Fred Saigh Distinguished Professor at Washington University.