INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL METHODS RESULTS DISCUSSION APPENDIX REFERENCES

For many years Hodgkin-Huxley models (Hodgkin and Huxley, 1952a, b) have been the standard for describing ionic current kinetics. However, with the development of better recording techniques, new data have shown that these models have significant limitations. First, many single channel behaviors such as mean open time and first latency cannot be described using traditional Hodgkin-Huxley models. These behaviors can be estimated by expanding the Hodgkin-Huxley models to have multiple resting and inactivated states, but it is controversial as to how well these expanded models can predict single channel experimental data (Horn and Vandenberg, 1984; Chay, 1991). Second, while a Hodgkin-Huxley model can reproduce ionic currents, it does not necessarily correctly reproduce the underlying single channel kinetics. For example, Aldrich and co-workers found that for neuroblastoma sodium channels, activation has very slow components, while inactivation is fast (Aldrich et al., 1983). This finding contradicts Hodgkin's and Huxley's assumption that activation is rapid and inactivation is slow (Hodgkin and Huxley, 1952a, b). Even though Hodgkin-Huxley models can reproduce this current, they do not correctly reproduce the underlying channel kinetics. Although single channel recordings of cardiac sodium channels indicate that activation is rapid relative to inactivation (Yue et al., 1989), it is questionable whether Hodgkin-Huxley models are sufficient for reproducing behaviors that may be critically state-dependent, such as how ionic channels interact with drugs and toxins (Irvine and Winslow, 1996; Liu and Rasmusson, 1997). In addition, since much more is now known about the structure of the sodium channel (Noda et al., 1984, 1986; Guy, 1988), it is desirable to incorporate this information into a description of channel function. Thus, future channel models should be biophysically detailed kinetic models, consistent with current generalizations of channel structure, capable of reproducing single channel behavior.

For the cardiac sodium channel, models that describe channel gating as a Markov process (Benndorf, 1988; Berman et al., 1989; Scanley et al., 1990) are a step in this direction. Existing models for the cardiac sodium channel are, however, incomplete in that they describe only certain features of channel behavior. Specifically, each of these models lacks rate constants with explicit voltage and temperature dependence. In addition, these models treat inactivation as an absorbing state, so that once a channel inactivates, there is no pathway by which it can recover. Thus, they can only be used to simulate certain channel behaviors in response to a single voltage clamp stimulus. They do not reproduce channel activity to the same extent as Hodgkin-Huxley models and therefore have not been as widely used.

More comprehensive Markov models exist for sodium channels of the squid giant axon (Patlak, 1991; Vandenberg and Bezanilla, 1991). Vandenberg and Bezanilla and Patlak have been able to develop these models by using a wide variety of both whole cell and single channel data simultaneously. Unfortunately, because of the many differences in channel kinetics between cardiac and neuronal tissue (Kirsch and Brown, 1989; Kuo and Bean, 1994; Hanck and Sheets, 1995; Fozzard and Hanck, 1996), these models cannot be used directly to model cardiac sodium channels. Nevertheless, the neuronal models and the techniques used to develop them are a starting point from which to develop a more complete model of the cardiac sodium channel.

Although Markov models exist from which cardiac sodium channel Markov models can be developed for a single temperature, no models exist that can reproduce ensemble-average and single channel behaviors for a range of temperatures. The models of Vandenberg and Patlak have a temperature-dependent rate constant coefficient and a temperature-dependent voltage term (Patlak, 1991; Vandenberg and Bezanilla, 1991). Changing the temperature in these terms, however, does not yield the correct channel activity at multiple temperatures. Each term in the model's rate constants needs to have its own temperature dependence or its own Q₁₀ factor (Kohlhardt, 1990). Temperature-dependent closed-closed and closed-open transitions have been incorporated into a partial model of neuronal sodium channels by formulating the rate constants as exponential functions of enthalpy and entropy (Correa et al., 1992). The same rate constant formulation can be used in a model of cardiac sodium channels to reproduce ensemble-average and single channel behaviors for a range of temperatures.

The goal of this study is to use neuronal models as a framework for developing a Markov model of the cardiac sodium channel. The model should exhibit correct macroscopic and single channel behavior, including recovery from inactivation, for a voltage range of 150 mV to 20 mV and a temperature range of 10°C to 25°C. Such a model would improve on existing Hodgkin-Huxley and Markov models significantly and may yield more insight into the molecular basis of channel function. In addition, such a model could be used as the basis for studies of antiarrhythmic drug action.

	MODEL

TOP ABSTRACT INTRODUCTION MODEL METHODS RESULTS DISCUSSION APPENDIX REFERENCES

The cardiac sodium channel Markov model is patterned after that by Kuo and Bean for sodium channels in CA1 hippocampal neurons (Kuo and Bean, 1994). This model is chosen as a starting point because it is consistent with current generalizations of channel structure, but uses symmetry and cooperative movement of the voltage sensors to reduce the number of free parameters. As shown in Fig. 1, the channel can occupy any of 13 states. The top row of states corresponds to zero to four voltage sensors being activated (C₀ through C₄) plus an additional conformational change required for opening (C₄  O₁ and C₄  O₂). The bottom row of states corresponds to the inactivation particle blocking the pore at each position of the voltage sensors. As in Kuo and Bean's model, the affinity of the inactivation particle binding site is hypothesized to increase by a scaling factor (a) as the channel activates and to decrease by the same factor as the channel deactivates. Closed-closed and closed-open transitions (horizontal transitions) are voltage-dependent, and closed-inactivated transitions (vertical transitions) are voltage-independent (Kuo and Bean, 1994).

