SIMULATION AND CALCULATION OF THE CONTRIBUTION OF HYPERPOLARIZATION-ACTIVATED CYCLIC NUCLEOTIDE-GATED CHANNELS TO ACTION POTENTIALS

The hyperpolarization-activated cyclic nucleotide-gated (HCN) channel, which mediates the influx of cations, has an important role in action potential generation. In this article, we describe the contribution of the HCN channel to action potential generation. We simulated several common ion channels in neuron membranes based on data from rat dorsal root ganglion cells and modeled the action potential. The ion channel models employed in this paper were based on the Markov model. After modifying and calibrating these models, we compared the simulated action potential curves under the presence and absence of an HCN channel and calculated that the proportional contribution of the HCN channel in the potential recovery phase was 33.39%. This result indicates that the HCN channel is critical in assisting membrane potential recovery from a hyperpolarized state to a resting state. Furthermore, we showed how the HCN channel modifies the firing of the action potential using mathematic modeling. Our results indicated that although the loss of the HCN channel made recovery of the membrane potential more difficult from the most negative point to resting in comparison with the control, the firing rate of the action potential increased in certain circumstances. We present a novel explanation for the HCN channels’ mechanism in neuron action potential generation using mathematical models.


INTRODUCTION
The action potential in a dorsal root ganglion cell is formed under the contributions of a number of ion channels, and these ion channels seek a dynamic equilibrium, which forms a sine-wave-like signal (Schulz, 2006;Turrigiano, 2008).Hyperpolarization-activated cyclic nucleotide-gated cation (HCN) channels are open when the membrane potential becomes hyperpolarized.HCN channels help the membrane potential recover from the hyperpolarized phase to resting as a "pacemaker", and the neuron easily accepts another stimulation and cycles through action potential firing (Magee, 1998;Wahl-Schott and Biel, 2009).However, it is complicated to measure the precise con-tribution of HCN channels using a patch clamp due to the mixed channel signal in a living cell.Here, we used a mathematical model to calculate the contribution of HCN channels to the action potential in dorsal root ganglion (DRG) cells.Mathematical simulations are very powerful in research as they can help us to avoid the defects of in vivo experiments.Also, researchers can use simulations to predict possible results and to create a reference for upcoming work.Our results showed that HCN channels provided approximately 1/3 of inward nerve currents.In addition, by applying this model we found a reasonable explanation for the priming effect of the downregulation of HCN expression in action potential generation.

Animals and electrophysiology
Animal studies conformed to the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health and were approved by the Institutional Animal Care and Use Committee at the South-Central University for Nationalities.

Simulation and calculation
Cell simulation software (Cel) was provided by Prof. Ding's group (Institute of Biochemistry and Biophysics, Huazhong University of Science and Technology) to for all simulation processes (Wang et al., 2012).The Cel software has three main modules: (i) voltage clamp for kinetic simulation of ion channels; (ii) current clamp for simulating the action potentials of cells, and (iii) fit-kinetics for automatically calculating the parameters of kinetic models.
To find the contribution of the HCN channels, we set up models of potassium channels, sodium channels, Ito (the transient outward K + current) channels and HCN channels in the Cel software.To simulate the type of channel, we began with a 2-state model, including one open form and then we expanded it to more open form, closed form and inactivation form.The expansion included two parts: the first one was to search for the appropriate number of states, and the second was to calculate the proper transformation probability among the different states.The initial transformation parameter value was set to 1; we used the Cel software to adjust this value within a range of 0.01-100 to simulate the ion channel current.We then input the actual experimental data into the Cel software and fitted the data with the parameters to test whether they were appropriate.While establishing our model, we used an MWC/KNF model to describe each ion channel.This model can only be used on the cells with soma only, because all the data were gathered from DRG cells without synapses.After models of each ion channel were built, we combined them to generate an action potential wave using the Cel software.Then, we contrasted the parameters between the artificial action potential and the recorded action potential.According to these results, we adjusted our model to make the generated curve more coincident to the real curve obtained by patch clamping.We drew the simulated and the real action potential curve in one graph and then calculated the area between these two curves.The smaller the area means a better simulation (Cui and Mao, 2012).Supplementary Fig. 2 shows the simulated and the recorded current curves.In general, the critical parameters were the local conductance, reversal potential and capacitance.The minimization of difference between recorded and modeled AP was the goal of our process.We selected mean absolute error and mean relative error as evaluation standards (Cui and Mao, 2012).
The mean absolute error (MAE) describes the average deviation between simulated and true values.The mean relative error (MRE) describes the average deviation degree between them, and can better reflect the advantages and disadvantages of the model.
The expressions of the MAE and the MRE are as follows: Obviously, if MAE and MRE are smaller, then the error is smaller, and the simulation accuracy higher.We performed the computation in our model and the MAE and MRE are presented as follows: which shows that the value of the MAE and MRE is in a reasonable range.The value of MRE is small enough for us to believe that the simulated curve is acceptable (Supplementary Fig. 1 shows the calculated area between curves).

