The Importance of the Dissociation Rate in Ion Channel Blocking.

Understanding the relationships between the rates and dynamics of current wave forms under voltage clamp conditions is essential for understanding phenomena such as state-dependence and use-dependence, which are fundamental for the action of drugs used as anti-epileptics, anti-arrhythmics, and anesthetics. In the present study, we mathematically analyze models of blocking mechanisms. In previous experimental studies of potassium channels we have shown that the effect of local anesthetics can be explained by binding to channels in the open state. We therefore here examine models that describe the effect of a blocking drug that binds to a non-inactivating channel in its open state. Such binding induces an inactivation-like current decay at higher potential steps. The amplitude of the induced peak depends on voltage and concentration of blocking drug. In the present study, using analytical methods, we (i) derive a criterion for the existence of a peak in the open probability time evolution for a model with an arbitrary number of closed states, (ii) derive formula for the relative height of the peak amplitude, and (iii) determine the voltage dependence of the relative peak height. Two findings are apparent: (1) the dissociation (unbinding) rate constant is important for the existence of a peak in the current waveform, while the association (binding) rate constant is not, and (2) for a peak to exist it suffices that the dissociation rate must be smaller than the absolute value of all eigenvalues to the kinetic matrix describing the model.

Introduction

Understanding the relationships between the rates and dynamics of current wave forms under voltage clamp conditions is essential for understanding phenomena such as state-dependence and use-dependence, which are fundamental for the action of drugs used as anti-epileptics, anti-arrhythmics, and anesthetics. In the present study, we mathematically analyze models of open state blocking mechanisms previously suggested for the local anesthetic bupivacaine action on Kv channels (Longobardo et al., 2001; Nilsson et al., 2003, 2008).

The dynamics of ion channels are generally considered to be memory-less (i.e., they possess the Markov property) and can be analyzed in terms of Markov chains (Colquhoun and Hawkes, 1995). We, therefore, explore Markov-chain type kinetic schemes, describing open state dependent drug-binding. Both analytical and Monte Carlo methods are used. We derive criteria for the existence of an induced current peak and formula for the peak height and its dependence on membrane potential. It should be noted that similar serial Markov chains can be used to describe a number of other pharmacological processes, such as competitive binding of antagonists or agonists, and, consequently, insights from the present analysis can also be of value in these cases.

Two findings are apparent from this study: (i) the dissociation (unbinding) rate constant is important for the existence of a peak in the current waveform, while the association (binding) rate constant is not, and (ii) for a peak to exist it suffices that the dissociation rate must be smaller than the absolute value of all eigenvalues to the kinetic matrix describing the model.

That the criterion is independent of the association rate (as long as it is greater than zero) has the implication that the concentration of the local anesthetic does not influence the existence of a peaked waveform under voltage clamp.

Methods

All ion channels, either voltage-gated or ligand-gated, open and close randomly, and an accurate explanation of their behavior must, therefore, be of a probabilistic character (Colquhoun and Hawkes, 1995; Johnston and Wu, 1995). In this work, we use analytical and Monte Carlo methods to analyze Markov chain models of ion channels, describing the probability of the channel being in each state. We investigate the following three-state Markov scheme using analytical methods

where C, O, and B denote closed, open and blocked states, respectively, where α, β, γ, and δ denote rate constants, and where γ includes the concentration of the blocking drug. We used the terms κ x L and λ for γ and δ, respectively. We also analyze extended versions of Scheme 1, namely

This extended kinetic scheme describes a channel system with m equal and independent gates in accordance with the Hodgkin-Huxley formalism for Kv channels (as first noted by Richard Fitzhugh, 1961). With this notation, the number of states in the scheme is n = m + 2.

Results

The existence of a current peak during voltage clamp: analysing a three-state scheme