View larger version (12K): [in this window] [in a new window]   FIGURE 1   State diagram for the cardiac sodium channel Markov model. C0-C4 are closed states, O1 and O2 are open states, C0I-C4I are closed-inactivated states, and I is the inactivated state. All rate constants are voltage- and temperature-dependent except for those governing transitions between closed and closed-inactivated states, which are only temperature-dependent.

In order to represent the cardiac sodium channel more accurately, two modifications are made to the Kuo and Bean model. The first modification is that an additional open state (O₂), with the same conductance as the first, is added. Transitions between the two open states are voltage-independent. The addition of a second open state provides another pathway by which the channel can open (C₄  O₂) and improves the fit to the decay of the ionic currents. Two arguments can be made for the existence of more than one open state. First, although single channel open time distributions are generally well fit by a single exponential (Patlak and Oritz, 1985; Berman et al., 1989; Scanley et al., 1990), they can also be fit well by multiple exponentials, particularly in the presence of toxins (Kunze et al., 1985; Nagy, 1987; Schreibmayer and Jeglitsch, 1992; Correa and Bezanilla, 1994). Second, sodium channels from many tissues, including cardiac tissue, produce tail currents with two exponential components (Goldman and Hahin, 1978; Dubois and Schneider, 1982; Hanck and Sheets, 1995; Elinder and Arhem, 1997). Erlinder and Arhem suggest that, in the absence of a two-step deactivation, in which the first transition is a fast equilibrium and the second is slow, this biexponential decay can only be produced by two open states connected by different pathways to a common closed state (Elinder and Arhem, 1997).

The second modification of the Kuo and Bean model is that open-inactivated transitions are made voltage-dependent. The change in the open-inactivated rate constants is supported by Yue, Lawrence, and colleagues' finding that a voltage-dependent open-to-inactivated transition is necessary to produce the correct voltage dependence of channel reopenings and mean open times (Yue et al., 1989; Lawrence et al., 1991). It is also supported by Sheets' and Hanck's measurement of a significant component of gating current due to this transition (Sheets and Hanck, 1995).

Rate constants are of the form from Eyring rate theory (Hille, 1992)

(1)

where k is the Boltzmann constant, T is the absolute temperature, h is the Planck constant, R is the gas constant, F is Faraday's constant, H is the change in enthalpy, S is the change in entropy, z is the effective valence (i.e., the charge moved times the fractional distance the charge is moved through the membrane), and V is the membrane potential in volts. By convention, along the top row, all transitions toward an open state have positive valences because they are favored by depolarization, while those away from an open state have negative valences because they are favored by repolarization. The same convention is used along the bottom row; transitions toward the inactivated state have positive valences, while those away from the inactivated state have negative valences.

There are several loops in the model that must satisfy microscopic reversibility. Microscopic reversibility is derived from the law of conservation of energy and states that the product of rate constants when traversing a loop clockwise must be equal to the product when traversing the same loop counterclockwise (Hille, 1992). For the closed-closed-inactivated loops, satisfying microscopic reversibility requires that the transitions among the closed-inactivated states be scaled by a, the same factor used to scale the transitions between rows. Microscopic reversibility is preserved around the closed-open-inactivated loop by isolating the H, S, and z terms in the product and satisfying each term separately using the following equations:

(2)

(3)

(4)

Similarly, microscopic reversibility is preserved around the closed-open-open loop using the following equations for H, S, and z:

(5)

(6)

(7)

	METHODS

TOP ABSTRACT INTRODUCTION MODEL METHODS RESULTS DISCUSSION APPENDIX REFERENCES

Model development

The probability of occupying any particular channel state is described mathematically by a set of ordinary differential equations, written in matrix notation as

(8)

where P(t) is a vector describing the probabilities of occupying each state and W is the transition matrix. In general, W will be a function of voltage and thus time. For voltage-clamped conditions, however, W is time-independent; thus Eq. 8 has the analytic solution

(9)

Equation 9 is solved on a Silicon Graphics computer using linear algebra subroutines from the Silicon Graphics mathematics library (complib.sgimath).

Parameters of the model are determined using a simulated annealing algorithm (Corana et al., 1987). This algorithm minimizes the cost function, which is the weighted sum of the least-squared errors between model responses and experimental data, by randomly searching the parameter space and incrementally decreasing the search radius. Whereas many minimization algorithms accept only downhill moves and tend to converge on local minima, the simulated annealing algorithm accepts uphill moves as well, and thus is more likely to find the global minimum. Uphill moves are accepted based on the Metropolis criterion, a probabilistic function determined from the difference between the new and old errors and the annealing temperature. The annealing temperature controls the rate of convergence by influencing what uphill moves are accepted and by limiting the search radius. In order to reach a minimum, the annealing temperature is decreased by 5% per 50 N function evaluations, where N is the number of parameters to be determined, as the algorithm converges on a solution. The algorithm is terminated when there is no more than 0.1% change in error since the last temperature reduction.

In order to limit the number of free parameters to be determined during each minimization, the fitting procedure is done in parts. First, the enthalpy and entropy terms are collapsed into a single Gibbs free energy term and the Gibbs free energies and effective valences are determined for a temperature of 13°C. Then, holding the effective valences constant, the enthalpy and entropy terms are determined. The entropy terms are written in terms of the enthalpy and the Gibbs free energy (G):

