Abstract
Voltage gated potassium channels can be composed of either four identical, or different, pore-forming protein subunits. While the voltage gated channels with identical subunits have been extensively studied both physiologically and mathematically, those with multiple subunit types, termed heteromeric channels, have not been. Here we construct, and explore the predictive outputs of, mechanistic models for heteromeric voltage gated potassium channels that possess either N-type or C-type inactivation kinetics. For both types of inactivation, we first build Markov models of four identical pore-forming inactivating subunits. Combining this with previous results regarding non-inactivating heteromeric channels, we are able to define models for heteromeric channels containing both non-inactivating and inactivating subunits of any ratio. We simulate each model through three unique voltage clamp protocols to identify steady state properties. In doing so, we generate predictions about the impact of adding additional inactivating subunits on a total channel's kinetics. We show that while N-type inactivating subunits appear to have a non-linear impact on the level of inactivation the channel experiences, the effect of C-type inactivating subunits is almost linear. Finally, to combat the computational issues of working with a large number of state variables we define model reductions for both types of heteromeric channels. For the N-type heteromers we derive a quasi-steady-state approximation and indicate where the approximation is appropriate. With the C-type heteromers we are able to write an explicit model reduction bringing models of greater than 10 dimensions down to 2.
Introduction
Potassium channels are an evolutionarily conserved set of membrane spanning proteins (Harding et al., ; MacKinnon, ). While the singular role of these channels is to permit the flux of potassium ions across the impermeable cell membrane, this seemingly simple function is responsible for a wide range of complex physiological processes. In excitable cells, such as neurons, muscle cells and pancreatic β-cells, the primary purpose of potassium channels is to change the shape of, or terminate, action potentials. For non-excitable cells, potassium channels are responsible for functions like volume regulation, maintaining cell shape and electrolyte balance, and neurotransmitter release (Barfield et al., ; Ghatta et al., ; Hoshi et al., ; Miller, ).
Although some of these channels are constitutively open, the majority are gated. Gated potassium channels are divided into three main subsets: ligand gated (Kir), calcium activated (KCa), and voltage gated (KV) channels (Coetzee et al., ; Harding et al., ). Contained within each of these classes are a number of unique channels which all have a distinct set of kinetic properties. The KV family has roughly 40 genes, which each encode a single α-protein subunit. A total of 4 of these subunits are required to form a functional tetrameric KV channel (Harding et al., ; Cordeiro et al., ). Experimentally it has been shown that these 4 subunits need not be identical. When they are identical, the channels are referred to as homomeric channels, but in the cases when more than 1 subunit type is present, a heteromeric channel is formed (Cordeiro et al., ; Al-Sabi et al., ). The existence of these heteromeric KV channels is a relatively recent discovery, and has lead researchers to start investigating heteromeric channel properties in relation to the homomeric channels comprised of similar subunits.
Through a variety of experimental techniques, including mathematical modeling, researchers have begun developing a clearer picture of the role and function of heteromeric KV channels. Working with cone snail toxins, Cordeiro et al. () demonstrated that these snails have naturally occurring compounds designed to target heteromeric KV channels with specific combinations of subunits in their prey. Al-Sabi et al. () showed using concatemeric KV1.1/KV1.2 constructs that heteromeric channels have activation open probability curves that lie between those of the related homomeric channels. Furthermore, they found a nonlinear shift in these open probability curves as additional KV1.2 subunits were included in the heteromeric concatemers. In previous modeling work, we replicated the Al-Sabi et al. () observations about heteromeric channels with a mechanistic Markov model. We also showed the model's predictive ability in relation to a number of cDNA expression experiments (McGahan and Keener, ). While these studies have elucidated certain features of heteromeric channels, the emphasis has been on studying KV channels with little to no inactivation.
Voltage gated potassium channel inactivation has been shown to take two forms: N-type and C-type inactivation. N-type inactivation is believed to be a voltage insensitive process by which a single “tethered ball” segment at the N-terminal of a subunit binds to, and occludes, the open channel pore (Holmgren et al., ; Sukomon et al., ; Hoshi et al., ; Miller, ). Meanwhile, C-type inactivation is believed to be an internal conformational change happening either near the mouth of the channel, or in the selectivity filter (Tan et al., ; Cuello et al., ; Hoshi et al., ; Miller, ). These inactivation mechanisms have been studied in homomeric KV channels using techniques such as cryo–electron microscopy, mathematical modeling, and deletion/mutation experiments (MacKinnon et al., ; Hoshi et al., ; Miller, ; Al-Sabi et al., ). One prominent modeling study conducted in 2011 by Bett et al. (), used mutated KV1.4 channels to derive and fit Markov models for N and C type inactivating homomeric channels.
In spite of this work, and the plethora of experimental techniques and computing power available, the understanding of N-type and C-type inactivation in heteromeric channels is incomplete. What is known regarding heteromeric channel inactivation, is that a single subunit from an inactivating homomeric channel is sufficient to confer inactivation kinetics in a heteromer. By combining specific mutations with the effects of a scorpion toxin, MacKinnon et al. () were able to demonstrate the existence of this property in Shaker potassium channels. Related to this work, Hashimoto et al. () created mutated KV1.4 channels with different numbers of tethered inactivation N-Balls to explore the impact of each additional ball on the homomeric channel. Much later, in 2011, Al-Sabi et al. () showed that the N-type inactivation prevention domains of KV1.6 subunits could be used to reduce, or eliminate, the effects of a single KV1.4 subunit's N-type inactivation in a heteromeric concatemer. Yet to be addressed experimentally, or via modeling, are questions surrounding the impact of including additional N-type or C-type inactivating subunits and the resulting heteromeric channel's overall kinetic properties. In particular:
What mathematical model structure for the homomeric channels is necessary to ensure a heteromeric channel with single inactivating subunit will inactivate?
Do more N-type/C-type subunits cause greater level of channel inactivation and slow the rate of recovery?
Is the relationship between number of N-type/C-type subunits and the level of inactivation linear, or is it nonlinear as was observed with the activation curves of KV1.1/KV1.2 (Al-Sabi et al., )?
Are there distinguishable differences between the impact that C-type versus N-type inactivating subunits have on a heteromeric channel?
Building upon previous experimental and modeling results, with the goal of addressing these types of questions, we propose novel, biophysically detailed mathematical models for heteromeric channels with either N-type or C-type inactivating subunits. We begin by briefly outlining the model used for non-inactivating, homomeric and heteromeric voltage gated potassium channels (McGahan and Keener, ). Then, using the (Bett et al., ) models as a guideline, we construct homomeric N-type and C-type inactivating channel models. Finally, for both N-type and C-type inactivation, we describe how to extend the homomeric models to heteromeric channels, show the response of every heteromeric and homomeric model to a series of complex voltage clamp protocols, and outline any relevant model reductions.
Methods
Voltage gated potassium channels: no inactivation
To model a non-inactivating, homomeric, voltage gated potassium channel, we assume there are 4 identical and independent subunits, each of which can be in either the open or closed conformation. Furthermore, all 4 subunits must be in the open state to have a conducting channel, and only one subunit can change state at a time. With these assumption, the probability of having i subunits in the open state can be described with the master equation:
where α and β are the rates of a single subunit transitioning from the closed to the open state and vice versa. This model has classically been used to describe channels of this type, dating as far back as the first mathematical description of potassium channels by Hodgkin and Huxley (). The model detailed in the groundbreaking work by Hodgkin and Huxley has been shown to be an exact one dimensional model reduction of this more general Markov scheme (Keener, , ). Using the same assumptions made for a homomeric channel, in earlier work we detailed how to model a heteromeric channel lacking significant inactivation kinetics (McGahan and Keener, ). In that work, we described both a general Markov model and detailed a 2 dimensional model reduction for each possible subunit ratio. These reduced systems of differential equations were shown to be globally attracting stable invariant manifolds, meaning the reductions were exact and replicate the behavior of the full system as time is taken to infinity (Keener, , ).
The work below addresses homomeric and heteromeric channels possessing inactivation kinetics. For both N-type and C-type inactivation, we present and justify a homomeric channel model, describe how to extend these models to heteromeric channels, show the response of every model to a series of different voltage clamp protocols, and outline any relevant model reductions.
Voltage clamp simulations
With each homomeric and heteromeric model, for fitting and analysis, we simulate the model using the three different voltage clamp protocols whose schematics are shown in Figure 1.
Figure 1
In order, the voltage clamp protocols are referred to as the activation protocol (Figure 1A), the inactivation protocol (Figure 1B), and the recovery protocol (Figure 1C). The activation protocol fixes the membrane potential at −90 mV and then, by multiples of 10 mV, steps the membrane potential to some new voltage value. From this protocol, 2 pieces of information are extracted: the maximum open probability attained, and the time constant, τ. By tracking the probability of the channel being in the open state, O, during this simulation, the maximum open probability is the maximum value O attains during the activation protocol in response to the value the voltage was stepped to. To remain consistent with previous experimental work, and the relevant experimental data used for model fitting, the second piece of information extracted from the activation protocol, for only the homomeric channels, is a time constant τ. Both the methodology for calculating τ, and the justification behind only calculating τ for homomeric channel model fitting, can be found in McGahan and Keener (). As more experimental data is made available, the fitting techniques and analyses could easily be adjusted to incorporate the entire time course traces seen during these protocols.
The inactivation protocol (Figure 1B) provides steady state inactivation information about the channel. In this protocol, the channel is exposed to 2 different pulses P1 and P2. The P1 pulse is a replica of the activation protocol, but is made to last for 5 seconds. Then, regardless of the voltage the channel is set to, during P2 the voltage is raised to 50 mV. The recorded output of this protocol divides the maximum value of O during P2, by the maximum value found during P1. This resulting ratio is then plotted in response to the value the voltage was stepped to during P1.
The last protocol is the recovery protocol (Figure 1C). This experiment begins with the first pulse P1 by stepping the fixed potential from −90 mV up to 50 mV and holding the voltage there for 5 seconds. After these 5 seconds, the voltage is returned to −90 mV for a variable time duration Δt, until it is forced back to 50 mV during P2. Similar to the inactivation protocol, the value of interest is the maximum probability of the open state O seen during P2 divided by the value seen during P1. This ratio is then plotted in response to the varying time duration Δt to get a metric of how long it takes the channel to return to the closed, non-inactive, state.
Results
Homomeric N-type inactivating channels
Here we detail a Markov model for N-type inactivating channels. Work by Bett et al. () mutated a KV1.4 channel to reduce C-type inactivation capabilities resulting in a channel with primarily N-type inactivation. Using this mutated channel, a 6 state Markov model was proposed and fit for a homomeric N-type inactivating KV1.4 channel by Bett et al. (). This Markov model, with a minor adjustment for clarity as we extend to heteromeric channels, is given below:
Similar to the homomeric channel lacking inactivation described above, there are assumed to be four independently gated subunits all of which must be in the open conformation to have a conducting channel. From the open state O, the subunits can close again, or a single lipophilic N-terminal “ball” region from one of the subunits can bind to the channel opening, thereby transitioning the channel to an inactive state (Holmgren et al., ; Sukomon et al., ). As there are 4 subunits, each containing an N-terminus, we write this transition rate from O to I as 4aI. However, as the unbinding rate of the single bound N-terminal ball is independent of the number of subunits, the transition rate from I to O is written as bI. The differential equations describing this model are:
Heteromeric N-type inactivating channels
Following in the spirit of the proposed model structure by McGahan and Keener (), we assume that a heteromeric channel still requires 4 activating gates, one per subunit type. To account for the N-terminal “ball” binding to the channel pore in a heteromeric channel, there is still only a singular inactive state as in Figure 2. However, the transition rate from O to I now becomes naI, where n is the number of N-type inactivating subunits in the heteromer with 4aI the transition rate in the homomeric channel. A general example Markov model scheme for a channel with 3 N-type inactivating subunits, and 1 non-inactivating subunit is given in Figure 3.
Figure 2
Figure 3
To examine the predictive outputs of this model scheme, we consider heteromeric complexes forming with KV1.1 and Kv1.4 subunits (Harding et al., ; Bett et al., ). In an attempt to accurately reflect the proper channel (and subunit) kinetic properties, forward and backward activating and inactivating rates are taken from Markov models that have been experimentally fit using voltage clamp data. The transition rates for the KV1.1 subunits are drawn from Masoli et al. (). For KV1.4 channel kinetics, we used the model fit by Bett et al. () to voltage clamp data for Kv1.4[K532Y] mutated channels which are known to possess N-type, but lack C-type, inactivation. Parameterizing the subunits in this manner, we took each possible subunit combination and exposed the corresponding heteromeric and homomeric channels to the voltage clamp protocols described in Figure 1. The results of these simulated experiments are shown in Figure 4. We emphasize here that although the simulations here are centered around KV1.4 and KV1.1 heteromeric channels, the methodology and proposed model scheme is widely applicable provided the channels follow the identical and independently operating subunit assumptions.
Figure 4
As response of the KV1.4[K532Y] and KV1.1 homomeric channels to the activating protocol are similar, this experiment does little to distinguish the heteromeric channels. However, in both the inactivation protocol and the recovery from inactivation protocol a trend becomes clear. The first observation to be made is that any number of KV1.4 subunits confer inactivating kinetics upon the channel. Furthermore, as the number of Kv1.4[K532Y] subunits increases, the level of inactivation increases and the rate of recovery decreases. Additionally, these trends are nonlinear. That is, at higher voltage values, there is a noticeably larger change in probability curves moving from the 3:1 heteromer to the 2:2 heteromer, than transitioning between the probability curves of the 1:3 and 0:4 channels. This phenomenon is replicated in the recovery from inactivation curves as well. By increasing the number of Kv1.4[K532Y] subunits, the channel takes longer to recover from inactivation, with the first Kv1.4[K532Y] subunits having a larger impact than the last.
N-type inactivating model: QSS approximation
The largest of these heteromeric models is the 2:2 heteromeric channel consisting of 10 ordinary differential equations. For models of this scale it is useful to find ways to reduce dimension, particularly if the models are to be incorporated into larger neuron models or used to simulate more complex protocols for extended time intervals. It has been shown that for homomeric and heteromeric channels without inactivation kinetics there exists a one (or two in the case of heteromeric channels) dimensional, globally attracting stable invariant manifold (Keener, , ). In other words, the channel kinetics of the full system can be described by one (or two) differential equations. However, the introduction of the N-type inactivating state I, to both the heteromeric and homomeric channel models, complicates the process of finding an invariant manifold to simplify model complexity.
That being said, a second common technique for dimension reduction, is a quasi-steady-state (QSS) approximation (Keener and Sneyd, ). This technique of model simplification relies on assuming, or knowing, that some of the kinetic transitions occur much faster than others. Utilizing the kinetic rates for the homomeric KV1.1 and Kv1.4[K532Y] models that were used to generate the results of Figure 4, we show below in Figure 5, that for these Markov models with N-type inactivation, this may be a reasonable assumption.
Figure 5
In Figure 5, it can be seen that for a significant portion of the physiologically relevant membrane potentials, the rates a1, a2, b1, b2 >> aI, bI, i.e., transitioning between closed states is much quicker than transitioning to the inactive state. Using the KV1.4 homomeric model as an example, we outline the process for using this relationship between the kinetic rates to construct a one dimensional ODE approximation of both the full heteromeric and homomeric KV1.1:Kv1.4[K532Y] Markov models.
Working with the model depicted in Figure 2, the first step in the QSS approximation is to define a new variable X = C0 + C1 + C2 + C3 + O. By summing together the ODEs for C0, C1, C2, C3, O we arrive at the ODEs:
To get these ODEs to be functions of only X and I we must solve for O in terms of X. This is done by making use of the fast equilibrium assumption that a1, a2, b1, b2 >> aI, bI to solve the ODEs of Equation 2 at steady state. Starting with , using the fast equilibrium assumption, the bII and −aIO terms go to 0, and we are left with . The remainder of the ODEs for C0, C1, C2, C3 can be solved in iteration at steady state to get an expression in terms of O. This gives X = κO, where, for the homomeric channel, . Noting that X and I are probabilities of being in a particular state, we have X + I = 1. Therefore, we can fully describe the system with the single ODE:
Here N is the number of N type subunits, which for the homomeric channel case gives N = 4. Using the relationships X = κO and X = C0 + C1 + C2 + C3 + O, it is possible to simulate the one-dimensional system, and then back out the value of any of the desired states C0, C1, C2, C3, O.
For any of the KV1.1:KV1.4[K532Y] heteromeric models, this technique of model dimension reduction with a fast timescale is applied in an identical manner. In all cases, a new variable X is defined as the sum of every closed and open state, and then enforcing the fast equilibrium assumption results in an ODE of the form seen in Equation 4. The only differences between the QSS reductions for each of the different channels is the value of N and κ. We note here that this technique applies generally to models resembling any of these Markov schemes provided the assumptions about the kinetic rates (a1, a2, b1, b2 >> aI, bI) hold. Taking the QSS and full model versions of the KV1.1:Kv1.4[K532Y] heteromeric and homomeric channels, we can compare the approximate and exact model solutions. This comparison is highlighted in Figure 6.
Figure 6
In Figure 6 we see that for the inactivation and recovery protocols, the QSS model performs almost identically to the full model. The only noticeable difference between the models is in the response to the activation protocol. At voltage values greater than −30 mV, the QSS homomeric KV1.1:Kv1.4[K532Y] models consistently predict a greater probability of opening than the full KV1.1:Kv1.4[K532Y] models. The difference in models is best illustrated by looking at the full time simulated traces for the homomeric Kv1.4[K532Y] channels in response to the activation protocol shown in Figure 7.
Figure 7
From Figure 7 it can be seen that the QSS model is immediately able to transition to its maximum probability of opening before inactivation takes place, which is forced by the QSS assumption. Whereas for the full model, the time it takes to reach its maximum probability of opening is long enough to allow for an increase in the probability of channel inactivation. The performance of this approximate solution on a more detailed voltage clamp protocol is provide in the Appendix. The detailed protocol is based on work by Fink and Noble () and is also depicted in the Appendix.
Homomeric C-type inactivating channels
The second form of inactivation to be modeled here, inherent to the α-protein subunits of KV channels, is C-type inactivation. C-type inactivation, while less understood than N-type inactivation, is believed to result from conformational changes happening near the mouth of the pore or in the selectivity filter (Bett et al., ; Tan et al., ; Cuello et al., ). In the same 2011 study modeling N-type inactivation, Bett et al. () also mutated K1.4 channels to delete the N-terminus ball (termed a Kv1.4ΔN channel) allowing them study and model C-type inactivation. This work led to a detailed Markov model that was able to replicate the responses of the Kv1.4ΔN channel to the protocols of Figure 1. Although this model was able to fit the experimental data and reflects the four transitions of each α-protein subunit opening, the biophysical mechanism of C-type channel inactivation is unclear based on the model structure.
We propose a new model structure where C-type inactivation is directly modeled as a conformational change of any of the four subunits from their open state. In particular, each of the four subunits can transition from a closed state C to an open state O, and only once it has reached an open state can it transition to an inactive state I (Figure 8B). Combining 4 identical, independent subunits that adhere to these kinetic transitions gives the full channel Markov model seen in Figure 8A, where the indices denote the number of subunits in the open and inactive states respectively.
Figure 8
C-type inactivating model fitting
To justify working with this new model scheme, we must show that it is comparable in its ability to fit data to the (Bett et al., ) Kv1.4ΔN model. To do this, we first generated simulated data using the published Bett model and parameters in response to the three experimental protocols. Then for the scheme presented in Figure 8A, we create a system of ODEs according to the master equation given in Equation 5.
As C-type inactivation is relatively voltage insensitive (Bett et al., ), the inactivation rates aI, bI were set to be constant parameters, while a1, b1 were given the forms below in Equation 6
for parameters m1, n1, m2, n2. To fit this model to the simulated (Bett et al., ) model data, we employed a bounded global parameter search on the simulated data. By minimizing the sum of squared errors between our model output and the Bett model output, we arrived at a parameter set that generated the model fit seen in Figure 9. For interested readers, more explicit parameter fitting details regarding the optimisation method, error function construction and parameter bounds are provided in the Appendix.
Figure 9
Figure 9 shows that with this parameter set our model can replicate the (Bett et al., ) model outputs. As there were a number of parameter sets found during the global parameter search that had a similar sum of squared errors, it is reasonable to hypothesize that a different parameter set for our model could fit the full experimental data as well.
Heteromeric C-type inactivating channels
Assuming that C-type inactivating subunits are modeled with the scheme presented in Figure 10A and non-inactivating subunits with the scheme shown in Figure 10B, we can construct any heteromeric combination of subunits. As with heteromeric non-inactivating channels and N-type inactivating channels; subunits of a specific type are assumed to be identical, all of a channel's subunits act independently, and all subunit kinetics are fit to the homomeric channels' voltage clamp experimental data. Again, provided these model assumptions hold, the framework is broadly applicable to other C-type inactivating subunits that are characterized mathematically in this manner. Presented in Figure 10C is the sample Markov diagram for a 3:1 heteromeric combination of 3 non-inactivating subunits and 1 C-type inactivating subunit. The Markov diagrams for the 2:2 and 1:3 heteromeric channels are provided in the Appendix in Figures 13, 14 for interested readers.
Figure 10
The corresponding master equations for any heteromeric combination of N non-inactivating subunits and M C-type inactivating subunits is given in Equation 7
As with the N-type models, we explore the impact of adding in additional inactivating subunits by simulating each heteromeric channel model through the three voltage clamp protocols depicted in Figure 1. For the homomeric non-inactivating kinetics we again work with a KV1.1 channel parameterized with transition rates from Masoli et al. (). For C-type inactivation, we use the parameter set generated by the fitting methodology outlined in “Section: C-Type Inactivating Model Fitting” for the KV1.4ΔN channels. The simulation results for the inactivation and recovery voltage clamp protocols across all heteromeric KV1.1:KV1.4ΔN channels are summarized in Figure 11.
Figure 11
In Figure 11, there are two immediate observations to make. The first observation is that like the non-inactivating:N-type heteromers, a single KV1.4ΔN subunit is sufficient to supply some, albeit small, amount of channel inactivation (Figures 11B, C). The second noticeable feature is that as the number of KV1.4ΔN subunits is increased, the level of inactivation is increased and the rate of recovery is decreased. However, unlike with the non-inactivating:N-type heteromeric channels, this relationship between KV1.4ΔN subunits and inactivation appears to be almost linear. Each additional KV1.4ΔN subunit roughly confers an equal amount of inactivation and has a similar impact on the time to recovery.
C-type inactivating model: invariant manifold
As was noted in the previous section, the dimension of the non-inactivating:C-type heteromeric channel models becomes quite large (smallest model is 12 ODEs). For heteromeric channel models described by the master equation seen in Equation 7, we show there is a 3 dimensional invariant manifold which stems from the assumption of identical and independent subunits. Using some results from Keener (, ) we detail the proof for these invariant manifold reductions using the 2:2 heteromer as an example.
Based on Keener (), the binomial distribution gives a vector Pk describing the probability of having k open, non-inactivating subunits. This vector
is an invariant manifold for the set of equations which describe the interactions between 2 non-inactivating subunits. This holds provided the parameter q, the probability of a single subunit being in the open state, satisfies the ODE:
We can find an invariant manifold Pi,j, for 2 interacting C-type subunits, which describes the probability of having i C-type subunits having reached the open state and j C-type subunits in the inactive state. This invariant manifold is given by the multinomial distribution:
so long as the parameters n (probability of a single subunit being open), and h (probability of an open subunit being inactive) satisfy
Note that for either type of subunit, the invariant manifold dimension is 1 less than the number of subunit states as a result of working with probabilities. For example, with the non-inactivating subunits, if q describes the probability of a single subunit being open, then (1 − q) gives the probability of the subunit being closed. Similarly, for C-type inactivating subunits, if we track the probability of a single subunit being open (n) or inactive (h), then the probability that subunit is closed is (1 − n − h).
If Equations 8, 10 give invariant manifolds for 2 non-inactivating subunits and 2 C-type inactivating subunits respectively, we should have that is an invariant manifold for Equation 7 with N = M = 2. To prove this, we begin by rewriting Equation 7 as
where P is an 18 dimensional state vector, and R is the corresponding 18 × 18 rate transition matrix. We show that Pi,j,k, with the parameter dynamics of Equations 9, 11, is an invariant manifold for Equation 7 by demonstrating that
The left hand side of this equality should be thought of as trajectories in the entire space being evaluated on the manifold, while the right hand side gives trajectories of the reduced space. Showing equality implies that once the full 18 dimensional system reaches Pi,j,k, it is restricted to movement along Pi,j,k, giving us an invariant manifold. Proof of this equality is a matter of matrix multiplication which is left to the reader. However as a guide, here we show a single row's worth of calculations using the 18th row of R.
Explicit expressions for Pi,j,k and the rate transition matrix R are provided in the Appendix.
It is reemphasized here that reducing the full Markov scheme in this manner provides an exact, and not approximate, simplified model. For this reason, no comparison simulations between the full and reduced models, as was done in Figure 6 for the N-type QSS model, are required. All future simulations and work can be conducted with the new reduced model and will provide identical results in tracking the open probability. However, the full Markov scheme still provides useful understanding by framing the model in terms of its physiological context as it relates to the opening, closing and inactivating of the α-subunits.
Discussion
In our previous work (McGahan and Keener, ), we showed how to translate from the standard Markov model for a homomeric channel with 4 identical, non-inactivating subunits to a heteromeric channel with 2 different types of non-inactivating subunits. Here we built off of this work by first modeling channel inactivation in a homomeric channel, and then by showing how subunits from an inactivating channel could be incorporated into a heteromeric channel model. In this study we constructed homomeric and heteromeric models for channels possessing either of the two primary inactivation mechanisms: N-type and C-type inactivation. In both cases, the homomeric models were based on experimental data, and the heteromeric models were then a result of the model assumptions with no further fitting. In spite of the limited heteromeric experimental work within the literature to date, we show (and discuss) which known properties our models reproduce and then highlight each of the testable predictions we were able to produce from our models.
N-Type inactivation
The model we used for a homomeric N-type inactivating channel, shown in Figure 2, was identical to the model presented by Bett et al. (). This model structure nicely reflects the known experimental observations that there are 4 identical subunits which must change to the open conformation to enter the conducting state, that only a single N-terminal “ball” can bind to the open channel state, and that this binding is a voltage insensitive transition.
Working with this homomeric N-type inactivating channel model, we described the manner by which subunits of this type would form heteromeric complexes with subunits of a non-inactivating channel. The heteromeric channels were similarly required to have 4 activation steps, with number and type identical to subunit number and type, and with a number of N-balls equal to the number of N-type subunits (Figure 3). To examine the predictions made by our proposed heteromeric model structure, we found fully parameterized homomeric models, for both N-type inactivating and non-inactivating channels that are known to form heteromers (Harding et al., ; Masoli et al., ; Bett et al., ). By simulating each homomeric and heteromeric channel through the three different voltage clamp protocols shown in Figure 1, we observed two key properties of this model.
The first result of this model, which is in agreement with known experimental evidence, is that a single N-type subunit is sufficient to confer some channel inactivation (Figures 4B, C) (Coetzee et al., ; MacKinnon et al., ). The second predictive result of this model, is that while increasing the number of N-type subunits does increase the level of inactivation, each additional N-type subunit has a diminished impact (Figure 4). This non-linear effect is similar to what has been observed and modeled with the changes in open probability curves of non-inactivating KV1.1:KV1.2 heteromers in response to increasing the number of KV1.2 subunits (McGahan and Keener, ; Al-Sabi et al., ). In the case of the non-inactivating KV1.1:KV1.2 channels, the non-linearity is a result of needing all four subunits to be in the open state. This means the channel opening is rate limited by the last subunit, in this case KV1.2, to open. With the N-type:non-inactivating heteromers, the non-linearity observed in the inactivation and recovery protocols can also be explained as a rate-limited effect. With this model, as the number of N-type inactivating subunits is increased, so is the number of N-balls. However, given that only one ball can bind to the pore opening at a given time, there is a competitive binding site saturation for the N-type balls, where adding additional N-balls eventually stops increasing the overall rate at which the first one of them can bind to the pore.
Despite the current experimental complications with studying inactivating heteromeric channels, and thus limited literature on the topic, there is one particularly relevant study to our N-type model done by Hashimoto et al. () looking at modified KV1.4 channels. In their work they were able to construct versions of a KV1.4 channel possessing either 1, 2, or 4 tethered N-balls. The key result, as it relates to our model, was an identical appearance of a non-linear effect of increasing the numbers of N-balls. More precisely, they show that the inactivation rate of a channel with 2 balls is more than half that of a 4 balled channel, and the inactivation rate of a 1 balled channel is more than one fourth that of the full homomeric KV1.4 channel. Additionally, their experiments reveal that as the number of inactivation balls is decreased the level of channel inactivation became less complete. In tandem, these results create a nice qualitative match with our model output given in Figure 4. Finally, Hashimoto et al. () also provides a biological mechanism behind the non-linear impact of each N-ball and the incomplete inactivation by hypothesizing that an electrostatic repulsive force felt between the balls is contributing to the decreased effectiveness of N-ball binding.
C-Type inactivation
In their 2011 work, Bett et al. () published a homomeric C-type inactivating model that was fit to experimental voltage clamp data. While their model was able to fit the data well, the chosen model structure does not help capture or explain the biophysical mechanisms of C-type inactivation that are at play. We proposed a new homomeric C-type model that encodes the C-type conformational change into the model structure, by assuming each of the 4 subunits has a closed, open, and inactive state (Figure 8). To justify using this new model, we first generated data by simulating the published (Bett et al., ) model for KV1.4ΔN channels through the three voltage clamp protocols from Figure 1. Then by performing a bounded, global parameter search we showed that there is at least one parameter set which qualitatively has the same properties as the (Bett et al., ) model (Figure 9). We note that while this may not be the optimal fitting parameter regime, with the limited available data we were still able to qualitatively match the channel behavior.
Utilizing our C-type inactivating homomeric model, we were able to construct heteromeric channels with C-type inactivating KV1.4ΔN channels and non-inactivating KV1.1 channels. The only assumption we needed to make was that each heteromeric channel has 4 total independent subunits, which individually are modeled with the ODEs corresponding to Figure 10A (C-type subunits), or Figure 9B (non-inactivating subunits). Taking each possible combination of these subunit types, we simulated the channels through the voltage clamp protocols and highlighted the response of these channels to the inactivation and recovery protocols (Figure 11). Like the N-type:non-inactivating model, these simulations revealed that a single KV1.4ΔN subunit was enough to make the heteromeric channel inactivate. Furthermore, as the number of C-type subunits is increased, the level of inactivation and time to recovery are increased.
However, in contrast to the N-type:non-inactivating model, this change in inactivation and rate of recovery scales linearly with the number of C-type inactivating subunits. Unlike the inactivation of heteromers with N-type inactivating subunits, there is no binding saturation issue since the inactivation is a conformational change inherent to the individual subunits, which are operating independently of one another. This independence of inactivation is what gives the linear relationship between number of C-type inactivating subunits and level of inactivation.
Implications of inactivation differences
A unique capability of this modeling study is that it allows us to compare between heteromeric channels with N versus C type inactivating subunits. One way to illuminate the similarities and differences is by framing them in the context of ion channel viability. In particular, which combinations of subunits would one expect to see naturally expressed, and which of these combinations could actually be functional channels? One possible hypothesis we propose based on these results is that heteromeric channel viability could be a result of how distinct their kinetic properties, and thus influence on spike timing and dynamics, are from their homomeric counterparts.
The first point reemphasized here is that all heteromeric channels with an inactivating subunit will inactivate thereby distinguishing them from a non-inactivating channel. That being said, there is one glaring difference between subunit types. Due to the linear effects of adding additional C-type inactivating subunits discussed in the previous section, there would be a wider degree of variability among the heteromeric properties than for N-type inactivating heteromers. If all heteromers with inactivating subunits are viable, then those with C-type subunits would cover a wider kinetic range thus markedly changing the window of cellular excitability from a cell containing only homomeric channels. Meanwhile, the N-type inactivating heteromers look remarkably similar to each other in response to the voltage clamp protocols regardless of the number of N-type subunits. Therefore, assuming the N-type inactivating homomer is functional, one might expect these heteromeric channels to be equally viable since interchanging N-type subunits with non-inactivating subunits likely produces indistinguishable neuronal dynamics. Future concatameric and co-expression experiments like the studies by Al-Sabi et al. (, ) and Cordeiro et al. () will be critical for answering these questions surrounding the functionality of heteromeric channels of different subunit ratios.
Model dimension reduction
The final outcome of this work stemmed from examining the model dimension of our various heteromeric channels. Although the smallest model, the N-type homomeric channel, consists of only 6 ordinary differential equations, the largest C-type:non-inactivating heteromeric channel model has 20 differential equations. Individually, any one of these channel models is small enough to include in a full neuron model for simulation. That being said, a common experimental procedure where mathematical modeling has been found useful, is one where two different subunit DNAs are injected into a cell and the resulting distribution of the heteromeric and homomeric channels that form is unknown (McGahan and Keener, ; D'Adamo et al., ; Miceli et al., , ; Hasan et al., ). To compare to experiments of this type, and gain an understanding of the types of heteromeric channels that are forming as a result of this procedure, requires modeling every possible subunit combination. Doing so would necessitate having 39 ODEs for every N-type:non-inactivating channel and 70 ODEs to include every C-type:non-inactivating channel.
Recall that the form and size of these models is enforced by the physiologically grounded model assumptions and not due to phenomenological construction. Therefore, our best option for making these models more computationally tractable, while continuing to adhere to the biological assumptions that require these model structures, is to consider possible model dimension reduction techniques.
For the N-type:non-inactivating heteromeric channels, since the presence of the single inactive state eliminates the possibility of finding a nice globally stable invariant manifold, we found Quasi-Steady State approximations of the full models. By assuming the transitions between closed states are faster than the transition between the open and inactive state, which we validated by looking at experimentally fit rate constants (Figure 5), we showed that each heteromeric N-type:non-inactivating model could be reduced to 2 ODEs (Equation 4). Comparing the results of the full and reduced models showed that the Quasi-Steady State approximation almost exactly replicates the full models performance for the inactivation and recovery protocols (Figure 6). The only noticeable distinction between our full and reduced models is their predictions of open probability curves in response to the activation protocol. The QSS model consistently overestimated the probability of being open in comparison to the full model (Figure 7). This is a consequence of assuming the activation is happening on a fast enough timescale that channels are fully able to open before any inactivation takes place.
For the C-type:non-inactivating heteromeric channels we showed that each model has a globally stable invariant manifold of reduced dimension. Based on a result of Keener (, ), the homomeric C-type channel model (Equation 5) has a two dimensional invariant manifold given by a multinomial distribution. Additionally, work by Keener (, ), has shown that the non-inactivating channel (Equation 1) has an invariant manifold given by a binomial distribution, where the open probability is exactly described by the equations of the Hodgkin-Huxley potassium channel (Hodgkin and Huxley, ). By using the multiplicative structure inherent to invariant manifolds, we demonstrate that each heteromeric C-type:non-inactivating channel has a 3 dimensional invariant manifold which is the product of the corresponding multinomial and a binomial distributions.
Future model uses
While one limitation of this work was a lack of data and existing experimental literature to compare with, this limitation also highlights our model's strengths and at possible future applications. Our model provides a biologically relevant story behind the entire class of less characterized, but highly prevalent, heteromeric ion channels. As lab techniques for future study of heteromeric channels become available; this study can be used to guide predictions, probe experimental outputs, and frame these outputs in terms of the intuition gained from our biophysical model. One immediate application this model is suited for is the examination of coexpression experiments (D'Adamo et al., ; Miceli et al., , ; Hasan et al., ) to see if there are preferred functional stoichiometries of inactivating and non-inactivating subunits. As detailed in McGahan and Keener (), this model formulation, especially the reduced versions, provides an easy pathway to examine the distribution of subunits in inactivating KV channels to see if there is either heteromeric or homomeric preference, or a tendency toward assembly according to a binomial distribution.
A second direction for future exploration and modeling of heteromeric channels involves the existence of heteromeric channels formed with auxiliary beta subunits and pore-forming alpha subunits. While we emphasize here that it is crucial to understand the types of heteromers formed from only alpha subunits, it is also well known that in some KV families there are auxiliary beta subunits that can alter channel properties. The strength of these alterations are dependent on the number of beta subunits, and can take the form of introducing inactivation to an otherwise non-inactivating channel or shifting the activation curves of the pore-forming subunits (England et al., ; Pongs and Schwarz, ; Heinemann et al., ; Xu et al., ). In the case of introducing an auxiliary subunit which confers inactivation upon the pore-forming subunits (Xu et al., ), there is a direct parallel to our Heteromeric N-type inactivating model. For example, a Markov model with 4 non-inactivating alpha subunits and 4 inactivating beta subunits would mimic (Figure 2) where the transition to the I state would now represent the auxiliary beta subunits instead of the N-balls. On the other hand, auxiliary subunits known for shifting activation curves of the pore-forming channel are likely a more complex phenomenon to model as this type of interaction necessarily implies some form of dependence between subunit state transitions (Pongs and Schwarz, ; England et al., ). Once the assumption of independence is violated, these Markov style models would no longer apply, and would therefore require a new model framework.
Conclusion
Here we have proposed a new model for heteromeric channels which are constructed from C/N-type inactivating subunits and non-inactivating subunits. This hypothesized model structure allowed us to make novel testable predictions about the relationship between the number of subunits that confer inactivating kinetics and the level of channel inactivation. These models demonstrated a noticeable difference between the types of inactivation (N versus C type) and the strength of heteromeric channel inactivation. Moving forward, these models can help frame the approach future researchers take to investigate heteromeric channel kinetics, particularly as it relates to channel inactivation.
Statements
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/keesmcgahan/Heteromeric-Inactivating-KChannels.
Author contributions
KM: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. JK: Conceptualization, Project administration, Writing – review & editing.
Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.
Conflict of interest
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.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2024.1418125/full#supplementary-material
References
1
Al-SabiA.KazaS.Le BerreM.O'HaraL.BodekerM.WangJ.et al. (2011). Position-dependent attenuation by kv1.6 of n-type inactivation of kv1.4-containing channels. Biochem. J. 438, 389–396. 10.1042/BJ20102169
2
Al-SabiA.KazaS. K.DollyJ. O.WangJ. (2013). Pharmacological characteristics of kv1. 1-and kv1. 2-containing channels are influenced by the stoichiometry and positioning of their α subunits. Biochem. J. 454, 101–108. 10.1042/BJ20130297
3
BarfieldJ.YeungC.CooperT. (2005). Characterization of potassium channels involved in volume regulation of human spermatozoa. Mol. Hum. Reprod. 11, 891–897. 10.1093/molehr/gah208
4
BettG. C.Dinga-MadouI.ZhouQ.BondarenkoV. E.RasmussonR. L. (2011). A model of the interaction between n-type and c-type inactivation in kv1.4 channels. Biophys. J. 100, 11–21. 10.1016/j.bpj.2010.11.011
5
CoetzeeW. A.AmarilloY.ChiuJ.ChowA.LauD.McCormackT.et al. (1999). Molecular diversity of k+ channels. Ann. N. Y. Acad. Sci. 868, 233–255. 10.1111/j.1749-6632.1999.tb11293.x
6
CordeiroS.Finol-UrdanetaR. K.KöpferD.MarkushinaA.SongJ.FrenchR. J.et al. (2019). Conotoxin κm-riiij, a tool targeting asymmetric heteromeric kv1 channels. Proc. Nat. Acad. Sci. 116, 1059–1064. 10.1073/pnas.1813161116
7
CuelloL. G.JoginiV.CortesD. M.PerozoE. (2010). Structural mechanism of c-type inactivation in k+ channels. Nature466, 203–208. 10.1038/nature09153
8
D'AdamoM. C.ImbriciP.SponcichettiF.PessiaM. (1999). Mutations in the kcna1 gene associated with episodic ataxia type-1 syndrome impair heteromeric voltage-gated k+ channel function. FASEB J. 13, 1335–1345. 10.1096/fasebj.13.11.1335
9
EnglandS. K.UebeleV. N.ShearH.KodaliJ.BennettP. B.TamkunM. M. (1995). Characterization of a voltage-gated k+ channel beta subunit expressed in human heart. Proc. Nat. Acad. Sci. 92, 6309–6313. 10.1073/pnas.92.14.6309
10
FinkM.NobleD. (2009). Markov models for ion channels: versatility versus identifiability and speed. Philosop. Trans. R. Soc. A367, 2161–2179. 10.1098/rsta.2008.0301
11
GhattaS.NimmagaddaD.XuX.O'RourkeS. T. (2006). Large-conductance, calcium-activated potassium channels: structural and functional implications. Pharmacol. Therap. 110, 103–116. 10.1016/j.pharmthera.2005.10.007
12
HardingS. D.ArmstrongJ. F.FaccendaE.SouthanC.AlexanderS. P.DavenportA. P.et al. (2022). The iuphar/bps guide to pharmacology in 2022: curating pharmacology for covid-19, malaria and antibacterials. Nucleic Acids Res. 50, D1282–D1294. 10.1093/nar/gkab1010
13
HasanS.BoveC.SilvestriG.MantuanoE.ModoniA.VenezianoL.et al. (2017). A channelopathy mutation in the voltage-sensor discloses contributions of a conserved phenylalanine to gating properties of kv1. 1 channels and ataxia. Sci. Rep. 7, 1–13. 10.1038/s41598-017-03041-z
14
HashimotoY.NunokiK.KudoH.IshiiK.TairaN.YanagisawaT. (2000). Changes in the inactivation of rat kv1.4 k+ channels induced by varying the number of inactivation particles. J. Biol. Chem. 275, 9358–9362. 10.1074/jbc.275.13.9358
15
HeinemannS. H.RettigJ.GraackH.-R.PongsO. (1996). Functional characterization of kv channel beta-subunits from rat brain. J. Physiol. 493, 625–633. 10.1113/jphysiol.1996.sp021409
16
HodgkinA. L.HuxleyA. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117:500. 10.1113/jphysiol.1952.sp004764
17
HolmgrenM.JurmanM. E.YellenG. (1996). N-type inactivation and the s4-s5 region of the shaker k+ channel. J. Gen. Physiol. 108, 195–206. 10.1085/jgp.108.3.195
18
HoshiT.ZagottaW. N.AldrichR. W. (1990). Biophysical and molecular mechanisms of shaker potassium channel inactivation. Science250, 533–538. 10.1126/science.2122519
19
KeenerJ.SneydJ. (2009). Mathematical Physiology 1: Cellular Physiology, volume 2. New York, NY, USA: Springer. 10.1007/978-0-387-75847-3
20
KeenerJ. P. (2009). Invariant manifold reductions for markovian ion channel dynamics. J. Math. Biol. 58, 447–457. 10.1007/s00285-008-0199-6
21
KeenerJ. P. (2010). Exact reductions of markovian dynamics for ion channels with a single permissive state. J. Math. Biol. 60, 473–479. 10.1007/s00285-009-0271-x
22
MacKinnonR. (2003). Potassium channels. FEBS Lett. 555, 62–65. 10.1016/S0014-5793(03)01104-9
23
MacKinnonR.AldrichR. W.LeeA. W. (1993). Functional stoichiometry of shaker potassium channel inactivation. Science262, 757–759. 10.1126/science.7694359
24
MasoliS.SolinasS.D'AngeloE. (2015). Action potential processing in a detailed purkinje cell model reveals a critical role for axonal compartmentalization. Front. Cell. Neurosci. 9:47. 10.3389/fncel.2015.00047
25
McGahanK.KeenerJ. (2022). Modeling the kinetics of heteromeric potassium channels. Front. Cell. Neurosci. 16:1036813. 10.3389/fncel.2022.1036813
26
MiceliF.SoldovieriM. V.AmbrosinoP.BarreseV.MiglioreM.CilioM. R.et al. (2013). Genotype-phenotype correlations in neonatal epilepsies caused by mutations in the voltage sensor of kv7.2 potassium channel subunits. Proc. Natl. Acad. Sci. 110, 4386–4391. 10.1073/pnas.1216867110
27
MiceliF.SoldovieriM. V.AmbrosinoP.De MariaM.MiglioreM.MiglioreR.et al. (2015). Early-onset epileptic encephalopathy caused by gain-of-function mutations in the voltage sensor of kv7.2 and kv7.3 potassium channel subunits. J. Neurosci. 35, 3782–3793. 10.1523/JNEUROSCI.4423-14.2015
28
MillerC. (2000). An overview of the potassium channel family. Genome Biol. 1, 1–5. 10.1186/gb-2000-1-4-reviews0004
29
PongsO.SchwarzJ. R. (2010). Ancillary subunits associated with voltage-dependent k+ channels. Physiol. Rev. 90, 755–796. 10.1152/physrev.00020.2009
30
SukomonN.FanC.NimigeanC. M. (2022). Ball-and-chain inactivation in potassium channels. Annu. Rev. Biophys. 52, 91–111. 10.1146/annurev-biophys-100322-072921
31
TanX.-F.BaeC.StixR.Fernández-MariñoA. I.HufferK.ChangT.-H.et al. (2022). Structure of the shaker kv channel and mechanism of slow c-type inactivation. Biophys. J. 121, 294a–295a. 10.1016/j.bpj.2021.11.1274
32
XuJ.YuW.WrightJ. M.RaabR. W.LiM. (1998). Distinct functional stoichiometry of potassium channel β subunits. Proc. Nat. Acad. Sci. 95, 1846–1851. 10.1073/pnas.95.4.1846
Summary
Keywords
mathematical modeling, ion channel kinetics, heteromeric potassium ion channels, KV channels, Kv1 channels
Citation
McGahan K and Keener J (2024) Mathematical models of C-type and N-type inactivating heteromeric voltage gated potassium channels. Front. Cell. Neurosci. 18:1418125. doi: 10.3389/fncel.2024.1418125
Received
16 April 2024
Accepted
03 September 2024
Published
08 October 2024
Volume
18 - 2024
Edited by
Richard Anthony DeFazio, University of Michigan, United States
Reviewed by
Gary Richard Mirams, University of Nottingham, United Kingdom
Chon Lok Lei, University of Macau, China
Updates
Copyright
© 2024 McGahan and Keener.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kees McGahan kmcgahan@bu.edu
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.