Time-dependent effects of bupivacaine on different Kv channels at + 60 mV at different concentrations have been analyzed in several studies, using voltage clamp technique (see Gonzalez et al., 2001; Nilsson et al., 2003). As can be inferred from these studies, a current peak can be observed at higher voltages. The voltage dependence of this induced peak, however, has been surprisingly little studied (but see Gonzalez et al., 2001).

We seek to find sufficient and necessary prerequisites for the existence of such a peak. We will first derive this for a three-state scheme, then use similar techniques to derive criteria for the four-state scheme, and in the Appendix for an n-state scheme. The techniques are generally those of dynamical systems and matrix algebra, which we (Gouwens et al., 2010; Zeberg et al., 2010, 2015; Sahlholm et al., 2016) and others (for an overview see Koch, 2004; Izhikevich, 2007) have used extensively for analyzing these systems. The differential equation for Scheme 1 written in matrix form is x′ = Ax where

with the general solution

and where V and r1 are the eigenvectors and eigenvalues of the transition matrix A. c are constants dependent on the initial conditions. Experimentally, voltage clamp experiments are typically done with a strongly negative potential as an initial condition, meaning that all ion channels are in the (first) closed state, i.e.,

We will now use the Putzer algorithm to solve O(t). The Putzer Algorithm is a method for analytically evaluating matrix exponentials using only eigenvalues and components in the solution of a relatively simple linear system (Putzer, 1966). This approach might seem a little bit cumbersome for solving the three-state scheme but will later enable us to solve not just the four-state model but also the general case (complete proof given in the Appendix). By the Putzer algorithm, the solution to x′(t) = Ax(t), where A is a 3 × 3 matrix, can be written on the form

Where p(t) and M1 are defined as follows. Define 3 × 3 matrices M1, M2, and M3 by the formula

and let the functions p1(t), p2(t), and p3(t) be given by solutions to the differential system

Note that here the eigenvalues can be in any given order, and we are not assuming that they are ordered in any particular way. We will now investigate p3. Solving the equation system above reveals that

Since is only to be found in p3(t) we have that

To get the pre-exponential factor c3V3 we write

The result above will be a 3 × 1 vector. Since we are solving for O(t), we are only interested in the second element of this vector. After matrix multiplication, the pre-exponential factor to in the expression for O(t) is

Since the sum of the eigenvalues of a matrix equals the trace, the expression be further reduced as in Equation (12). One of the eigenvalues is always zero so let r1 = 0, and assume without loss of generality that r3 > r2. r3 is then the slowest decaying term (note that all eigenvalues are non-positive). For some positive value of t the following inequalities must hold

why c3V3 will be the dominant term as t → ∞. Thus, a peak will exist for O(t) if the slowest decaying term r3 has a positive pre-exponential factor, i.e.,

By the assumption that r1 > r3 > r2 and that α is always positive it follows that this is true exactly when δ + r3 < 0. Using Vieta's relationships in conjunction with the border condition, r3 = −δ we can reduce δ + r3 < 0 to a simpler form. This is due to the fact that there are three equations of the Vieta relationships for det(A−rI) = 0 if A is a 3 × 3 matrix. In conjunction with r3 = −δ we have four equations and three eigenvalues, why it is possible to fully eliminate the eigenvalues. Some straight forward algebra reveals that

Thus a sufficient condition for a peak to exists in a three-state scheme is that the opening rate constant α must be greater than the dissociation rate constant δ.

The existence of a current peak during voltage clamp for a four-state and n-state scheme

For the four-state model the solution to p4(t) is obtained by solving the natural extended version of Equation (8), and by the Putzer algorithm we have

To get the pre-exponential factor c4V4 we write

where U is some vector with a zero on the third row (due to the tri-diagonal nature of the matrix A, the number of non-zero diagonals increase with two for each multiplication and we can omit lower terms than A2). After matrix multiplication of Equation (17) and taking the third element of the resulting vector we obtain the pre-exponential factor to in the expression for O(t),