(10)

where T is the temperature, 286 K. Since the Gibbs free energies and valences are known from the previous minimization, substitution of Eq. 10 into Eq. 1 leaves only the enthalpies to be determined. The enthalpies are determined by fitting experimental data at 21°C using the simulated annealing algorithm.

As shown by Vandenberg and Bezanilla (1991) and Patlak (1991) in developing models of sodium channels in squid giant axon, a variety of experimental data sets are needed to fully determine the model parameters. In this study, the experimental data for 13°C include ionic currents (provided by Hanck and Sheets, similar to Sheets et al., 1996), gating charge accumulation (Hanck and Sheets as above), the steady-state inactivation curve (Hanck and Sheets as above), the rate of tail current relaxation (Hanck and Sheets, 1995), the time course of recovery from inactivation (Sakakibara et al., 1993), and single channel open times (Sheets and Hanck, 1995). The majority of data are taken from hH1 sodium channels or, where these data are unavailable, from canine sodium channels. Recovery data at 13°C are unavailable, so data at 17°C are used to approximate the data at 13°C. This approximation is acceptable because the difference in recovery rate between 13°C and 17°C is probably similar to the variation in recovery rate among cells at a single temperature. Ionic currents are calculated as

(11)

where I_Na is the sodium current, G_Na is the maximal channel conductance, P_open is the probability of occupying the open states (O₁ + O₂), V is the membrane potential, and E_Na is the reversal potential for sodium. G_Na is a function of temperature and thus is a parameter to be determined at both 13°C and 21°C. E_Na is dependent upon the experimental solutions, which are different for the data at 13°C and 21°C, and so is set accordingly at each temperature. Gating current is calculated according to the formula (Vandenberg and Bezanilla, 1991):

(12)

where n is the number of channels, e is the elementary charge unit, z is the effective valence, P_j is the probability of occupying state j, and _jk is the rate constant for the transition from state j to state k. Gating charge is found by integrating the gating current. Ionic currents and gating charge are computed using the following protocol. The membrane potential is held at 150 mV and then stepped for 20 ms to a potential between 70 mV and 20 mV inclusive in 10-mV increments. To eliminate convergence problems introduced by experimental error for potentials 20 mV, gating charge accumulation curves are fit with a single exponential function and the curve fit values, instead of the experimental data, are used for the plateau portion of the gating charge accumulation curves. Tail currents are computed by stepping from 150 mV to 40 mV until the current reaches its maximal value (after 0.94 ms) and then stepping down to potentials between 150 mV and 90 mV inclusive in 10-mV increments for 5 ms. Recovery from inactivation is assessed using a double-pulse protocol. The membrane potential is held at 140 mV and then stepped to 20 mV for 1 s. The potential is then stepped down to either 100 mV or 140 mV for lengths of time varying from 5 to 600 ms. Current is then measured during a 4-ms step to 0 mV to assess the amount of recovery. Open time distributions are calculated for potentials between 90 and 10 mV inclusive in 10-mV increments using a simplified model in which all transitions out of the open states are to an absorbing non-open state. The simplified model and the equation for the open time distributions are shown in the Appendix. Each data set is weighted so that all sets have approximately the same influence on the cost function and so that the parameters determined by the algorithm are those parameters which best reproduce all of the channel kinetics. The weights for ionic current, gating charge, steady-state inactivation, tail current, recovery from inactivation, and open time cost function terms are 1, 250, 1, 500, 1000, and 5000, respectively.

The experimental data for 21°C include ionic currents (Wang et al., 1996), gating charge accumulation (Josephson and Sperelakis, 1992), the steady-state inactivation curve (Wang et al., 1996), the time course of recovery from inactivation (Wang et al., 1996), and single channel open times (Benndorf, 1988). To compute the ionic currents, the membrane potential is held at 120 mV and then stepped for 15 ms to a potential between 60 mV and 20 mV inclusive in 10-mV increments. The same protocol is used to compute gating charge accumulation except that the holding potential is 150 mV. From measurements of gating charge in squid giant axon, the maximum charge displaced at each potential does not vary with temperature (Jonas, 1989). Therefore, the model's computed maximum charge values for 13°C were used as the experimental charge values for 21°C. Recovery from inactivation is again measured using a double-pulse protocol. The holding potential is 120 mV and the test potential is 20 mV; recovery times are varied between 10 ms and 250 ms. Open time distributions are calculated for potentials between 70 and 20 mV inclusive in 10-mV increments. The weights for ionic current, gating charge, steady-state inactivation, recovery from inactivation, and open time cost function terms are 100, 10, 5000, 500, and 500, respectively.

Model testing

The single channel behavior of the model was tested, which required the state transitions to be determined using a stochastic approach (Clay and DeFelice, 1983). In this method, the length of time a channel stays in its current state (i.e., its dwell time, denoted as T_j) is calculated according to the formula

(13)

where r is a random number from the uniform distribution [0, 1] and _jk is the transition rate from state j to state k. The sum is over the x pathways out of state j. At the end of the dwell time, the new state of the channel is determined by assigning random numbers to a portion of the interval [0, 1] based on the probabilities of changing to neighboring states. These probabilities are equal to the rate constant for a particular transition divided by the sum of the rate constants for all possible transitions. For example, a channel in state C₁ can transit to C₂ or C₁I. The probability of changing to C₂ is 4/(C_n + 4), where 4 is the rate constant for C₁  C₂ and C_n is the rate constant for C₁  C₁I. Once the new state is determined, another random number is used to calculate the dwell time in the new state. At an instantaneous voltage step, channels remain in their current state, but the dwell times are recalculated.

	RESULTS

