Kinetic Model of Nav1.5 Channel Provides a Subtle Insight into Slow Inactivation Associated Excitability in Cardiac Cells

Voltage-gated sodium channel Nav1.5 has been linked to the cardiac cell excitability and a variety of arrhythmic syndromes including long QT, Brugada, and conduction abnormalities. Nav1.5 exhibits a slow inactivation, corresponding to a duration-dependent bi-exponential recovery, which is often associated with various arrhythmia syndromes. However, the gating mechanism of Nav1.5 and the physiological role of slow inactivation in cardiac cells remain elusive. Here a 12-state two-step inactivation Markov model was successfully developed to depict the gating kinetics of Nav1.5. This model can simulate the Nav1.5 channel in not only steady state processes, but also various transient processes. Compared with the simpler 8-state model, this 12-state model is well-behaved in simulating and explaining the processes of slow inactivation and slow recovery. This model provides a good framework for further studying the gating mechanism and physiological role of sodium channel in excitable cells.


Introduction
Nav1.5, the cardiac isoform of the voltage-dependent sodium channel a subunit, is encoded by the human SCN5A gene which locates on chromosome 3p21 [1]. In the firing process of cardiac action potentials (APs), Nav1.5 channel carries a rapid inward Na + current (I Na ) in response to a depolarization and then goes into ''a fast inactivated state'' with a small fraction of remnant currents throughout the whole process of depolarization. This persistent current termed as a late I Na plays a critical role in the heart [2][3][4]. After that, inactivated channels recover from the process of repolarization and prepare for the next process of depolarization. The generalized morphology of myocardial APs in humans and other large mammalian species contains a rapid upstroke followed by a depolarized plateau potential lasting for more than 100 ms [5,6], whereas the AP in neurons usually consists of a rapid upstroke and an immediate repolarization [7]. With the prolonged depolarization, Nav1.5 channels progressively enter ''a slow inactivated state'', corresponding to a slow recovery process with the time constants ranging from hundreds of milliseconds to several seconds [8,9]. Slow inactivation substantially suppresses Na + currents to control the cell excitability. Mutants which cause enhanced slow inactivation are often associated with several clinical heart diseases [10][11][12][13]. However, the gating mechanism of slow inactivation of Nav1.5 remains elusive.
Over the past 15 years, numerous mutations in SCN5A have been reported to be associated with various rare arrhythmia syndromes, such as congenital Long QT syndrome type 3 (LQTS3), Brugada syndrome (BrS), progressive cardiac conduction defect (PCCD), sick sinus syndrome (SSS) and arterial standstill [10][11][12][13][14]. To better understand the linkage between the gating of the Nav1.5 channels and heart diseases, we think that one feasible way is to use kinetic models. The dominant paradigm for modeling voltage-gated ion channel kinetics over the past 60 years has been dependent on the giant squid axons experiments of Hodgkin and Huxley [15]. Since then, the H-H models have been extensively used in data analysis of cellular electrophysiology. However, with the availability of high resolution data, many ion channels exhibit features beyond the traditional H-H models [16,17]. As a consequence, more complicated Markov models have been proposed for analyzing the ion-channel kinetics [18][19][20][21]. Such models produce more precise description to the ionchannel kinetics, which can be ultimately used to understand the firing properties of APs in excitable cells.
In this study, we proposed a two-step inactivation Markov model for simulating the Nav1.5 currents including the slow inactivation and bi-exponential recovery.
This work provided a solid basis for studying the detailed gating mechanism and the electrophysiological role of sodium channel in excitable cells.

Electrophysiology
The whole-cell mode was only used in all experiments. Patch pipettes were pulled from borosilicate glass capillaries with a resistance of 1.5-2.5 MV, after filled with pipette solution. The series resistances were compensated with 80%-90%. All the experiments were performed with a patch clamp amplifier (Axopatch 200B, Axon Instruments, Union City, Calif., USA) with its software (Clampex) at room temperature (23-25uC). The currents were typically digitized at 10-100 KHz and filtered at 5 kHz (Fig. S1)

Data Analysis
Patch clamp recording data were analyzed with Clampfit (Axon Instruments, Inc.) and Sigmaplot (SPSS, Inc.) software. Unless otherwise stated, the data are presented as mean 6 S.D..
The conductance-voltage (G-V) curves for activation were fitted to the a Boltzmann equation as below, where V 50 is the half maximal voltage, G the conductance, G max the maximum conductance and k the slope factor. The steady-state inactivation was fitted to a Boltzmann equation as below, where V 50 is a half availability voltage, and k is the slope factor. Development of slow inactivation was fitted to a single exponential function as below, where t is time constant, a and b are the partition coefficients. Recovery curves were fitted to the mono-exponentail or biexponential equations as below, where I is the peak current, I max the maximal peak current, A 1 and A 2 the proportional coefficients, t the time, t 1 and t 2 the fast and slow recovery time constants, respectively.