Again, the expression is reduced using the relationship between the sum of eigenvalues and the trace. Now we can assume without loss of generality that r4 > r3 > r2 and r1 = 0. As for the three-state model, a peak will exist in the waveform if the pre-exponential factor to the slowest decaying eigenvalue is positive. This is true for a model with an arbitrary number of closed states, since all eigenvalues are real. Then by a sign analysis, the pre-exponential factor to is positive if, and only if,

Again using the Vieta's relationships this can be reduced to

For the interested reader the solution to the general case is given in the Appendix. The technique used in the Appendix is the same as for the three and four state model. It transpires that the pre-exponential factor to is

Again, let r > r > ⋯ >r2 and let r1 = 0. A similar sign analysis can be performed, yielding that the pre-exponential factor to the slowest decaying eigenvalue is positive exactly when

If β = 0 this criterion is reduced to α > δ. We observe that the influence of β seems to increase with the number of closed states, whereas γ is not involved in the criterion for a peak. The appendix includes a proof that the association rate γ does not affect the criterion for a peak. Figure 1 shows regions associated with the existence of a peak in the β − δ plane, scaled for α.

Figure 1

Regions associated with a peak in the β−δ plane (α = 1). The existence of a peak was independent of the value of γ, as long as γ > 0, for all cases. (A) One closed state (α > δ). (B) Two closed states. (C) Three closed states. (D) Four closed states.

Peak height as a function of the rate constants

To analyze the role of the rate constants in the height of the peak, we introduce a factor ψ such that the peak probability O can be expressed as

where ψ is dependent on the rate constants and o is the steady state open probability in the presence of a blocking agent (see Nilsson et al., 2003). For the three-state model the peak will have its maximum at time

This formula was obtained by taking the time derivate of the solution for O(t) (see Appendix, A18) and equating this expression with zero. Taking O(t) yields

Using our analytical solution it is possible to analyze how extreme values of the rate constants affect the peak height. Extreme values are shown in Table 1.

Table 1

Limits for ψ for various extremes.

	i= ∞	i= 0
