MATLAB-based Simulation of Whole-Cell and Single-Channel Currents.

Mathematical models of electrophysiological data are used to investigate biophysical mechanisms that underlie electrical excitability. Although the resources and time required for obtaining experimental data to create these models may not be available to undergraduates enrolled in a biophysics course, computational tools that simulate cellular or single-channel responses to electrophysiological stimuli can be utilized to provide these data. We have developed two MATLAB-based simulation packages that are being used in a cellular electrophysiology course for upper-level undergraduate engineering students to demonstrate the design of electrophysiological stimuli, and the analysis and modeling of ionic currents in excitable tissues. The first package simulates a Hodgkin-Huxley style voltage-gated current elicited during voltage-clamp experiments. Users specify the duration and magnitude of a voltage waveform; the model returns a simulated whole-cell current traces with superimposed noise, and various measurements including peak current, steady state current, and time constants from exponential fits of the current time course. The second package simulates a voltage- or ligand-gated single-channel current as a stochastic process using a state transition matrix. Users specify the membrane voltage, ligand concentration, and number of trials; the model returns simulated single-channel current traces with superimposed noise, and various measurements including amplitude and dwell time histograms. This software has been used during lectures to demonstrate various principles in class, and for class projects in which students derive kinetic models that underlie currents obtained during whole-cell and single-channel recordings. These software packages are freely available and can be downloaded at www.eng.utoledo.edu/∼smolitor/download.htm.

The analysis of data generated by biophysical experiments is an important component for teaching the principles that underlie the electrophysiological responses of excitable tissues. However, many factors such as available facilities and resources, class size, or time constraints prevent such experiments from being feasible. Computational simulations provide a valuable surrogate for these experiments and can expose students to the design of biophysical experiments, the acquisition and analysis of data from these experiments, and the development of mathematical models of the underlying biophysical phenomena.

Existing simulation programs such as NEURON (Hines and Carnevale, 1997; Hines and Carnevale, 2001) and QuB (Qin et al., 1996; Qin et al., 2000a; Qin et al., 2000b) provide a powerful platform for simulating whole-cell or single-channel currents in response to electrophysiological stimuli and pharmacological stimuli. Although these and other packages are available without cost, they may not be part of the standardized computer desktop available to students in engineering computing facilities. In addition, these packages may require substantial training to use, or may not allow the instructor to hide biophysical parameters that would ordinarily be unknown to an experimenter during whole-cell and single-channel recordings.