Data analysis
Electrophysiology data were recorded by patch clamp with voltage clamp or current clamp.Data were analyzed by Clampfit (V 10.0, Molecular Devices), Sig-maPlot 10.0 and GraphPad Prism 4. Current tracing was approximated using the Runge-Kutta method.

Simulation of ion channel model
As shown in Fig. 1a, the overall model of the membrane channel was established.Here, we only considered the most important ion channels for action potential production, including potassium channels, sodium channels, Ito channels and, our main focus, HCN channels.Leakage was characterized as a consistently open channel on the cell membrane.In this model, two global parameters were determined first: the capacitance of the cell membrane and the resting potential (Table 1a).The mimic stimulation protocol is shown in Fig. 1b, for which the sampling frequency was 10 KHz.
In the process of simulating the action potential, it was necessary to adjust the local conductance and reverse potential of each ion channel model to bring the simulated curve closer to the real data.The final parameters used in this article are shown in Table 1b.

Sodium channel model
The sodium channel model was based on an existing model (Xiong, 2009) and was modified to make it more suitable to the actual data.We assumed that a sodium channel had 5 closed states, 5 inactivated states, an open state and a blocking state (Kuo and Bean, 1994;Raman and Bean, 2001).The rate of transition between states of similar types was proportional to sodium channel current (Douglas et al., 1967) (Fig. 2-a).The values of the parameters in Fig. 2a are shown in Table 2a.

Potassium channel model
Similar to the method of simulating the sodium channels and in accordance with the characteristics (Conforti and Millhorn, 1997;Thompson and Nurse, 1998;Conforti et al., 2000;Fearon et al., 2002) of potassium channels, we defined a model for the potassium channels.This model consisted of 5 states, including 4 closed states and 1 open state.The values of the parameters in the expression shown above and the figures are presented in Table 2b.

Transient outward potassium channel model
For transient outward potassium channels, we used a kinetic model described by Patel and Campbell (2005), from whom we obtained the specific models and parameters (Fig. 2c, Table 2c).

HCN channel model
With respect to kinetics, HCN channels and potassium channels have relatively high similarity.Thus, we first built the HCN channel model with reference to the potassium channel model.According to the characteristics of HCN channels, we designed a 5-state model containing 4 closed states and 1 open state (Fig. 2b).This model reflects the kinetics of HCN channels suitably because the simulated curve can fit the actual data well.
We then tried the model described by Altomare et al. (2001) and obtained a better simulation result.This model reflected the different states of the tetramer structure of HCN channels, including 5 closed states and 5 open states (Fig. 2d).After fitting the experimental data, we obtained the specific parameters (Table 2d).Analyzing these parameters, we found that h1 and i4 had the largest conversion probabilities among all the data, which were several times or even tens of thousands of times the values of other parameters.Hence, the conformations of the HCN channel tetramer were mainly restricted in C4 and O1 states in our model.

Leakage
Leaks were modeled as consistently open channels; thus, it was not necessary to build a specific leak model.

Simulation of action potential
As previously noted, we successfully simulated sodium channels, potassium channels, Ito channels and HCN channels and acquired several critical parameters.We then input these parameters into the Cel software to generate an action potential curve.After adjustment, we obtained a simulated action potential close to the real one.This simulated curve was suitable for testing the electrophysiological function of HCN channels.For example, we could selectively remove the HCN channel model to identify its effect on action potential generation.Comparing the result with experiment data, we could determine whether our model and parameters were reasonable.We analyzed the action potential data generated from adult rat DRG cells.Under the same current injection stimulation (300pA), blocking of HCN channels increased the action potential firing number (Fig. 3c).This result demonstrated that the inhibition of HCN channels facilitated the generation of nerve signals

HYPERPOLARIZATION-ACTIVATED CYCLIC NUCLEOTIDE-GATED CHANNELS AND ACTION POTENTIALS
and was consistent with observations in epileptic animals (Chung et al., 2009).Although this phenomenon was difficult to explain on the cellular level, our model successfully simulated it.As shown in Fig. 3c in the model without HCN channels, we simulated a 300pA current injection for 150 ms, and obtained 4 action potentials.After we set the contribution of HCN channels to 0, the number of action potentials increased to 6.Although the net influx of cations was reduced, the action potential could be elicited by a lower potential (Fig. 3b).