Mathematical modeling and simulation
The differential equations for the kinetic modeling were solved numerically, using a QMatrix or Five-order Runge-Kutta integration method. The fitting procedure is bad on a PSO-GSS algorithm for direct estimation of rate constants from macroscopic currents [22]. The integrating routines were written and executed with software CeL (HUST, Wuhan, Hubei, China), compiled with the C++ compiler to run under Windows XP [23]. Kinetic parameters were optimized with CeL as previously described. [22] Results

Models for Nav1.5 channels
To understand the gating mechanism of Nav1.5, the currents of activation, steady-state inactivation and deactivation are absolutely necessary to provide us the detailed information of rate constants required for a kinetic model. In this study, we are to provide a detailed description on the behavior of Nav1.5 currents arising from expression of SCN5A a-subunits in HEK293 cells, to present a kinetic model that appears to account for the observed currents of Nav1.5, and finally to construct a cardiac model cell with the built-in Nav1.5 channels to evaluate its physiological role in cardiac cells.
In models, the activation current depends on the forward (activation) rates from the closed (C) to the open (O) to the inactivated (I) states; the deactivation current depends on the backward (deactivation) rates from O to C; the steady-state inactivation current depends on the rates from O to I or C to I. Although currents rely on all the rates in model, each pathway is predominantly influenced by a certain combination of rates. Therefore, the model with all the rates can be finally determined from those currents [22].
Kinetic properties of Nav1.5 channels A set of the whole-cell experiments was thus performed on Nav1.5-transfected HEK293 cells for collecting the kinetics of activation (Fig. 1A), steady-state inactivation (Fig. 1B) and deactivation (Fig. 1C). In Fig. 1A, Nav1.5 currents exhibit a rapid activation and then a completed inactivation by a 20 ms depolarizing voltage steps ranging from290 to+60 mV in 5 mV increments after a holding potential of2120 mV to remove the possible inactivation. The activation and inactivation processes of Nav1.5 in Fig. 1A can be described by CRORI. Thus, the forward rates can be determined by fitting it to the activation currents. Deactivation currents were acquired by a 0.25 ms pulse to210 mV following a holding potential to2120 mV and ending with a 20 ms voltage steps from2100 to230 mV in 10 mV increments (Fig. 1C). Similarly, the backward rate of CrO can be determined by fitting it to deactivation currents. The voltage dependence of steady-state inactivation currents, standing for the availability of channels, was obtained by applying a set of conditioning voltages ranging from2120 to 0 mV for 500 ms in 10 mV increments and then measured at210 mV with a 20 ms test pulse as indicated in Fig. 1B. The backward rates of CrOrI can be determined by fitting them to the corresponding currents.
Alternatively, the channel availability can be obtained by directly measuring the recoveries of Nav1.5. In this study, a twopulse (prepulse P1 and test pulse P2) protocol was used for all of recovery experiments. Here we define the fractional currents as a ratio of I i (P2, V, t i )/I(P1, 2120 mV). Here I is the current, P1 the prepulse duration, P2 the test pulse duration, V the recovery voltage and t i the ith time interval between P1 and P2. Fig. 2A displays a set of fractional recovery curves arising from a two-pulse protocol (P1 = 30 ms, P2 = 20 ms, V = 2120,2110,2100 and290 mV). The recovery time constants are 5.160.9, 12.562.1, 26.163.8 and 47.963.4 ms for2120,2110,2100 and290 mV, respectively. A typical experiment exhibited a mono-exponential recovery (Fig. 2B). Especially, when t i = ', a set of the fractional ratio I i (P2, V, t i )/I(P1,2120 mV) will tend to their steady-state values, which can confer a steady-state inactivation curve. When P1 = 1000 ms, however, it shows a biexponential recovery (Fig. 2C), suggesting that a secondary (slow) inactivation state exists. Compared with the single-recovery curve (2120 mV), the double-recovery curve (2120 mV) is plotted in Fig. 2D. For P1 = 30 ms (empty circle), the recovery time constant is 5.0 ms; for P1 = 1000 ms (solid circle), the fast recovery time constant (t r-fast ) is 5.2 ms (78%) and the slow one (t r-slow ) is 596.3 ms (22%). To further explore the effect of P1 duration, a development of slow inactivation was executed. Cells were depolarized by a220 mV prepulse P1 with various durations to elicit inactivation, followed by a2120 mV interpulse with only a 30 ms duration, presumably to remove the recovery from fast inactivation (Fig. 2E). The remaining peak currents of the test pulse P2 measured at220 mV confers a time constant of 1.7960.11 s (Fig. 2F), suggesting that there is a slow inactivation component in Nav1.5 channels.