TOP ABSTRACT INTRODUCTION MODEL METHODS RESULTS DISCUSSION APPENDIX REFERENCES

The model parameters determined to give the best total fit of model responses to experimental data for ionic currents, gating charge, steady-state inactivation, tail currents, recovery from inactivation, and open times are listed in Table 1. The rate constants governing the O₁  C₄ (deactivation) and the O₁  I (inactivation) transitions have been measured experimentally. Benndorf and Koopmann found the enthalpies, entropies, and effective valences to be 129 kJ mol¹, 0.23 kJ mol¹ K¹, and 1.54 for deactivation and 79 kJ mol¹, 0.10 kJ mol¹ K¹, and 0.68 for inactivation, respectively (Benndorf and Koopmann, 1993). The model parameters are 128 kJ mol¹, 0.229 kJ mol¹ K¹, and 1.33 for deactivation and 62 kJ mol¹, 0.039 kJ mol¹ K¹, and 0.66 for inactivation, respectively. The model parameters for deactivation are very similar to the experimental data, whereas the parameters for inactivation differ slightly from the experimental data. The Gibbs free energies for inactivation for the experimental data and the model parameters, however, are very similar. Therefore, the discrepancy in inactivation parameters probably results from the minimization algorithm not being able to discriminate between several pairs of enthalpy and entropy terms yielding the same Gibbs free energy.

To assess the sensitivity of the model parameters, each parameter is varied by ±1% of its value and the change in the cost function is computed. A change in the value of the Gibbs free energy produces a much larger change in the cost function than does a change in the value of the corresponding valence. Thus, the Gibbs free energies are much more sensitive to a change in their values than are the valences. This sensitivity difference can be attributed to the different importance Gibbs free energies and valences have in determining the rate constants. The Gibbs free energies are the larger of the two terms in the exponential function and therefore are mainly responsible for determining the rate constant. In contrast, the voltage-dependent terms serve only to slightly modify the basic rates set by the Gibbs free energies. Thus, by changing the Gibbs free energies, one can produce a much larger change in the rate constant and a much larger change in the cost function.

Changes in the enthalpy and entropy terms of , , and C_n produce the largest changes in the cost function. The total error increases by up to 10 times for a 1% change in these parameters. The large sensitivity of these parameters is probably due to their role in providing temperature dependence. Parameters  and  must have precise enthalpy and entropy terms in order to accurately describe the increased rate of channel activation and rate of recovery from inactivation with temperature. C_n requires precise enthalpy and entropy terms in order to describe the shift in the steady-state inactivation curve with temperature. At 13°C, a change in the enthalpy of  also produces a large change in the cost function. However, at 21°C, the same change produces much less change in the cost function, because the probability of occupying the second open state is much lower at this temperature.

In contrast to the large sensitivity of , , C_n, and , a 1% change in the enthalpy and entropy terms of , , and O_f, as well as in the entropy term of C_f, produces almost no change in the total error. Two different explanations account for this small sensitivity. First, the entropies of O_f and C_f are a much smaller fraction of their respective enthalpies than are the entropies of other parameters. Thus, based on their relative size alone, one would expect a 1% change in their values to have little influence on the total error. Second, a change in the enthalpy and entropy terms of , , and O_f most likely produces little effect on the error because there is not enough information in the experimental data to adequately constrain two separate terms. This observation is in contrast to changing the Gibbs free energy term for these parameters, which does produce a measurable change in the total error. These results suggest that while the entropy and enthalpy terms of , , and O_f may not be well determined, the Gibbs free energy they define is. In other words, there may be several combinations of entropy and enthalpy terms that produce an appropriate Gibbs free energy for , , and O_f.

The parameter values can be used to assess the amount of charge moved with each transition in the model. The model parameters suggest that activation requires the movement of 6.8 charges and inactivation requires the movement of 0.66 charges. Estimates of the charge associated with activation usually range from 4 to 7 (Hodgkin and Huxley, 1952a; Oxford, 1981; Sheets and Hanck, 1995), although some researchers have found that at least 12 charges are needed (Hirschberg et al., 1995). Estimates of the charge associated with inactivation range from 0.75 to 1.9 charges (Horn et al., 1984; Vandenberg and Horn, 1984; Yue et al., 1989; Lawrence et al., 1991; Sheets and Hanck, 1995). Thus, the model's estimates of charges required for activation and inactivation are similar to values measured experimentally.

Most of the charge movement in the activation pathway is concentrated in the last transition (C₄  O₁ or C₄  O₂). This finding seemingly contradicts the hypothesis that the final transition in the activation pathway is voltage-independent for all voltage-gated channels. However, these transitions in the model probably represent several steps lumped together so that, in reality, the final step may really be voltage-independent (Kuo and Bean, 1994). Furthermore, gating currents, which depend heavily on the voltage-dependence of each transition, have only been used in developing a few models. In another sodium channel model developed using gating currents, the closed-to-open transition also has the greatest voltage dependence (Vandenberg and Bezanilla, 1991).