Other packages that demonstrate how ionic currents contribute to electrical excitability, such as the NERVE program (http://nerve.bsd.uchicago.edu/) or HHsim (www.cs.cmu.edu/∼dst/HHsim/), are freely available and are easy to install and use. These packages can be used to examine the role of individual conductances in action potential generation, and are very useful for introducing students to the biophysical mechanisms that underlie electrophysiological responses. However, we feel it is important for our students not only to use computational models to understand underlying mechanisms, but to learn how these computational models are created from experimental data. To our knowledge, no tools exist that simulate electrophysiologic responses with the objective of creating computational models of the underlying biophysical mechanisms.

MATLAB (MathWorks, Inc.) is a commonly used software package in many engineering courses, and is generally available to engineering students at on-campus computing facilities. MATLAB combines powerful computational processing capabilities with an extensive programming language similar to C++. The MATLAB programming environment also provides for the creation of graphical user interfaces (GUIs) for user-friendly programs. In addition, programmers can create new MATLAB commands through the use of compiled FORTRAN or C++ programs known as MEX files. Unlike standard MATLAB programs that are text files containing source code any user can view, MEX files are compiled binary files in which the original source code remains hidden from the user unless otherwise provided.

To help engineering students understand the process of modeling biophysical phenomena, we have developed two MATLAB-based software packages that simulate recordings of whole-cell and single-channel currents. Students utilize a MATLAB-based GUI to design experiment stimuli, to simulate whole-cell or single-channel currents evoked by these stimuli, to measure parameters from these simulated currents, and to create biophysical models of the underlying phenomena. To prevent students from knowing model parameters a priori, the portion of the code that simulates whole-cell or single-channel currents are contained within compiled MEX files whose contents remain hidden. We have successfully used these software packages in an introductory electrophysiology and biophysics course for upper-level undergraduate students in a bioengineering program. This software, developed to run using MATLAB 6.5 or higher under Windows 2000 or XP, is freely available and can be downloaded as a ZIP file from www.eng.utoledo.edu/∼smolitor/download.htm.

WHOLE-CELL CURRENTS

To simulate ionic currents obtained from whole-cell voltage-clamp recordings, a MEX file was written to implement a computational model of voltage-gated currents. This file contains C++ code to simulate a non-inactivating and an inactivating voltage-gated ionic current based on the Hodgkin-Huxley model of ionic currents from the squid giant axon (Hodgkin and Huxley, 1952). The equations for non-inactivating and inactivating currents are of the form: The activation gates n and m and inactivation gate h vary with time and voltage; these variables represent the probability that an ion channel would be permissive to ion flow at any given instant of time. The values of these activation gates as functions of time are calculated by the following differential equations: Dependence upon the membrane potential Vm is imparted by the steady-state parameters n∞(Vm), m∞(Vm) and h∞(Vm) and the time constants τn(Vm), τm(Vm) and τh(Vm). The equations for these parameters are easily accessed in the code and can be modified by the instructor. Constants for the maximal conductance g̅max, the reversal potential Erev and the number of activation gates N and M are also accessible and can be easily changed. The instructor can also specify the amount of Gaussian noise that is added to the simulated currents. A text file accompanies the software that provides detailed instructions for changing model parameters.

Once the model was implemented and compiled as a MEX file, a GUI was constructed using the MATLAB programming environment to provide a straightforward interface for students to simulate voltage-gated currents (Figure 1). The main graph window displays simulated current traces (black lines) and corresponding curve fits (red lines) with voltage traces below; a bounded Marquardt-Levenberg minimization routine (Press et al., 1992) was ported to MATLAB to obtain model parameters using curve fits. Smaller graph windows along the bottom display peak current amplitudes, steady-state current amplitudes, activation time constants ( ) and inactivation time constants ( ) from curve fits of the simulated current traces as a function of the magnitude or duration of command voltage steps.

Figure 1:

MATLAB-based graphical user interface (GUI) to simulate whole-cell recording using the Hodgkin-Huxley model of the squid giant axon Na+ current. Current traces (black lines), curve fits (red lines), and command voltage waveforms are displayed on the upper left; parameters from curve fits including peak current, steady-state current, and time constants are displayed along the bottom. The user can specify various parameters for command voltage waveforms and analysis of simulated currents using edit boxes and buttons displayed on the upper right; see text for more details. The simulation parameters were specified to evoke Na+ currents from a holding potential of −95 mV using a 6 ms command step that varies from −95 mV to −5 mV in 10 mV increments. The analysis parameters were specified to fit current traces during the 6 ms command step assuming the presence of three activation gates.

To introduce concepts required for the design of voltage-clamp experiments, students are required to select various parameters for the acquisition of whole-cell currents (Figure 1). As with real experiments, the duration (Tstep) and magnitude (Vstep) of voltage steps in the command voltage waveform can adversely affect the results from the fit algorithm. Voltage steps with durations that are much shorter or much longer than the activation time constants may produce erroneous results: extremely short voltage steps may not elicit enough activating current to obtain an activation time constant; and extremely long voltage steps may not provide enough data points sampled during the activation time course before the steady-state current is reached. A control button is provided for users to simulate a single current trace (Simulate One) with or without a curve fit (Curve Fit = yes or no) during experiment setup; students can also design voltage-clamp protocols to obtain a sequence of simulated currents for simultaneous analysis (Simulate Seq). Parameters for these protocols include whether to vary the duration or magnitude of a step during a sequence (Vary Type); which voltage step to vary during a sequence (Vary Step); and the values of the magnitude or duration during the sequence in the form min : step : max (Sequence). These parameters can be modified to design a tail current protocol to obtain the reversal potential and deactivation time constants of a non-inactivating current (Figure 2), or a two-pulse protocol to examine the recovery of an inactivating current (Figure 3).

Figure 2:

Tail current protocol to obtain the deactivation time constants and reversal potential for the Hodgkin-Huxley model of the squid giant axon K+ current. K+ currents are activated by an 8 ms command step to +20 mV from a holding potential of −70 mV; deactivating tail currents are subsequently evoked by a 6 ms tail step that varies from −100 mV to −50 mV in 5 mV increments. The reversal potential and K+ conductance at +20 mV is obtained from the linear peak current vs. tail voltage graph on the lower left; deactivation time constants as a function of tail voltage are obtained from the graph on the lower right.

Figure 3.

Two-pulse protocol to measure the inactivation recovery time constant for the Hodgkin-Huxley model of the squid giant axon Na+ current. Na+ currents are activated and fully inactivated by a 6.5 ms command step to +20 mV from a holding potential of −120 mV. A second command step with the same magnitude and duration is used to measure the recovery of inactivated Na+ currents following a repolarizing step to −90 mV whose duration varies from 0.5 ms to 5 ms in 0.5 ms increments, and 6 ms to 10 ms in 1 ms increments. The recovery time course is obtained from the exponential peak current vs. repolarizing step duration graph on the lower left.

Students are required to specify additional parameters for the analysis of simulated currents. Initial simulations are performed to determine the minimal number of activation gates (Act Gate) needed to obtain an accurate fit of simulated currents. In addition, the user must specify which voltage step to obtain curve fits of simulated currents (Measure Step). For example, obtaining the inactivation recovery time course (Figure 4) requires varying the duration (Vary Type = duration) of the repolarizing step (Vary Step = 3) while measuring the peak current evoked by the second command step (Measure Step = 4). Obtaining the steady-state inactivation (Figure 4) requires varying the magnitude (Vary Type = level) of the holding potential (Vary Step = 1) while measuring the peak current evoked by the command step (Measure Step = 2).

Figure 4.

Measuring the steady-state inactivation for the Hodgkin-Huxley model of the squid giant axon Na+ current. Na+ currents are activated by a 6 ms command step to −10 mV from a holding potential that varies from −100 mV to −20 mV in 10 mV increments. The steady-state inactivation is obtained from the sigmoidal peak current vs. holding step level graph on the lower left.

Using this interface, students are assigned the task of finding all parameters required to completely model the voltage-dependence and kinetic properties of the simulated currents. Table 1 shows the parameters used to characterize the voltage-dependence and kinetic properties of the Hodgkin-Huxley squid giant axon K+ and Na+ currents; a control button is provided to display curve fit results in the MATLAB command window (Save Data) that can be exported to a spreadsheet program for further analysis and display. To increase the level of difficulty, an option is provided to simulate non-inactivating currents (Inactivate = no) or inactivating currents (Inactivate = yes). These options simulate experiments in which one type of current has been pharmacologically isolated; a third option can be specified in which both current types are evoked simultaneously (Inactivate = both). When both types of current are simulated simultaneously, the user is required to subtract currents evoked from a depolarizing holding potential from currents evoked from a hyperpolarizing holding potential to isolate the inactivating current. Buttons are provided to allow the user to simulate currents evoked from a depolarizing holding potential (Subtract One or Subtract Seq); these traces will be subtracted from the next set of currents evoked from a hyperpolarizing holding potential (Simulate One or Simulate Seq).

Table 1.

Parameters for voltage-clamp protocols to obtain kinetic parameters from simulated Hodgkin-Huxley squid giant axon K+ and Na+ currents. To simulate and fit K+ currents, set Inactivation = no and Act Gate = 4. To simulate and fit Na+ currents, set Inactivation = yes and Act Gate = 3. (1) To obtain an accurate measure of m∞ from the peak current, the calculations must correct for the amount of inactivation at the time where the peak current is measured. (2) This protocol should be repeated using a range of repolarizing voltage values from −120 mV to −55 mV in 5 mV increments following the first pulse; substitute the appropriate voltage value for “xx” under Vstep. The value of τh can then be obtained as the time constant of the exponential peak current vs. step duration.

Parameter		Tstep	Vstep	Vary Type	Vary Step	Sequence	Measure Step	Relevant Results
K+	GmaxEK	2 8 6	−70 50 −70	Level	3	−120 : 5 : 50	32	Peak current
	n∞	2 15 3	−70 −20 −70	Level	2	−120 : 5 : 50	2	Steady-state
	τn	2 8 62 15 3	−70 50 −70−70 −20 −70	LevelLevel	32	−120 : 5 : −55−50 : 5 : 50	32	Time constantTime constant
Na+	GmaxENa	0.05 0.3 0.05	−120 50 −120	Level	3	−120 : 5 : 50	3	Peak current
	m∞	1 6 1	−120 50 −120	Level	2	−120 : 5 : 50	2	Peak currentPeak time (1)
	τm	0.05 0.3 0.050.1 0.3 0.41 6 1	−120 50 −120−120 50 −120−120 50 −120	LevelLevelLevel	332	−120 : 5 : −100−95 : 5 : −55−50 : 5 : 50	332	Time constantTime constantTime constant
	h∞	1 5 2	−120 0 −120	Level	1	−150 : 5 : −5	2	Peak current
	τh	0.2 0.5 1:10 201 6 1	−120 0 xx 0 −120−120 50 −120	DurationLevel	32	0.2 0.5 1:10 20−50 : 5 : 50	42	Peak current (2)Time constant

SINGLE-CHANNEL CURRENTS

To simulate ionic currents obtained from single-channel recordings, a MEX file was written to implement a computational model of stochastic currents. This file contains C++ code to simulate the random openings and closings of an individual ion channel. The model consists of various channel states (C for closed, O for open) with transition rates between states: where the transition rates kij provide a stochastic rate constant for a channel in state i at time t to make a transition to state j before time t + Δt. It can be shown that the dwell time, or the time spent continuously in state i, is a random variable with an exponential distribution that has a rate constant equal to the sum of all transition rates exiting state i (Colquhoun and Hawkes, 1981). In addition, the probability that the next transition will be from state i to state j is equal to the transition rate from state i to state j divided by the sum of all transition rates exiting state i (Colquhoun and Hawkes, 1982). Thus, simulated single-channel currents can be obtained by repetitive drawing of random number pairs, one to determine to which state the next transition will be, and the second to determine the dwell time in this state before another transition occurs. As with the whole-cell current model, parameters for the simulated single-channel currents are easily accessed in the code and can be changed by the instructor. Constants that define the single-channel conductance gsc, the reversal potential Erev, the number of states n and the transition rates kij are also accessible and can be easily changed. To increase the level of difficulty, rate constants can vary with ligand concentration and/or membrane potential. The instructor can also specify the amount of Gaussian noise that is added to the simulated currents. A text file accompanies the software that provides detailed instructions for changing model parameters.

Once the model was implemented and compiled as a MEX file, a GUI was constructed using the MATLAB programming environment to provide a straightforward interface for students to simulate single-channel currents (Figure 5). The main graph window displays simulated current traces from all trials from the most recent simulation. Currents from up to ten trials can be displayed simultaneously; a scroll bar is provided for viewing traces from remaining trials. To detect channel openings, an amplitude histogram of the current traces is constructed and Gaussian curve fits are obtained to determine the mean current amplitude for closed and open channels (solid green lines, lower left graph). The threshold for channel openings is initially calculated to be halfway between the mean closed and open channel current amplitude (dashed green line, lower left graph) and can be changed by the user. Regions of the current traces that were detected as channel openings are highlighted in red, and the number of openings and closings is displayed above the main graph window.

Figure 5.

MATLAB-based GUI to simulate single-channel recording using a stochastic model of a non-inactivating current. Current traces are displayed in groups of 10 within a scrolling window on the upper left; portions of the current trace that are detected as channel openings are highlighted in red. Channel openings are detected when the current crosses the amplitude threshold (dotted green line) that is halfway between the mean open and closed current amplitudes (solid green lines) obtained from Gaussian fits of the amplitude histogram on the lower left. Exponential fits of open and closed time histograms displayed on the lower middle and lower right are used to obtain mean open and closed times. The user can specify various parameters for the generation and analysis of simulated currents using edit boxes and buttons displayed on the upper right; see text for more details. The simulation parameters were specified to evoke currents over 100 trials of 500 ms each at a holding potential of −70 mV using a 6 ms command step that varies from −95 mV to −5 mV in 10 mV increments; the analysis parameters were specified to fit one exponential to the open time histogram and two exponentials to the closed time histogram.

To prevent erroneously long openings and closings from appearing in the open and closed channel histograms (lower middle and lower right graphs), users can visually inspect the traces to find short openings and closings that were not detected as threshold crossings. Figure 6 shows an example of a brief closing (middle trace, right) that was missed by the threshold detection algorithm. Because this event was not detected, a dwell time of 15 ms was added to the open time histogram. However, the user can click on this missed event (marked by ), and the 15 ms dwell time is removed from the open time histogram and replaced by two events of 7 ms and 8 ms in duration.

Figure 6.

Correcting for brief openings and closings that escaped threshold detection. Top three current traces from Figure 5 where magnified between 90 – 180 ms (icon for zoom tool not shown). A brief opening (green cross, middle trace, left) and closing (green cross, middle trace, right) were manually detected by the user clicking on these data points; event counts and dwell time histograms were adjusted accordingly. Also shown is the use of the square root – log axes for display of dwell time histograms; histogram peaks correspond to mean open and closed times obtained by exponential fits. See text for more details.

As with simulations of whole-cell currents, students are required to select various parameters for the acquisition and analysis of single-channel currents. For simulating single-channel currents, the user must provide the duration (Duration) and number (Trial Count) of trials, in addition to the ligand concentration (Ligand), membrane potential (Voltage), and whether the rate constants have ligand and/or voltage dependence (Dependence). For analysis of the results, the user must specify the number of exponentials for curve fits of dwell time histograms (Open exp, Closed exp). The user can also specify whether to display dwell time histograms on linear axes or on square-root vs. logarithmic axes that produce histogram peaks at the mean dwell times (Colquhoun and Sigworth, 1995). Users can also select the cutoff frequency of a low-pass filter for the current traces using a slider control. As with real experiments, low cutoff frequencies will attenuate fast openings and closings to prevent threshold detection, and high cutoff frequencies will obscure small amplitude openings with unfiltered noise. A control button is provided to simulate the specified number of trials (Simulate). Openings and closings for subsequent simulations will be added to the previous results unless model parameters are changed or the Reset button is pressed.

Using this interface, students are assigned the task of finding the rate constants required to completely model the simulated currents, which can be obtained from the amplitudes and time constants of the dwell time histogram curve fits. In addition, students must plot the magnitude of single-channel currents as a function of membrane potential to obtain the single-channel conductance gsc and the reversal potential Erev. To increase the level of difficulty, an option is provided to simulate fixed rate constants (Dependence = none), or rate constants that depend on ligand concentration (Dependence = ligand), membrane potential (Dependence = voltage) or both (Dependence = both). In the case of rate constants that vary with ligand or membrane potentials, students must try to identify which rate constants are changing and which remain constant. Students then try to reconstruct the state transition model of the simulated channel using the rate constants obtained from the analysis of the channel openings and closings.

DISCUSSION

The software described in this paper provides a useful tool to introduce students to experimental techniques used for obtaining biophysical measurements from excitable cells. Over the past three years, we have used this software for a cellular electrophysiology course developed for upper-level undergraduates in bioengineering. In addition, this software has also been used by a number of non-engineering students from a neurobiology program also enrolled in this course.

The primary role of this software is for the completion of two class projects. In the first project, students must obtain parameters to completely describe a Hodgkin-Huxley style model of whole-cell currents; in the second project, students must obtain parameters to completely describe a stochastic model of single-channel behavior. These projects could be adapted to various levels of difficulty and effort, such as running a few simulations to obtain the single-channel conductance or the steady-state inactivation, or a complete characterization of non-inactivating and inactivating currents obtained during the same recordings. Although there is no control group for comparison because this software has been used every year the electrophysiology course has been offered, we have found that the students have performed very well in homework and exam questions related to the mechanisms and characterization of electrophysiological phenomena. Such questions include an explanation of why K+ channel activation time constants must be slower than Na+ channel time constants for a cell to produce an action potential, or an estimation of the minimal number of states needed to explain data from representative single-channel traces without the use of dwell time histograms.

In addition, we have found this software provides a powerful tool for in-class lectures. Instead of relying on textbooks and lecture notes, we have started using the whole-cell current simulation software to teach students about the Hodgkin-Huxley model of squid giant axon currents. Following these lectures, model parameters are changed from the Hodgkin-Huxley non-inactivating K+ current and inactivating Na+ current to a non-inactivating Ca2+ current (Molitor and Manis, 1999) and an inactivating K+ current (Kanold and Manis, 1999) for class projects. Therefore, students learn the entire process of kinetic characterization for a standard dataset, and then apply this knowledge to a new dataset that had not been encountered previously. This approach has facilitated student understanding of the modeling process: whereas undergraduates were only able to analyze data from non-inactivating current in previous years, these students are now able to comprehend the kinetic processes and analyze data from inactivating currents as well. This approach has also improved the ability of students to understand the role of different classes of ion channels, and how these channels contribute to the electrophysiological responses of excitable cells.

Similarly, this software has been used for lectures on the analysis and modeling of single-channel currents. Initially, the software is used to present the appearance and characteristics of single-channel currents. After defining transition rates and deriving exponential dwell times using probability theory, the software is used to confirm the results of these derivations. Once the dwell times are obtained from the exponential fits, a kinetic model is constructed to explain the results. Following these lectures, students are assigned a class project in which they obtain transition rates that are either ligand dependent and/or voltage-dependent. As with the analysis and modeling of whole-cell currents, the use of this software for a sequence of lectures followed by a class project has facilitated student understanding of theoretical principles and data analysis techniques utilized for the modeling single-channel currents.