The kinetic model of Nav1.5 channels
To better understand the gating mechanism of Nav1.5, it is necessary to construct a Markov model for precisely matching all of the Na + currents shown in this study. After that, we aim to further explore the physiological role of Nav1.5 in cardiac cells.
At first, a typical 12-state model, composed of five closed states, six inactivated states and one open state, was used as the Nav1.5 model [25], However we failed to get a good fit to the present data. Alternatively, we considered the formalized Nav model of m 3 h proposed by Hodgkin and Huxley [15], where m is an activation factor and h an inactivation factor. Actually, m 3 h has an open probability equal to that of an 8-state Markov kinetic model [21]. Mimicking the Hodgkin-Huxley (H-H) model, an 8-state model, composed of three closed (C 1 -C 3 ) states, four inactivated (I 11 -I 14 ) states and one open (O) state, was chosen to simulate the Nav1.5 currents (Fig. 3A, Model I). Differing from the previous models, all rates are voltage-dependent in this model. The reason for this change was to make model more flexible and obtain the highest score evaluated by the software CeL [22]. For slow inactivation, we constructed a new model (Model II) by adding the four slow inactivated states I 21 -I 24 to I 11 -I 14 in Model I as shown in the boxed region (Fig. 3B). In this model, we assume that the occupancies of I 11 -I 14 transits into I 21 -I 24 so slowly as to lead a slow inactivation. The parameter values of two models are listed in Table 1.
In Fig. 4A-C, the Model I replicates the activation, deactivation and steady-state inactivation of Nav1.5, indicating that it is capable to be a Nav1.5 model. The time constants of data and  (Fig. 4D-F). The V 50 of activation is234.5 mV for data and234.0 mV for simulation (Fig. 4G); the V 50 of steady-state inactivation is289.1 mV and289.2 mV (Fig. 4H). Recovery for the short P1 of 30 ms could also be replicated by Model I (Fig. 4I, J). Time constants from2120 to290 mV are 5.1 and 5.6, 12.5 and 11.9, 26.1 and 26.6, 47.9 and 49.6 ms for data and simulation, respectively.
However, Model I failed to account for the slow recovery process. We realized that this insufficiency was related to the structural constraint of Model I. According to this 8-states model, it can not account for the bi-exponential recovery process, which has the time constants t r-fast = 5.2 ms and t r-slow = 596.3 ms with a nearly 100 times difference (Fig. 2). Interestingly, the time constant (t r-fast ) of fast component remains no change, compared with that of single recovery, suggesting that it is better to insert a secondary (slow) inactivation pathway to model I. The simplest way is obviously to directly add one more inactivated state I s (or I slow ) to open state (O) [7,26]. Essentially, this I s is to share a certain amount of probability with the original I f (or I fast ) state. It will slowly go and back between I s and O states as a slow component of recovery. Unfortunately, this modified model was incapable to account for the duration-dependent inactivation or bi-exponential recovery in our trials. After trying several different models, we finally determined to use the 12-state two-step inactivation model named as Model II (Fig. 3B). Based on the structure of Model II, the slow inactivation of channels should undergo two steps during the depolarization process, for simplicity, first from O to I 14 (fast) and then to I 24 (slow). Upon repolarization, the fractional occupancy residing in the fast inactivation state I 14 transfers quickly from I 14 to C 1 mainly via I 13 , I 12 and I 11 , while the fraction residing in the slow inactivation state I 24 transfers slowly from I 24 to C 1 mainly via I 24 , I 23 , I 22 , I 21 and I 11 . Prolonging the duration of depolarization, more channels go to I 24 from I 14 and more channels have to recovery via a slow pathway. The rates between I 1x and I 2x are so small that Model II fitted nearly the same results to Model I in short-duration processes (Fig 4A, C-G, I-J). There is only a slight difference on fits, comparing the Model II with Model I (Fig 4B  and H). However, Model II was also well-behaved in simulating the long P1 duration case. More kinetics of Model II is shown in Fig 5. The development curve of slow inactivation is shown in Fig 5A. Time constants are 1.79 s for data and 1.58 s for simulation. Moreover, all the fits (red) generated by the Model II in Fig. 5B-C are coincident with the recovery traces (black) shown previously in Fig 2B-C, suggesting that the Model II can also well replicate the whole processes including the slow recovery. For the short P1 of 30 ms, recovery time constants are 5.0 ms and 5.6 ms for data and simulation, respectively; for the long P1 of 1000 ms, are 5.2 ms (fast, 78%) and 596.3 ms (slow, 22%) for data, 5.6 ms (fast, 77%) and 557.2 ms (slow, 23%) for simulation (Fig. 5D).  Each of rate constants has the expression x = k * exp V/n except the r 1 = g/ (1+exp(2(v+a)/f)) and the cyclic balancing rates b 3 , calculated from microscopic reversibility of cycles. Here the rate x and the pre-exponential factors k and g are in ms 21 and the exponential factor n, a, f and voltage V (trans-membrane voltage) in mV 21 . The parameter c is non-dimensional. doi:10.1371/journal.pone.0064286.t001