The model is able to reproduce a wide range of experimental data. Fig. 2 shows representative traces of the model-derived ionic current in comparison to the experimental data at 13°C (provided by Hanck and Sheets, similar to Sheets et al., 1996) and at 21°C (provided by Bennett, similar to Wang et al., 1996) for clamp voltages of 50 mV, 30 mV, 10 mV, and 10 mV. Although the peak current values deviate slightly from the experimental values, the time courses of activation and inactivation are generally well fit by the model at both temperatures. For potentials >40 mV, the model-derived currents decay to zero within 50 ms, whereas the experimental currents do not. The model decays to zero faster than the experimental data at 13°C because the experimental data have both fast and slow components of inactivation, whereas the model has only fast inactivation.

View larger version (14K): [in this window] [in a new window]   FIGURE 2   Comparison of voltage-clamped sodium current tracings for clamp voltages of 50 mV, 30 mV, 10 mV, and 10 mV for experimental data (dashed lines) and the model (solid lines). (A) At 13°C (experimental data provided by Hanck and Sheets similar to Sheets et al., 1996). (B) At 21°C (experimental data provided by Bennett similar to Wang et al., 1996).

Fig. 3 compares model-derived current-voltage relationships at 13°C (Fig. 3 A) and 21°C (Fig. 3 C) with experimental data used in the fitting process. Fig. 3 B compares the current-voltage relation predicted at 17°C with experimental data (provided by Wasserstrom, similar to Sakakibara et al., 1993). As shown by the current-voltage relationships, the model reproduces peak current well throughout the voltage and temperature range tested. At 13°C, the current magnitudes deviate most significantly from the corresponding experimental data over the voltage range 60 mV to 45 mV. In this range, the model produces less current than the experimental data. This reduction in current probably results from inactivation that is too fast, which reduces the number of channel reopenings. At 17°C, peak currents also differ over this voltage range. Considering, however, that data at this temperature were not fit and were obtained from a different experimental preparation than any of the other data, the model reproduces the data well. At 21°C, peak currents are well fit except at very depolarized potentials. Both the model and experimental current-voltage curves are fit using a modified Boltzmann function

(14)

where G_Na is the conductance, V is the voltage, E_Na is the reversal potential, V_0.5 is the potential at which the current is half-maximal, and s is the slope factor. The resulting parameters are listed in Table 2. The conductances and slope factors are similar with the model having a slightly shallower slope for 13°C and 17°C and a slightly steeper slope for 21°C. The model's half-maximal voltages differ from the corresponding experimental data by 3.1 mV, 6.1 mV, and 1.2 mV for 13°C, 17°C, and 21°C, respectively. This difference is reflected in the slightly rightward shift of the model's current-voltage curve for 13°C and the slightly leftward shift of the model's current-voltage curve for 17°C. The model shows a rightward shift of the half-maximal potential (11.2 mV per 8°C) as the temperature is increased. Rightward shifts in the half-maximal potentials of the steady-state activation (8 mV per 10°C) and inactivation curves (7 mV per 10°C) as the temperature is increased have been measured experimentally (Murray et al., 1990). These shifts would produce a shift in the current-voltage relationship similar to that produced by the model.

View larger version (19K): [in this window] [in a new window]   FIGURE 3   Normalized current-voltage curves for experimental data and the model. Curves are the best fits to Eq. 14 in the text. See Table 2 for fitted parameters. (A) At 13°C experimental data () and the model () (experimental data provided by Hanck and Sheets similar to Sheets et al., 1996). (B) At 17°C experimental data () and the model () (experimental data provided by Wasserstrom similar to Sakakibara et al., 1993). (C) At 21°C experimental data () and the model () (experimental data provided by Bennett similar to Wang et al., 1996).

Fig. 4, A and B show the time to peak current and the time constants of inactivation for 13°C, 17°C, and 21°C. The model's time to peak current is very similar to the corresponding experimental data for all the voltages and temperatures tested. As the temperature is increased, the time to peak current is reduced, as expected from data in the literature (Colatsky, 1980; Murray et al., 1990). The time constants of inactivation are estimated by fitting an exponential function to the current decay. For potentials of 30 mV and greater, the model's ionic current decay is better fit with two exponentials. The experimental data also have a second exponential at these potentials. However, the time constant of this second exponential in the experimental data is too large to accurately determine with a voltage clamp of 35 ms. Therefore, in order to compare the model's time constants with those of the experimental data, a single exponential fit is used in Fig. 4 B. At 13°C, the inactivation time constants predicted by the model are larger than those of the experimental data at both ends of the voltage range, while they are similar to those of the experimental data in the middle of the voltage range. The discrepancy in the time constants at negative potentials is probably due to error in the fitting procedure as a result of having only two or three time constants' worth of data. In addition, the model's time constants increase for very depolarized potentials, whereas the experimental data show that the time constants decrease monotonically as voltage is increased. This apparent discrepancy is an artifact of fitting the model's decay with a single exponential. As the potential is increased, the model's decay switches from that described by a large fast component and a small slow component to one described by a small fast component and a large slow component. As the slower component becomes larger, the time constant determined by fitting a single exponential to the decay increases. If just the fast or slow time constant from a biexponential fit of the model's decay is plotted, the time constants do decrease monotonically with increasing voltage. At 17°C and 21°C, the model's time constants of inactivation are similar to the corresponding experimental data. As the temperature is increased, the time constants become faster, as expected from data in the literature (Colatsky, 1980). Thus, the model reproduces well the activation and inactivation properties of the channel for a large voltage and temperature range.