limα→iψ(α)	γδ	0
limβ→iψ(β)	0	{γδ(γα−δ)^γ+δα−γ−δ α≤γ+δγδ(α−δγ)^−γ+δα−γ−δ α≥γ+δ
limγ→iψ(γ)	αδ-1	0
limδ→iψ(δ)	0	∞

In contrast to the criterion for the existence of a peak current, which involves only two rate constants (i.e., α and δ), the peak height depends on all rate constants of Equation (1) (i.e., α, β, γ and δ). To facilitate the use of the derived relationships, we can choose simplifying conditions for the system. Thus, assume that we investigate the blocking effect at high voltage steps (i.e., assuming β = 0 and α > γ + δ) and at a concentration equal to the Kd-value of the blocking agent (i.e., assuming γ = δ). Algebraically, the following formula is obtained

Peak height as a function of the voltage

To analyze how the voltage affects the amplitude of the induced peak we used Monte Carlo simulations. Let α and β be described according to Eyring-Polanyi rate theory (Evans and Polanyi, 1935; Eyring, 1935) as

where k is the pre-exponential factor (or the characteristic frequency factor), V1/2 is the potential for which α = β and s is a slope factor describing the influence of the potential. Again, let the concentration of blocking agent be at Kd-value (i.e., γ = δ).

Using Monte Carlo simulations in conjunction with analytical tools, we found three different regions in the δ−V plane: (i) One area, which we call A1, in which there is no open probability peak; (ii) a second area, A2, in which there is a peak and the relative amplitude (peak open probability relative to steady state with no blocking agent) decreases with potential and (iii) a third area, A3, in which there is a peak, but the relative amplitude increases with potential, approaching a value of one. Area A1 was defined using the criteria obtained above, i.e. α > δ. The border between A2 and A3 was obtained in the following way. Using custom written software in Mathematica 11.0 (Wolfram Inc.), we evenly assigned values using the built-in function RandomReal in parameter space of the right half of Figure 2 and investigated the voltage dependence (see Equations 26–28). The border between A2 and A3 were drawn using the built-in function ContourPlot. Figure 2 shows the three topological regions.

Figure 2

Topologically equivalent regions in the rate-potential plane. The three-state model (Equation 1) at Kd-concentration of the blocking drug (i.e., γ = δ). A1 represents an area in which there is no open probability peak and the steady state value decreases with potential approaching a value of a half, A2 represents an area in which there is a peak and the relative amplitude decreases with potential and A3 represents an area in which there is a peak, but the relative amplitude increases with potential, approaching a value of one.

The findings in Figure 2 are to some extent congruent with the findings of the relatively few experimental investigations of local anesthetic effects on non-inactivating Kv channels (see e.g. the effect of bupivacaine on Kv1.5 channels, described by Gonzalez et al. (2001), where the relative peak amplitude goes from a potential region where there is no peak, i.e., corresponding to A1, and to a potential region where the relative peak increases with potential, i.e., corresponding to A3). Clearly this issue requires further experimental investigations to be clarified.

Discussion

The present analysis was undertaken to derive general principles of kinetics from the models of binding mechanisms. In previous studies, we showed that the action of the local anesthetic bupivacaine on a non-inactivating potassium channel could be described by Markov chain models, assuming that binding occurs exclusively to channels in open state (Nilsson et al., 2003, 2008).

In the present study, we mathematically analyzed Markov models, with special reference to the question of the existence of a peak and its voltage dependence. We (i) derived the criterion for the existence of a peak in the open probability time evolution for an open-state binding kinetic scheme, comprising one and two closed states (Schemes 1 and 2.1), (ii) derived the criterion for a peak in an open-state kinetic scheme with an arbitrary number of closed states, (iii) derived formula for the relative height and the block of the peak amplitude for Scheme 1, and (iv) determined (by Monte Carlo simulations) the voltage dependence of the relative peak block for Scheme 1.

Contemplating these findings, we note the important role that the dissociation rate constant plays for the open probability peak features. Intuitively, one would have expected that the on rate (γ) would be the key player in the existence of a peak. Nevertheless, we find that the association (binding) rate constant is the only rate constant that does not influence the existence of a peak. In a three-state model (Scheme 1), the existence of a peak is given by the simple relation α > δ (i.e., the activation rate constant is greater than the dissociation (unbinding) rate constant). If β = 0 (the deactivation rate constant), this relationship holds for models with a higher number of closed states. Generally, a peak exists if −r > δ where r is the eigenvalues closest to zero (i.e. the slowest decaying term).

Additionally, we could determine a topological region in the dissociation rate–voltage plane for Scheme 1 that is characterized by a decreasing block of the induced peak with potential.

Using analytical and Markov chain models that describe the action of blocking drugs on ion channels, we could derive general principles of linear three- and higher-state schemes. Thus, we could show that the existence of a current peak for open state binding schemes mainly depends on the dissociation (unbinding) rate constant δ, the criterion for a scheme with one closed state (Scheme 1) being α > δ (i.e., the activation rate constant should be greater than the dissociation (unbinding) rate constant). We could also show that different peak amplitude-voltage relations characterize specific topological regions of the dissociation rate–voltage plane for Scheme 1, a finding that still awaits experimental corroboration.

Understanding the relationships between the rates and the dynamics of the block of ion channels is essential for understanding how more drugs modulate neuronal firing patterns and thus how they function in pain, epilepsy, arrhythmia and anesthesia; these insights are crucial when developing anti-epileptic, anti-arrhythmic and anesthetic drugs. Furthermore, it should be noted that many pharmacological processes other than drug binding to ion channels are analysable in terms of serial Markov chains, and thus, these pharmacological processes are constrained by the mathematical expression derived from the present study.

Author contributions

HZ, JN, and PÅ: designed the research; HZ and PÅ: performed the mathematical work; HZ, JN, and PÅ: analyzed the results; and HZ, JN, and PÅ: wrote the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer AS and handling Editor declared their shared affiliation.