Abstract
In-silico simulation based on Markov chains is a powerful way to describe and predict the activity of many transport proteins including ion channels. However, modeling and simulation using realistic models of voltage- or ligand-gated ion channels exposed to a wide range of experimental conditions require building complex kinetic schemes and solving complicated differential equations. To circumvent these problems, we developed IonChannelLab a software tool that includes a user-friendly Graphical User Interface and a simulation library. This program supports channels with Ohmic or Goldman-Hodgkin-Katz behavior and can simulate the time-course of ionic and gating currents, single channel behavior and steady-state conditions. The program allows the simulation of experiments where voltage, ligand and ionic concentration are varied independently or simultaneously.
Key words: software, ion channels, simulation, kinetic models, Markov, IonChannelLab
Introduction
Ion channels are membrane-spanning proteins that control the flow of ions across the membranes of living cells.1 These proteins are essential to the generation of electrical signals in excitable cells and, therefore, are of fundamental importance for the functions of the brain, heart and muscle.1–3 Gating is the process that controls the opening and closing of the ion conducting pore. During gating, an ion channel may visit many conformational states while undergoing transitions controlled by transmembrane voltage or ligand binding. Discrete conformational states and the transitions that connect them are commonly represented by Markov chain models,1–4 which help understand the molecular conformational changes that occur during ion channel gating.5–9 General-purpose software packages, such as MATLAB and Mathematica, can be used to simulate gating kinetics of ion channels. Also, more specialized simulation programs, such as QUB Express, ChannelLab, IChSim, StateEditor and PulseSim, are available to perform this task. The main advantage of general-purpose programs is their flexibility; the user can perform almost any kind of simulation thanks to a large number of algorithms available in these packages. The main disadvantage is that implementation of the model and stimulus is in general time consuming and sometimes requires a high-level of programming experience. In contrast, specialized programs for kinetic simulation of ion channels are more popular because researchers can place more attention on the experimental work and analysis rather than model implementation. Because we consider that a successful program for kinetic simulation needs to be user-friendly and flexible enough for general applications across different fields we developed IonChannelLab. This program includes a Graphical User Interface (GUI), intuitive model-building dialogues and tools and a simulation library. IonChannelLab can be used to build kinetic models, stimulus protocols and perform simulations over a large range of experimental conditions.
Implementation
IonChannelLab uses Markov chain models to simulate the dynamic behavior of voltage- and ligand-gated ion channels. Accordingly, a kinetic model has states that represent discrete conformations and/or ligand occupancies of the channel and rate constants that govern transitions between states. Figure 1 provides a schematic overview of the IonChannelLab architecture. As mentioned before, the software includes a GUI and a simulation library. The GUI can be used to build kinetic models and experiments with voltage and ligand dependences. The simulation library includes algorithms to simulate time-courses, steady-state conditions and time constants. The time-course and steady-state simulations include: ionic current, gating current, conductance and custom functions of state probabilities. The simulations can be obtained using numerical integration, spectral decomposition of the transition rate matrix (Q-Matrix) and Monte Carlo.
Figure 1.
Simulation flow with IonChannelLab. User input properties and program processing are represented by rectangles, respectively. Typically, an IonChannelLab project contains an arbitrary number of experiments as well as the kinetic model. Projects can be saved as .ichl files. The simulated data and stimulus can be saved as .txt or .cvs files. From the kinetic model and the experimental conditions, the program builds a system of differential equations which describe the dynamics of the channel. The time course simulation can be obtained using numerical integration, spectral decomposition of the transition rate matrix (Q-Matrix) and the Monte Carlo method.
Ion Currents
Voltage-clamp experiments record either macroscopic currents resulting from an ensemble of ion channels or unitary currents through a single ion channel. Even though single channel current provides more detailed kinetic information, macroscopic current recordings are widely performed and used to build conceptual kinetic models. Thus, we describe a framework to study macroscopic recording in terms of kinetic models.
If the current through an ion channel in the jth state is represented by ij. The current through NC ion channels can be written as:
| Eqn. 1 |
where NS is the number of conformational states and pj is the probability of occupation of the j state. The single-channel current ij is commonly described using Ohm's law.1 For ion channels obeying this law, the single channel current is given by:
| Eqn. 2 |
where M is the number of permeant ions, V is the voltage, gS is the unitary conductance for the ion S, R is the gas constant, zs is the ion valence, F is the Faraday constant, T is the temperature (in K) and [S]o and [S]i are the extracellular and intracellular molar concentrations of permeant ions, respectively. However, for some open channels the current-voltage relationship is not linear, particularly when ion gradients are imposed, and therefore it is not well described by Equation 2. To take into account ionic gradients the Goldman-Hodgkin-Katz equation (GHK) is used to describe the open channel I–V curve:
| Eqn. 3 |
where PS is the permeability of ion S.
The kinetic behavior of an ion channel with NS states can be obtained by calculating the probability of visiting these states. This probability is determined by solving a system of differentials equations. If the probability of occupying the ith state is represented by pi, then pi is given by:
| Eqn. 4 |
where kij and kji are rate constants that describe the transition from state i to state j and vice versa. The set of differential equations for the probabilities can be expressed in terms of matrix notation as:
| Eqn. 5 |
where P(t) is the probability vector and Q is the transition rate matrix, the elements of the matrix Q are specified by
When voltage or ligand concentrations are changed in a stepwise manner (the most common protocol in voltage- and concentration-clamp experiments), the transition rate matrix is time-independent. Then, the probability of occupying the ith state can be written as:
| Eqn. 6 |
where the time constant τj=1/λj; λj is an Eigen value of the matrix −Q.4
Gating Currents
Gating current results from displacing charged parts of a voltage-gated ion channel within the membrane electric field but without transporting these charges across the membrane.10 If qG,i is the charge associated with the state i and qG is a vector with elements qG,i, then the total charge (QG) can be written as:
| Eqn. 7 |
The gating current is Ig(t) = dQg/dt. Thus, by taking the derivative of Equation 7 and using Equation 5, the gating currents can be written as:
| Eqn. 8 |
Results and Discussion
Implementation of a kinetic model using ionchannellab.
Figure 2A displays a hypothetical Ca2+-dependent ion channel in three different states: closed (C), closed with Ca2+ bound (CCa) and open with Ca2+ bound (OCa). A plausible kinetic model for this channel can be represented by the three-state model shown in Figure 2B. Because the Ca2+ binding site (green sphere) is located internally in the channel and the transition from the state CCa to OCa involves a movement of Ca2+, then the rate constants k+(V), k2(V), a(V) and b(V) are assumed to depend exponentially on membrane potential. A screenshot of the kinetic model editor of IonChannelLab that illustrates implementation of such kinetic model is shown in Figure 3A. In this editor, the rate constants for the transitions between states C and CCa are represented by k1 and k2. As illustrated in Figure 2, this transition is triggered by Ca2+, therefore, k1 = k+*[Ca2+]i. The rate constant k1 is defined in the rate constant editor (Fig. 3C). This editor includes a set of predefined voltage and ligand-dependent functions that can be modified to generate custom rate constants. Since k1 is Ca- and V-dependent, the rate constant k1 is represented by the voltage and ligand-dependent function F = P0[X]e−P1VF/RT, where the ligand concentration is represented for [X]. The generic function F is then adapted to the rate constant k1 by selecting internal Ca2+ as the ligand and specifying the values for Po and P1. The permeation model, permeant ion(s) and the single channel conductance are implemented using the permeation editor (Fig. 3D). In this example, the permeant ion is Cl− and the single-channel conductance is 1 pS. The implementation of complex kinetic models using IonChannelLab is illustrated in Figure 3B.
Figure 2.
Hypothetical Cl− channel activated by Ca2 and voltage. Scheme of a cell membrane containing three hypothetic ion channels in three different states: closed (C), closed with Ca2+ (CCa) and open with Ca2+ (OCa). The ion channel is permeable to Cl−. (B) Representation of the gating transitions using a sequential kinetic model. The kinetic model assumes that the ion channel requires binding of Ca2+ in order to open and that the rate constants k+(V), k2(V), a(V) and b(V) depend exponentially on membrane potential.
Figure 3.
Implementation of a kinetic model using IonChannelLab. (A) Implementation of the kinetic model displayed in the Figure 2B in IonChannelLab. Each state is labelled as in Figure 2 and interconnected by rate constants. The permeant ion is indicated on the upper left side. On the bottom left side, rate constants editor is shown. (B) Screenshots that illustrate the implementation of complex models, the upper and lower parts show a kinetic model of 16 and 13 states, respectively. (C) Screenshot of the rate constant editor with the parameters used for the rate constant k1 (=k+*[Ca2+]i). The rate constant editor includes a set of predefined voltage- and ligand-dependent functions that the user can modify to define the rate constants. (D) The permeation model editor is used to define the details of the permeation model, permeant ion(s) and the single channel conductance. In this example, the permeant ion is Cl− and the single-channel conductance is 1 pS.
The variable editor is also used with the example shown in Figure 2, where the open probability of the channel is modulated by both voltage and the internal concentration of Ca2+, and the current in the open state is dependent on the internal and external concentration of Cl−. In IonChannelLab these experimental variables can be controlled using an unlimited number of steps and/or ramps. For example, Figure 4A and B display the voltage and [Ca2+]i protocols used to simulate the currents shown in the Figure 5A. In each case the number and duration of sweeps can be indicated. Complex waveforms of different amplitude and duration for voltage and ion concentrations can be constructed in a similar fashion as constructing voltage clamp protocols when using the software pClamp. As detailed in Figure 4A, the voltage is changed from −100 mV to 100 mV using square 500 ms pulses from a holding voltage of −100 mV, then the voltage is returned to −100 mV during 400 ms and finally changed to +100 mV for the remainder of the time (1,000 ms). Figure 4B illustrates a basal [Ca]i of 50 nM during the first 1,000 ms of the experiment, then the concentration changes from 0 to 100 nM during 500 ms and then it is clamped at 0 nM during the last 500 ms. The holding values used for voltage and [Ca2+]i were −100 mV and 50 nM, respectively, whereas the internal and external Cl− concentration were held at 100 mM.
Figure 4.
Variables Editor. (A) Screenshot of voltage changes used to simulate the ionic current displayed in Figure 5A. In the holding values section, the experimentalist indicates the holding Vm and ion concentrations. In the variable section, we can indicate the variable that is changing, in this case Vm. We can also define the number of sweeps, duration of each steps and Vm value. In the illustrated example, the voltage was changed from −100 mV to 100 mV in 10 mV increments during 500 ms. (B) Screenshot of the internal Ca2+ concentration used during the simulation of the ionic current displayed in Figure 5A. The value of the [Ca2+]i was changed between 0 and 100 nM.
Figure 5.
Ionic current simulation. (A) Screenshot of the simulation tab of IonChannelLab, the figure shows the ionic current simulation using the kinetic model displayed in Figure 2B, when the voltage and the internal Ca2+ concentration were changed as indicated in Figure 4A and B, respectively. (B) Single channel simulation using the kinetic model displayed in Figure 2B. The membrane voltage was changed from −100 mV to +100 mV, the internal Ca2+ concentration was 50 nM. This simulation includes the addition of white noise. (C) Simulation of the steady-state conductance (upper parts) and time constants (bottom parts) as a function of voltage and [Ca2+]i,. (D) Dynamics of state occupancy is visually displayed along with numerical values representing the percentage of occupation for each state. The upper part shows the state occupation for the scheme shown in Figure 2B when the simulation time is 300 ms. The percentage of occupation was 49, 8 and 42 for the states C, CCa and OCa, respectively. The bottom left part displays the voltage protocol used to active the channel and the bottom right part shows the corresponding whole cell current.
Simulations using IonChannelLab.
Another feature of IonChannelLab is the simulation editor that allows the user to define the type of experiment to be simulated (time course of macroscopic ionic current, gating currents, conductance, time constants, open probability or single channel current) as well as to execute the simulation and obtain a modeled response. This editor permits the creation of an unlimited number of experiments in the same project. In the experiment part the user can define the simulation method, number of channels, sampling rate, number of points to simulate each current trace, number of sweeps, temperature, noise, etc. Figure 5A displays the ionic current simulated according to the kinetic model showed in Figure 2B, when the voltage and [Ca2+]i are changed as indicated in the stimulus part (Fig. 4A and B). Figure 5B shows the single channel current simulated using the voltage square pulse (from −100 to +100 and back to +50 mV during 2.9s) shown above the simulated current. Single channel current is easily seen due to long open time duration even though the current amplitude is only 0.1 pA. Figure 5C (upper left part) displays the steady-state conductance as a function of voltage at different [Ca2+]i, (0 to 200 nM). As [Ca2+]i increases the sigmoideal curves are shifted to the left on the voltage axis since at a single Vm, the conductance increased as the [Ca2+]i augmented. Our hypothetical channel (Fig. 2A) binds Ca2+ in a Vm-dependent manner. This property is illustrated in Figure 5C (upper right part) where the steady-state whole cell conductance is plotted as a function of [Ca2+]i (0 to 200 nM) at various voltage values ranging from −100 mV (small conductance) to +100 (larger conductance) mV in 20 mV steps. At the same [Ca2+]i, positive Vm facilitated channel opening, which is manifested as an increase in conductance. Figure 5C also display the time constants (Eqn. 6) as a function of voltage when the [Ca2+]i = 100 nM (lower left part), or as a function of [Ca2+]i (lower right part) when the voltage is +100 mV. Finally, an important feature of IonChannelLab is the capability to visualize the dynamics of state occupancy. This is particularly important for teaching purposes since the user can visually judge how the channel visits each state during the transition from the closed to the open conformation. Figure 5D shows the dynamics of state occupancy (for the model in Fig. 2B) both numerically and schematically—color-coded—when a 300 ms depolarization to +100 mV is applied from a holding voltage of −100 mV.
Numerical Methods
In IonChannelLab, the numerical calculations are based on the numerical library DotNumerics (www.dotnumerics.com). The solutions of Equation 5 can be calculated in IonChannelLab by numerical integration or by a spectral decomposition of the transition rate matrix Q.4 IonChannelLab includes two methods of numerical integration: the implicit Runge-Kutta method of order 5 with step size control and the Gear method or Backward Differentiation Formula.11,12
Conclusions
IonChannelLab is a flexible and user-friendly software for ion channel simulation that can be readily applied to research and education. It includes intuitive GUI that allows building simple or complex kinetic schemes and easy manipulation of model parameters and experimental conditions. These kinetic schemes are readily adapted to study either ligand- or voltage-gated ion channels both at the single and macroscopic levels, as well as gating currents for voltage-gated ion channels. In the research laboratory, IonChannelLab could be used to suggest plausible gating mechanisms after an easy evaluation of many kinetic models. In addition, this software can be used to design experiments to test the suggested gating mechanism. IonChannelLab can be easily implemented in the classroom to teach basic principles of kinetic theory, advanced topics of kinetic modeling and the fundaments of ion channel gating. The classical Hodgkin-Huxley and contemporary kinetic models for Na+ and K+ channels are already built into this software and can be used to demonstrate the mechanisms of voltage-dependent activation and inactivation, processes that are vital for the initiation and propagation of action potentials in nerve cells.
Acknowledgements
We thank Macrina Calixto for her valuable suggestions during the realization of this project and critical reading of the manuscript.
This work was supported by a Postdoctoral Fellowship Award from CONACyT Mexico to J.A.De S.-C. and grants from CONACyT, Mexico (79897 to J.A. and 105457 to P.P.C.) and the National Institutes of Health (R01 NS032337 to M.C.).
Abbreviations
- GUI
graphical user interface
- Vm
membrane potential
Footnotes
Previously published online: www.landesbioscience.com/journals/channels/article/13404
Availability and Requirements
Availability: The software is freely downloadable from:www.jadesantiago.com/Electrophysiology/IonChannelLab/
Demo videos: http://jadesantiago.com/Electrophysiology/IonChannelLab/Videos/
Operating system: Windows XP SP2, Windows Vista and Windows 7.
Dependencies: Microsoft .NET Framework 2.0 (or higher).
Programming Language: The core algorithm and the GUI were written in C#.
References
- 1.Hille B. Ionic Channels of Excitable Membranes. Second Edition. Sunderland MA: Sinauer Associates Inc; 1992. Classical biophysics of the squid giant axon; pp. 23–58. [Google Scholar]
- 2.Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117:500–544. doi: 10.1113/jphysiol.1952.sp004764. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 3.Luo CH, Rudy Y. A dynamic model of the cardiac ventricular action potential: I. simulations of ionic currents and concentration changes. Circ Res. 1994;74:107–196. doi: 10.1161/01.res.74.6.1071. [DOI] [PubMed] [Google Scholar]
- 4.Colquhoun D, Hawkes AG. A Q-matrix cookbook: how to write only one program to calculate the single-channel and macroscopic predictions for any kinetic mechanism. In: Sakmann B, Neher E, editors. Single-Channel Recording. Second Edition. New York: Springer; 1995. pp. 589–633. [Google Scholar]
- 5.Zagotta WN, Hoshi T, Aldrich RW. Shaker potassium channel gating. III: evaluation of kinetic models for activation. J Gen Physiol. 1994;103:321–362. doi: 10.1085/jgp.103.2.321. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 6.De Santiago JA, Nehrke K, Arreola J. Quantitative analysis of the voltage-dependent gating of mouse parotid ClC-2 chloride channel. J Gen Physiol. 2005;126:591–603. doi: 10.1085/jgp.200509310. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 7.Dougherty K, De Santiago-Castillo JA, Covarrubias M. Gating charge immobilization in KV4.2 channels: the basis of closed-state inactivation. J Gen Physiol. 2008;131:257–273. doi: 10.1085/jgp.200709938. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 8.Tang QY, Zeng XH, Lingle CJ. Closed-channel block of BK potassium channels by bbTBA requires partial activation. J Gen Physiol. 2009;134:409–436. doi: 10.1085/jgp.200910251. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 9.González-Pérez V, Neely A, Tapia C, González-Gutiérrez G, Contreras G, Orio P, et al. Slow inactivation in Shaker K channels is delayed by intracellular tetraethylammonium. J Gen Physiol. 2008;132:633–650. doi: 10.1085/jgp.200810057. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 10.Bezanilla F. How membrane proteins sense voltage. Nat Rev Mol Cell Biol. 2008;9:323–332. doi: 10.1038/nrm2376. [DOI] [PubMed] [Google Scholar]
- 11.Hairer E, Wanner G. Springer Series in Computational Mathematics. Second Edition. Berlin: Springer-Verlag; 1996. Solving Ordinary Differential Equations II: Stiff and differential-algebraic problems. [Google Scholar]
- 12.Hindmarsh AC. Odepack, a systematized collection of ode solvers. In: Stepleman RS, editor. Scientific Computing. Amsterdam: North-Holland Company; 1983. pp. 55–64. [Google Scholar]