View larger version (24K): [in this window] [in a new window]   FIGURE 4   Temperature dependence of activation and inactivation at 13°C (experimental data (), model ()), at 17°C (experimental data (), model ()), and at 21°C (experimental data (), model ()). Experimental data from sources listed in Fig. 3. (A) Time to peak current. (B) Time constants of inactivation determined by fitting a single exponential to the current decay.

The second data set reproduced by the model is gating charge movement in response to voltage-clamp stimuli. Fig. 5 A shows a comparison of the charge-voltage curves for the model and the experimental data (provided by Hanck and Sheets, similar to Sheets et al., 1996) at 13°C. The magnitudes plotted here for the model are the maximum charge accumulated at 30 ms. Some error is associated with these values because charge is still being accumulated at a very slow rate at 30 ms; that is, in the model, the plateau portion of the charge accumulation curves is not completely flat, but rather has a small slope. The nonequilibrium movement of charge is most likely due to constraints imposed by microscopic reversibility on the closed-open-inactivated loop. Nevertheless, the magnitude of charge movement is well reproduced by the model over most of the voltage range tested. Fitting a Boltzmann function to the experimental data at 13°C yields slope factor and half-maximal potential values of 15.7 mV and 63.3 mV. However, the experimental data for 70 mV and 60 mV appear to deviate significantly from the remainder of the data. Exclusion of these two points results in a better fit (correlation coefficient 0.996 versus 0.984). The slope factor and half-maximal potential values for this fit are 15.8 mV and 70 mV. The model's half-maximal potential is identical to that of the experimental data (70 mV) and its slope is slightly steeper (14 mV). Fig. 5 B compares the model's charge-voltage curves for 13°C and 21°C. At 21°C, maximum charge values are taken at 20 ms even though charge is still being accumulated at a very slow rate. Fitting a Boltzmann function to the 21°C curve yields slope factor and half-maximal potential values of 19.6 mV and 63.9 mV, respectively. Thus, with increased temperature, the charge-voltage curve shifts rightward by 6.1 mV per 8°C. Although data on the temperature dependence of the charge-voltage curve for heart tissue have not been published, Hanck and co-workers have shown that the charge-voltage curve, although having a steeper slope, is similar to the peak sodium conductance curve (Hanck et al., 1990). Thus, one would expect the charge-voltage curve to be shifted rightward with increasing temperature since the peak sodium conductance curve is shifted rightward with temperature (Murray et al., 1990). Therefore, the model can reproduce the charge-voltage relationship over a large range of voltages and temperatures.

View larger version (21K): [in this window] [in a new window]   FIGURE 5   Normalized gating charge-voltage curves. (A) At 13°C for experimental data () (provided by Hanck and Sheets similar to Sheets et al., 1996) and the model (). Curves are the best fits to a Boltzmann function where the slope factor and half-maximal potential values are 15.7 mV and 63.3 mV for the experimental data and 14.0 mV and 70 mV for the model, respectively. (B) Model's charge-voltage curves at 13°C () and at 21°C (). The slope factor and half-maximal potential values for 21°C are 19.6 mV and 63.9 mV, respectively.

The model can also approximate the rate of gating charge movement. Fig. 6 shows the gating charge accumulation time constants as a function of voltage for the model and experimental data (provided by Hanck and Sheets, similar to Sheets et al., 1996) at 13°C and for the model at 21°C. For potentials 50 mV, gating charge accumulation is well fit by a single exponential. For potentials <50 mV, the model's gating charge accumulation curves exhibit an initial fast decay followed by a much slower return to zero. Since the gating charge accumulation is thus not well fit by a single exponential, these potentials were excluded from Fig. 6. Note that the experimental time constants plotted in Fig. 6 are those obtained from a single measurement of gating charge from one cell and are from a different cell than the ionic currents. Taking these facts into consideration, at 13°C, although the model has larger time constants at each potential, it approximates the voltage dependence of these time constants reasonably well. At 21°C, one would expect the time constants of gating charge accumulation to be much faster (Josephson and Sperelakis, 1992) and the model meets this expectation with ~1 ms reduction in the time constants throughout the voltage range.

View larger version (23K): [in this window] [in a new window]   FIGURE 6   Time constants of gating charge accumulations (experimental data at 13°C (), model data at 13°C (), model data at 21°C ()). Experimental data from source listed in Fig. 5.

A third data set reproduced by the model is steady-state availability. Fig. 7 A shows the steady-state availability curves for the model and the experimental data (provided by Hanck and Sheets) at 13°C. The model's curve is nearly identical to that of the experimental data as reflected in the slope factor and half-maximal potential values of the respective Boltzmann functions. For the experimental data, the slope factor and half-maximal potential values are 9.9 mV and 106.1 mV. For the model, the respective values are 10.6 mV and 107 mV. Fig. 7 B compares the steady-state availability curves at 13°C and 21°C. At 21°C, there is a noticeable rightward shift of the curve. Fitting a Boltzmann function to the 21°C curve yields slope factor and half-maximal potential values of 15.2 mV and 101.7 mV. Thus, as the temperature is increased, the model produces a rightward shift of 5.3 mV per 8°C. This shift is similar to that measured experimentally, which is 7 mV per 10°C (Murray et al., 1990). Experimental data also predict that the slope factor increases slightly (0.5 mV per 10°C) over the temperature range 16°C to 26°C (Murray et al., 1990). The additional increase in slope factor over that predicted by the experimental data probably results from the exponential form of the rate constants, which prevents them from becoming constant except at the extremes of the voltage range. The nonsaturating rate constants prevent a true plateau of the curve at very negative potentials and thus, as temperature is increased, the curve cannot simply be shifted to the right. Nevertheless, the model is able to reproduce well the level of inactivation over a large voltage and temperature range.