The role of HCN channels in action potential hyperpolarization recovery
HCN channels are activated when the membrane potential is hyperpolarized, but their proportional contribution from the hyperpolarized phase to the resting phase remains unclear.Here, we used the slope of the recovery in the simulated curve to measure the HCN channels' contribution.If the slope reached 0 when the HCN channels were removed from the simulation model, then the contribution of the HCN channels to recovery is 100%.However, if the slope did not change after the HCN channels were removed, then the contribution of the HCN channels is 0%.As shown in Fig. 3b, removing or keeping HCN channels in our model changed the slope of recovery.Therefore, we obtained the following equation: In this equation, K is the slope of recovery; V rest is the resting potential; V min is the potential of the most negative point; ΔT is the time from the most negative point to the resting potential of the curve.Thus, we obtained slope K m from the control simulated curve and K mwHCN from the model curve in which the HCN channels were removed.The values of K m and K mwHCN were 1.12 and 0.746, respectively.Then, we calculated the contribution of the HCN channels from the slope ratio (G HCN ): 0.746 ( 1) 100% (1 ) 100% 33.39% 1.12

DISCUSSION
HCN channels carry a mixed sodium-potassium ion current inward and play an important role in excitable tissues, such as cardiac and brain.HCN channels are activated by membrane hyperpolarization.In the generation of an action potential firing, HCN channels lead an inward membrane current following hyperpolarization, depolarizing the membrane to the threshold for the next action potential (Brown et al., 1979;Brown and Difrancesco, 1980).However, the exact contribution of HCN channels in this depolarization progress is an open question for many reasons, such as experimental data uncertainty, interruption of the specific channel inhibitors of other channels and cell states.In this investigation, we found the slope of the curve in a mathematical model and found the contribution of the HCN channels to be 33.39%,which could be used as a reference in future experiments.
Due to their function in action potential generation, HCN channels are related to many neuron system diseases.HCN channels are key determinants of neuronal network activity (Baruscotti et al., 2010) and related to nerve system diseases, such as neuropathic pain (Emery et al., 2012), inflammatory pain (Emery et al., 2011) and epilepsy (Kase et al., 2012).Many HCN blockers or openers are being screened to cure such symptoms and one of them is already used (Postea and Biel, 2011).But there is a contradiction about the function of HCN channels.Generally speaking, at the cellular level, blocking HCN channels results in less cation flow into the cell, which makes the cell more difficult to excite (He et al., 2014).This mechanism is sufficient to explain the effect of HCN channel inhibition on neuropathic pain (Wells et al., 2007;Weng et al., 2012), but this role is contradicted by the currently accepted mechanism of epilepsy.Knocking out HCN channels from mice does not reduce the risk of seizure but causes the mice to be more prone to epilepsy (Kole et al., 2007;Kase et al., 2012).Hence, the effect of the regulation of HCN channels on cell excitability is distinct in different circumstances.Our model shows one such possibility: that the inhibition of HCN channels could make the cell more excitable and facilitate the firing of action potentials (Fig. 3).
In short, we successfully simulated action potential curves in our model, and we estimated the contribution of the action potential to hyperpolarization.We have also explained the function of HCN channels in regulating cell excitability.This model contains only sodium channels, potassium channels, Ito channels, leaks and HCN channels, which may lead to some deviation and should be expanded in future work.2.

Fig. 1 .Fig. 2 .
Fig. 1.Models and simulated protocol.a -model of membrane channels; b -the protocol used to generate the simulated action potential.

Fig. 3 .
Fig. 3. Simulation models for Ion channels.a -simulation model of a sodium channel.C i (i = 0,1,2,3,4), closed state of the sodium ion channel; O, open state of the sodium channel; I i (i = 0,1,2,3,4), inactivated state of the sodium channel; O B , blocked state of the sodium channel; α,β,λ,λ 4o ,θ,ξ u ( u= 0, 1, 2, 3, 4), ζ u ( u = 0, 1, 2, 3, 4),r, s, j, n, w,, the probability of a conversion factor; b -potassium channel model.C i (i = 0,1,2, 3), closed state of the potassium channel; O, open state of the potassium channel.a, base of the potassium channel transition probability from the closed state to the open state.b, base of the potassium channel transition probability from the open state to the closed state; c -transient outward potassium channel model.C i (i = 0,1,2, 3), closed state of the transient outward potassium channel; I i (i = 0, 1, 2, 3, 4) and I, inactivated state of the transient outward potassium channel; O, open state of the transient outward potassium channel; α, β, θ, λ, μ, ω, rate of transition between different states of the transient outward potassium channel; d -HCN channel model.C i (i = 0, 1, 2, 3, 4), closed state of the HCN channel; O i (i = 0, 1, 2, 3, 4), e, open state of the HCN channel.α, β, θ, λ, rate of transition between similar states of the HCN channel; h i (i = 0, 1, 2, 3, 4), rate of transition from the closed state to the open state of the HCN channel; i i (i = 0, 1, 2, 3, 4), rate of transition between from various open states to various closed states of the HCN ion channel.ν, variable in the model.The values of other parameters are given in Table2.

Table 2 .
Parameter values of expressions in Fig.2a-d.