Discussion
In this study, a novel 12-state two-step inactivation kinetic model was successfully developed to study the slow inactivation mechanism of Nav1.5 channels, especially to illustrate the mechanism of duration-dependent bi-exponential recovery based on a set of Nav1.5 currents recorded by patch-clamp experiments.
The proposed model II (or two-step inactivation model) can be regarded as a hierarchical structure. Firstly, Model I well duplicated the activation, deactivation, steady-state (fast) inactivation and fast recovery (Fig. 4). Three closed states in Model I correspond to the m 3 Hodgkin-Huxley formalism, previously suggested by Markov models [18], and FRET experiments [27,28]. Especially, the latter suggested that only three major protein conformations were required for ion permeation to occur. Secondly, Model II well described all the Nav1.5 kinetics including the slow inactivation and slow recovery, the mission impossible for Model I (Fig. 5).
The slow inactivation process in Model II can be briefly summarized as OrRI 1 rRI 2 , named as two-step inactivation model. Here I 1 and I 2 are the fast and slow inactivation states, respectively. In contrast, Jarecki et al. used an alternative two-state inactivation model I 2 rROrRI 1 to investigate the resurgentsodium currents [29]. We tested this model and found that it did not work well in explaining both the inactivation and the slow recovery at variety of voltages. Therefore, the two-step inactivation model should be more appropriate to our work.
HEK293 cell has endogenous voltage-gated Na channels such as Nav1.7 and voltage-gated Ca2+ channels [30]. But the endogenous currents of HKE293 in different labs were variable because of different generation, environment, or other factors. We used cesium instead of potassium in our pipette solution (see in Method) in order to block the endogenous potassium current of HEK293. Secondly, we used 100 nM TTX as blockers (Fig. S2A). The change of Nav1.5 current was very small, which means there was no significant Nav1.7 current (TTX sensitive) in our HEK293 cells. Then we tested the currents of un-transfected cells. A large proportion of results revealed that there was no significant endogenous currents. Other proportion of results showed endogenous inward current with amplitude of 20 to 40 pA (Fig. S2B), which was less than one percent of transfected Nav1.5 current described in our work.
It is important for this discussion to point out that sodium channel (including a and b subunit) functions are critically dependent on the particular heterologous expression system used [31]. Current data sometimes are conflicting, possibly due to differences in experimental conditions or, more importantly, species studied [32]. In Xenopus oocytes, for instance, injection of the a subunit cRNA shows an abnormally large component of a intermediate inactivation mode [33][34][35], which is substantially reduced by co-injection of cRNA encoding the b1 subunit [34]. In contrast, the alone expressed a subunit in HEK293 has shown a kinetic behaviour similar to that of the native preparations [36,37]. This is due to an abundant endogenous expression of mRNA encoding the b1A subunit in HEK293 [38]. The endogenous b1A subunit is sufficient for suppressing the intermediate inactivation of sodium currents by co-assembly with a-subunits [38]. Though the kinetics of our Nav1.5 model may have slight difference compared with other results obtained from variant tissues or species, the mechanism of slow inactivation should be the similar. Our model can be also applicable to the data from other tissues or species with proper adjustment in parameters. As mentioned previously, some mutations associated with heart diseases, e.g. BrS or PCCD, often occur from the enhanced slow inactivation of Nav1.5 [10,[39][40][41][42], as it suppresses the Na + current (loss of function) to decrease the excitability of cardiac cells. In other words, a loss of Na + current prevents the membrane potential from reaching the threshold of AP, thereby slowing conduction. Therefore, our work provides a useful tool to investigate the linkage between the slow inactivation of Nav1.5 and clinical heart diseases.
Considering that there is a high similarity in the Nav1.x family, we believe that the two-step inactivation model is applicable to the other sodium channels with their mutations, after making proper changes in parameters. All of quantitative analyses were made possible by recent advances in kinetic modeling algorithms and software (CeL) that allowed us to quickly construct and fit models with large amounts of states from the comprehensive data ranged from 5 ms to more than 10 s [22]. Therefore, this work provides a convenient platform for further investigating the detailed linkage between the sodium channels and diseases in the future.