View larger version (21K): [in this window] [in a new window]   FIGURE 7   Steady-state availability curves. (A) At 13°C for experimental data () (provided by Hanck and Sheets similar to Sheets et al., 1996) and the model (). Curves are the best fits of a Boltzmann function where the slope factor and half-maximal potential values are 9.9 mV and 106.1 mV for the experimental data and 10.6 mV and 107 mV for the model, respectively. (B) Model's steady-state availability curves at 13°C () and at 21°C (). The slope factor and half-maximal potential values for 21°C are 15.2 mV and 101.7 mV, respectively.

In addition to ionic currents, gating charge, and steady-state availability, the model can reproduce recovery from inactivation. Shown in Fig. 8 A is the rate at which the model recovers from inactivation at 13°C in comparison to experimental data (Sakakibara et al., 1993). The model recovers from inactivation at a similar rate as the experimental data for a holding potential of 100 mV and, as the holding potential is decreased, the model recovers faster. For 100 mV, both the experimental data and the model show a delay of 10 ms before recovery occurs. Sakakibara and co-workers fit the experimental data using the sum of two exponentials. However, their data can be fit just as well with a single exponential, and this form yields better fits to the model results. For 100 mV the time constants are 164.9 ms and 177.7 ms, for 120 mV they are 42.9 ms and 46.7 ms, and for 140 mV they are 6.9 ms and 21.9 ms for the experimental and model data, respectively. The model's time constants are generally larger than those of the experimental data and, as the holding potential is decreased, the difference between the time constants increases. One possible explanation for the model's recovery rate being too slow is a lack of voltage dependence of the rate constants governing transitions between the closed and closed-inactivated states. Adding voltage dependence here could allow the model to more correctly approximate the time constant of recovery as holding potentials are lowered, but such a change would increase model complexity substantially by adding additional loops for which microscopic reversibility must be satisfied. Fig. 8 B compares the rates of recovery from inactivation at 13°C and 21°C for a holding potential of 120 mV. As temperature is increased, the recovery rate is increased significantly. Fitting a single exponential to the recovery data at 21°C yields a time constant of 13.3 ms. This rate is similar to that found by refitting published data at 21°C (Wang et al., 1996) with a single exponential, which yields a time constant of 15 ms. The model presented here thus differs from existing Markov models of the cardiac sodium channel (Benndorf, 1988; Berman et al., 1989; Scanley et al., 1990) in that it recovers from inactivation with the correct voltage dependence.

View larger version (27K): [in this window] [in a new window]   FIGURE 8   Rates of recovery from inactivation. (A) At 13°C for experimental data (Sakakibara et al., 1993) and the model. Data are plotted for holding potentials of 100 mV (experimental data (), model ()), 120 mV (experimental data (), model ()), and 140 mV (experimental data (), model ()). Model curves are fit using a single exponential with time constants of 177.7 ms, 46.7 ms, and 21.9 ms for 100 mV, 120 mV, and 140 mV, respectively. (B) Model's recovery from inactivation curves for a holding potential of 120 mV at 13°C () and 21°C (). The time constant at 21°C is 13.3 ms.

The model also reproduces the rate at which an open channel deactivates. Plotted in Fig. 9 are the time constants from a single exponential fit of the current decay at 13°C upon stepping from 40 mV to the test potential. The model has two deactivation pathways (O₁  C₄ and O₂  C₄) and therefore will have a biexponential tail current. For potentials of 100 mV and below, the model's time constants are similar to those measured experimentally (Hanck and Sheets, 1995) and there is little difference between a monoexponential and a biexponential fit to the data. At these potentials, the O₁  C₄ pathway, which is faster, appears to dominate. At more depolarized potentials, the model predicts larger time constants than measured experimentally and the current decay is much better fit using two exponentials. Thus, as the test potential is increased, channels more readily exit the open states using both deactivation pathways. At 21°C, the current decay produced by the model is best fit using a single exponential at all potentials. At this temperature, the O₁  C₄ pathway dominates because of the low probability of occupying the second open state. As expected, as the temperature is increased, the rate of deactivation increases.

View larger version (22K): [in this window] [in a new window]   FIGURE 9   Time constants from a single exponential fit of the tail current relaxations at 13°C for experimental data (Hanck and Sheets, 1995) () and the model () and at 21°C for the model ().

The final data set used to determine the model parameters is single channel open durations. Fig. 10 shows the model's densities of single channel open durations for both 13°C and 21°C. Twelve hundred channels are simulated, as described in Methods, for a 40-ms sweep and their open durations measured. At 13°C, for a clamp voltage of 50 mV, a histogram of open durations shows a wide variation of open times including a significant fraction >2 ms. In contrast, for a clamp voltage of 15 mV at 13°C, most of the channels have open times of <2 ms. These densities are clearly biexponential because the model has two open states and the dwell times in each open state are significantly different. At 21°C, for both clamp potentials, almost all of the channels have open times of <2 ms, but the densities are still biexponential. Experimental data for the densities and distributions of open times are usually fit with a single exponential. In order to compare the model and the experimental data, the distributions are calculated from Eq. A14 and are fit with a single exponential. Fitting the densities calculated from stochastic channel simulations yields the same results.

View larger version (23K): [in this window] [in a new window]   FIGURE 10   Densities of single channel open durations. Twelve hundred channels are simulated for 40 ms as described in the text. Bin size is 0.5 ms. (A) At 13°C for 50 mV. (B) At 13°C for 15 mV. (C) At 21°C for 50 mV. (D) At 21°C for 15 mV.

Fig. 11 shows the open time distribution time constants versus voltage at 13°C and 21°C and the model's prediction at 17°C. Experimental data by Scanley et al. (1990) are plotted for 13°C and by Benndorf (1988) for 21°C. For the entire voltage and temperature range depicted, the model data agree well with the experimental data. As temperature is increased, the model's open times become shorter and the peak open time is shifted rightward, both of which are supported by the experimental data. It should be emphasized that time constants obtained by fitting a single exponential to the open time distributions are not equal to the mean open times of the model. (The mean open times can be calculated using the probability density function in Eq. A13.) The mean open time at each potential is larger, particularly at very depolarized potentials, due to rare long occupancies of the second open state (see Fig. 10). Despite its larger mean open times, the model reproduces the distributions of the open durations well for a large voltage and temperature range.

View larger version (30K): [in this window] [in a new window]   FIGURE 11   Model-predicted time constants for the open time distributions calculated from Eq. A14 at 13°C (solid line), 17°C (dashed line), and 21°C (dotted line). Experimental data by Scanley et al. (1990) () at 13°C and by Benndorf (1988) () at 21°C are also plotted.

The majority of the data presented to this point were used to determine the model's parameters. While it is important that the parameters adequately reproduce all the data used to determine them, it is also important that the parameters can be used to predict data not used in the fitting process. The ability of the model to fit data not used in determining the parameters is an independent test of how well the model approximates reality. In developing the model, ionic currents obtained with different voltage clamp protocols were used. The model was thus constructed so that it could reproduce the ensemble behavior of many sodium channels. In testing the model, therefore, measures of single channel behavior were chosen to see if the model could represent one channel as well as the average of many channels.

Fig. 12 shows the first latency densities for clamp potentials of 50 mV and 15 mV at 13°C and 21°C. Single channel simulations are done as described previously. At 13°C, for a clamp voltage of 50 mV, there is a wide variation of first latencies with many longer than 10 ms. In contrast, for a clamp voltage of 15 mV at 13°C, almost all of the channels have first latencies of <5 ms. The probability that a channel first opens after time t was computed from these histograms and plotted versus time. The plots have a nonzero plateau that describes the probability of not opening during the channel simulation as well as a single time constant with which the probability relaxes to this plateau value. The plateau and time constant values are 0.34 and 6.41 ms for 50 mV and 0.22 and 1.65 ms for 15 mV, respectively. At 21°C, for a clamp voltage of 50 mV, a much larger fraction of channels have first latencies <5 ms, although there is still a wide variation in latencies. For 15 mV, almost all of the latencies are <2 ms. The plateau and time constant values at 21°C are 0.76 and 3.69 ms for 50 mV and 0.52 and 0.63 ms for 15 mV, respectively. Experimental data at 21°C for a clamp potential of 50 mV yield a time constant of 1.15 ms (Berman et al., 1989). The corresponding plateau value is not available. The model probably has a larger time constant because of the wide variation in the latencies. The long latencies are due to channels inactivating from a closed state for a considerable time and then returning to a closed state from which the channel can then open. The long latencies might be eliminated by adding voltage dependence to the transitions between the closed and closed-inactivated states. However, as discussed previously, such a change would greatly increase the complexity of the model. Although experimental data with which to compare the model's data are limited, first latencies are related to the time to peak current and thus, like these values, should decrease with temperature at all voltages, as the model predicts.

View larger version (24K): [in this window] [in a new window]   FIGURE 12   Densities of single channel first latencies. Twelve hundred channels are simulated for 40 ms as described in the text. Bin size is 0.5 ms. (A) At 13°C for 50 mV. (B) At 13°C for 15 mV. (C) At 21°C for 50 mV. (D) At 21°C for 15 mV.

Finally, the model is used to predict the fraction of channels that do not open and the number of channels that reopen during a voltage clamp of 40 ms as additional tests of the model's ability to predict single channel behavior. Fig. 13 A shows the probability of not opening as a function of voltage at 13°C and 21°C. The probability of not opening is high at very negative potentials, while it is much lower at depolarized potentials. Even at these depolarized potentials, though, there is still a significant fraction of channels that do not open at both temperatures. The model's predictions for 13°C are generally lower than the experimental data reported by Scanley et al. (1990), although the overall trend of the curve is similar. The tendency for the model to predict slightly lower probabilities than those measured experimentally may be due to brief and missed openings in the experimental data. Experimentally, openings shorter than 178 µs cannot be detected (Scanley et al., 1990). If these openings are excluded from the model simulation, the resulting probabilities of a null sweep are increased by 5-20% depending on the clamp potential. The probability a channel does not open increases as the temperature is increased, as expected from data in the literature (Correa et al., 1992). This result suggests that as the temperature is increased, a larger fraction of channels are inactivated before they reach the open state.

View larger version (24K): [in this window] [in a new window]   FIGURE 13   Single channel probabilities for experimental data () at 13°C (Scanley et al., 1990), model data at 13°C (), and model data at 21°C (). (A) The probability of not opening versus clamp voltage. (B) The number of channel openings, normalized by the number of channels that open, versus clamp voltage.