Ion Channel Modeling beyond State of the Art: A Comparison with a System Theory-Based Model of the Shaker-Related Voltage-Gated Potassium Channel Kv1.1

The mathematical modeling of ion channel kinetics is an important tool for studying the electrophysiological mechanisms of the nerves, heart, or cancer, from a single cell to an organ. Common approaches use either a Hodgkin–Huxley (HH) or a hidden Markov model (HMM) description, depending on the level of detail of the functionality and structural changes of the underlying channel gating, and taking into account the computational effort for model simulations. Here, we introduce for the first time a novel system theory-based approach for ion channel modeling based on the concept of transfer function characterization, without a priori knowledge of the biological system, using patch clamp measurements. Using the shaker-related voltage-gated potassium channel Kv1.1 (KCNA1) as an example, we compare the established approaches, HH and HMM, with the system theory-based concept in terms of model accuracy, computational effort, the degree of electrophysiological interpretability, and methodological limitations. This highly data-driven modeling concept offers a new opportunity for the phenomenological kinetic modeling of ion channels, exhibiting exceptional accuracy and computational efficiency compared to the conventional methods. The method has a high potential to further improve the quality and computational performance of complex cell and organ model simulations, and could provide a valuable new tool in the field of next-generation in silico electrophysiology.


Introduction
Mathematical models of individual ion channels form the building blocks of electrophysiological in silico approaches, allowing the investigation of biophysical mechanisms and the bioelectric activity of excitable and non-excitable cells [1,2]. A variety of whole-cell models of different levels of complexity and abstraction have been introduced for the simulation of ion current kinetics and action potential alterations in neural and cardiac cells, facilitating the prediction of disease processes and the development of therapeutic interventions, which have become an integral part in neuroscience and cardiac electrophysiology [1,[3][4][5][6][7]. Furthermore, single-channel models predicting emergent ion channel drug effects on both cellular and tissue levels are increasingly under consideration in pharmacological research, in conjunction with experimental investigations, opening up an innovative and 54 individual cell measurements for deactivation, inactivation, and ramp protocols (see Figure 1). CHO cells (Chinese hamster ovarian cells), stably expressing rat Kv1.1 channels without Kvβ1 and Kvβ2 subunit expression, were used for measurement of the Kv1.1 macroscopic currents. Electrophysiological recordings were performed with the automated patch clamp system Nanion NPC-16 Patchliner Quattro (Nanion Technologies, Munich, Germany), equipped with EPC-10 HEKA Quadro amplifiers (HEKA Elektronik, Reutlingen, Germany), PatchControlHT software (Nanion Technologies, Munich, Germany), and temperature control in whole-cell configuration. Basic quality criteria for measurements were met, showing an offset voltage of Voffset < 45 mV, seal resistance of Rseal > 200 MΩ (after whole-cell configuration), series resistance of Rseries < 15.5 MΩ, and membrane capacitance of Cslow < 35 pF. Data were further processed based on calculated activation index (AI), maximum currents, and a subsequent manual exclusion of measurements [15].
Macroscopic currents were recorded to activation protocols consisting of a 100 ms long initial-and re-pulse at −80 mV and pulses starting from −90 mV to 80 mV (in increments of 10 mV) of 500 ms duration (Figure 1a). The deactivation protocol applied consisted of an initial-and re-pulse of −80 mV for 100 ms, a depolarization pulse at 70 mV over 300 ms for activation, followed by 300 ms long deactivation pulses from −80 mV to +30 mV in 10 mV steps (Figure 1b). Inactivation characteristics were measured according to a voltage protocol of an initial-and re-pulse of −80 mV for 100 ms, depolarization pulses from −40 mV to 70 mV (increment 10 mV) of 1500 ms duration, followed by an activation pulse of 30 mV for 100 ms (Figure 1c). The ramp protocols considered comprised four intervals of de-and hyperpolarization ranging from −90 mV to 50 mV with varying pulse duration (400 ms, 200 ms, 100 ms, 50 ms) and 400 ms pulse breaks to allow the channels to recover (Figure 1d). Macroscopic currents were recorded to activation protocols consisting of a 100 ms long initial-and re-pulse at −80 mV and pulses starting from −90 mV to 80 mV (in increments of 10 mV) of 500 ms duration (Figure 1a). The deactivation protocol applied consisted of an initial-and re-pulse of −80 mV for 100 ms, a depolarization pulse at 70 mV over 300 ms for activation, followed by 300 ms long deactivation pulses from −80 mV to +30 mV in 10 mV steps (Figure 1b). Inactivation characteristics were measured according to a voltage protocol of an initial-and re-pulse of −80 mV for 100 ms, depolarization pulses from −40 mV to 70 mV (increment 10 mV) of 1500 ms duration, followed by an activation pulse of 30 mV for 100 ms (Figure 1c). The ramp protocols considered comprised four intervals of de-and hyperpolarization ranging from −90 mV to 50 mV with varying pulse duration (400 ms, 200 ms, 100 ms, 50 ms) and 400 ms pulse breaks to allow the channels to recover (Figure 1d).

Data and Data Pre-Processing Considered for HMM and STB Model Parameterization
Model parameterization was based on pre-processed data, excluding cell measurements with seal resistance R seal < 300 MΩ and cell measurements exhibiting a high noise level or seal instabilities, resulting in a sample size of n = 60 cells for the activation curves, n = 37 cells for the deactivation curves, n = 45 cells for inactivation curves, and n = 54 cells for the ramp curves. The measured voltage steps considered for parametrization of the HMM model were limited from −50 mV to 70 mV for the activation protocol and from −80 mV to −30 mV for the deactivation protocol, representing voltage levels at which deactivation occurs after channel activation.

Available HH Model and HMM of the Ion Channel Kv1.1
Several HH models [15,16] and HMM-based approaches were developed for the Kv1.1 ion channel family, modeling their native gating behavior as well as specific ion channeldrug interactions, such as the effect of fluoxetine or syntaxin on channel activation [17][18][19][20][21][22]. The channels were reported as non-or slowly inactivating at room temperature, but exhibited a fast inactivation when co-expressed with Kvβ1 or Kvβ3 subunits [23][24][25][26]. A comparably strong inactivation was similarly observed near physiological temperature even in the absence of β subunits [15].
As current HMMs only reflect the activation behavior of these channels at room temperature, while a possible inactivation is not or only insufficiently considered in the proposed Markov schemes, the currently available HMM approaches can, as a consequence, scarcely be adopted and applied for the simulation of other datasets, in particular at higher or physiological temperature levels. Hence, in order to subsequently provide a reliable juxtaposition of the different modeling approaches, here, we further developed an HMM for simulating the macroscopic current of Kv1.1, that also takes into account the slow and fast inactivation at physiological temperature. In contrast to traditional modeling concepts in computational electrophysiology such as Hodgkin-Huxley or hidden Markov-based models, system identification, a methodology known from the field of control engineering and system theory, deals with the characterization of linear or non-linear systems based on observed input and output data. This approach involves specification of the model structure, estimation of the unknown model parameters, and validation of the resulting model. As the kinetics of an ion channel can be considered as a non-linear system in which the output information is not proportional to the change in the input information, we pursued a non-linear system identification approach for modeling. After a detailed analysis of the measured Kv1.1 macroscopic currents to the given input voltage protocol, the Hammerstein-Wiener (HW) model, which is a block-structured system model, was selected. The HW model consists of a linear dynamic subsystem G(s) between two static nonlinear elements, as shown in Figure 2 [27,28].

Data and Data Pre-Processing Considered for HMM and STB Model Parameterization
Model parameterization was based on pre-processed data, excluding cell measurements with seal resistance Rseal < 300 MΩ and cell measurements exhibiting a high noise level or seal instabilities, resulting in a sample size of n = 60 cells for the activation curves, n = 37 cells for the deactivation curves, n = 45 cells for inactivation curves, and n = 54 cells for the ramp curves. The measured voltage steps considered for parametrization of the HMM model were limited from −50 mV to 70 mV for the activation protocol and from −80 mV to −30 mV for the deactivation protocol, representing voltage levels at which deactivation occurs after channel activation.

Available HH model and HMM of the Ion Channel Kv1.1
Several HH models [15,16] and HMM-based approaches were developed for the Kv1.1 ion channel family, modeling their native gating behavior as well as specific ion channel-drug interactions, such as the effect of fluoxetine or syntaxin on channel activation [17][18][19][20][21][22]. The channels were reported as non-or slowly inactivating at room temperature, but exhibited a fast inactivation when co-expressed with Kvβ1 or Kvβ3 subunits [23][24][25][26]. A comparably strong inactivation was similarly observed near physiological temperature even in the absence of β subunits [15].
As current HMMs only reflect the activation behavior of these channels at room temperature, while a possible inactivation is not or only insufficiently considered in the proposed Markov schemes, the currently available HMM approaches can, as a consequence, scarcely be adopted and applied for the simulation of other datasets, in particular at higher or physiological temperature levels. Hence, in order to subsequently provide a reliable juxtaposition of the different modeling approaches, here, we further developed an HMM for simulating the macroscopic current of Kv1.1, that also takes into account the slow and fast inactivation at physiological temperature. In contrast to traditional modeling concepts in computational electrophysiology such as Hodgkin-Huxley or hidden Markov-based models, system identification, a methodology known from the field of control engineering and system theory, deals with the characterization of linear or non-linear systems based on observed input and output data. This approach involves specification of the model structure, estimation of the unknown model parameters, and validation of the resulting model. As the kinetics of an ion channel can be considered as a non-linear system in which the output information is not proportional to the change in the input information, we pursued a non-linear system identification approach for modeling. After a detailed analysis of the measured Kv1.1 macroscopic currents to the given input voltage protocol, the Hammerstein-Wiener (HW) model, which is a block-structured system model, was selected. The HW model consists of a linear dynamic subsystem G(s) between two static nonlinear elements, as shown in Figure 2 [27,28].  Considering the common patch clamp recordings with voltage step and ramp protocols as system input functions and the measured macroscopic current as the system output  Figure 3, with the measured activation curves from a voltage step protocol as an example. Considering the common patch clamp recordings with voltage step and ramp protocols as system input functions and the measured macroscopic current as the system output function, the Kv1.1 channel model, according to the Hammerstein-Wiener model structure, is shown in Figure 3, with the measured activation curves from a voltage step protocol as an example. The Kv1.1 model with input v(t) and output i(t) nonlinearities of the HW-based Kv1.1 ion channel model was structured as piecewise linear (v'(t) and i'(t)) with two breakpoints between input and output nonlinearities. Note that input and output nonlinearities can also be defined as a sigmoid network, piecewise linear with more breakpoints, saturation, dead zone, wavelet network, one-dimensional polynomial, or other elements known from control engineering. Here, we considered different input and output nonlinearities and adopted a model with the same type of input and output nonlinearities.

Mathematical Concepts of Ion Channel Modelling
In general, the HW-based Kv1.1 ion channel model can be described as: For Kv1.1, the system input is defined as v(t), i.e., the voltage signal according to the applied protocol, and system output is i(t), i.e., the measured macroscopic current. According to the definition of the block-structured HW model, it is necessary to define intermediate input functions, vi(t), and intermediate output functions ii(t). The intermediate input is the output of the input nonlinear element and the input to the linear element G(s). Analogously, the intermediate output is the output of the linear element G(s) and the input of the output nonlinear element. The intermediate input and the output functions are defined in Equation (2) and Equation (3), respectively.
The linear element G(s) is the transfer function (TF), which represents the differential equation of the dynamic behavior of the system. The TF is a mathematical represen- The Kv1.1 model with input v(t) and output i(t) nonlinearities of the HW-based Kv1.1 ion channel model was structured as piecewise linear (v (t) and i (t)) with two breakpoints between input and output nonlinearities. Note that input and output nonlinearities can also be defined as a sigmoid network, piecewise linear with more breakpoints, saturation, dead zone, wavelet network, one-dimensional polynomial, or other elements known from control engineering. Here, we considered different input and output nonlinearities and adopted a model with the same type of input and output nonlinearities.
In general, the HW-based Kv1.1 ion channel model can be described as: For Kv1.1, the system input is defined as v(t), i.e., the voltage signal according to the applied protocol, and system output is i(t), i.e., the measured macroscopic current. According to the definition of the block-structured HW model, it is necessary to define intermediate input The linear element G(s) is the transfer function (TF), which represents the differential equation of the dynamic behavior of the system. The TF is a mathematical representation between an intermediate input and an intermediate output function of the system. Hence, the TF of a linear system is defined as the ratio of the Laplace transform of the output to the Laplace transform of the input, where all initial conditions are assumed to be equal to zero (see Equation (4)) [29,30].
According to Figure 3, the system output is represented by the measured Kv1. For channel activation, deactivation, and inactivation, the linear part of the Kv1.1. ion channel model G(s) with three poles and two zeros is given by Equation (5), representing the TF in Laplace domain.
Mathematical transformations can now be used to determine the differential equation of the system in the time domain, i.e., a third-order differential equation that describes the kinetic characteristics of the ion channel, representing the opening behavior of the channel at different voltage levels.
The transfer function TF in the time domain, thus, represents the so-called behavioral differential equation (BDE) and can be denoted as: where a k , k ∈ (0, 3) and b k , k ∈ (0, 2) are coefficients of TF and BDE, when all initial conditions are equal to zero. It is important to emphasize that all a k coefficients are positive, which can be explained by the transient response of the system, but also results from system identification. In terms of the former, we can conclude that the obtained ion channel model is "stable" according to the stability criteria in control theory. The model was obtained in MATLAB using the System Identification toolbox, [31] and estimated using PEM (prediction error minimization). Fitting results were determined by using RMSE values. Figure 4 shows the results of the Kv1.1 STB model parameterization.
Cells 2022, 11, x FOR PEER REVIEW 7 of 28 inactivation properties of the ion channel, ramp protocols in turn provide a continuous recording of the overall dynamic behavior over a large voltage range, which is essential for reliable system identification. Based on the available experiments, the voltage data of the ramp protocol was, thus, used as the input function for system identification of the Kv1.1 STB model. Figure 4 shows the corresponding results of the estimation of the STB model. Detailed information on model parameterization can be found in Appendix A.  If we look at the whole system, having the theoretical considerations in mind, the final Kv1.1 model is described by a nonlinear system with regularly coupled subsystems. i (t) + 1.596 · 10 9 i (t) = −6.268 · 10 8 One of the first and most important steps in system identification is the selection of the input based on prior knowledge and experiment design [27]. In comparison to voltage step protocols commonly used to determine the specific activation, deactivation, and inactivation properties of the ion channel, ramp protocols in turn provide a continuous recording of the overall dynamic behavior over a large voltage range, which is essential for reliable system identification. Based on the available experiments, the voltage data of the ramp protocol was, thus, used as the input function for system identification of the Kv1.1 STB model. Figure 4 shows the corresponding results of the estimation of the STB model. Detailed information on model parameterization can be found in Appendix A.

The HMM-Based Kv1.1 Model
In order to represent the possible conformational states and structural changes underlying channel gating adequately and sufficiently, kinetic schemes of HMMs are derived based on the specific protein structure and known functional properties of the ion channel are additionally considered. The specific structure and investigated kinetic characteristics of the Kv1.1 channel that form the basis for model derivation are briefly summarized below.
The conductivity of voltage-gated potassium (Kv) channels depends on protein conformational changes in response to membrane depolarization [32]. The Kv pore-forming protein consists of 4 α-subunits, where each subunit is composed of six transmembrane segments (S1-S6) and intracellular N-and C-terminal domains, responsible for inactivation of the channel. The first four segments comprise the voltage sensor domain (VSD), segment S5, and S6 form the ion-conducting pore (PD) of the channel, as shown in Figure 5. Positively charged amino acids within the S4 segment trigger movements of the sensor in response to changes in membrane potential, which are transmitted to the pore via the S4-S5 linker for controlling the opening and closing of the channel [32][33][34]. Inactivation occurs by both a rapid N-type inactivation caused by the cytoplasmatic N-terminal sequence occluding the channel pore in the open state and by C-type inactivation, which is a slower time-dependent conformational change, leading to a narrowing of the outer mouth of the channel pore [34]. The α-subunits of Kv1.1 channel of mammalian cells lack the N-terminal sequence, but the proteins, however, exhibit a fast inactivation when complexed with subunits or auxiliary proteins that contain this domain and substitute the functionality, such as Kv1.4 or Kvβ1 and Kvβ3 [35][36][37][38]. In vivo, the channels are typically assembled with peripheral β-subunits, which modify the surface expression of these channels in addition to the gating behavior [36,37]. As recently demonstrated, physiological temperature equally provokes a fast inactivation in Kv1.1 channels, even in the absence of Kvβ1 and Kvβ3 subunits, emphasizing the important role of temperature on channel kinetics and function [15]. As found in several optimization runs, a better fit of the data was obtained when assuming a transition possibility between the fast and slow inactivation, suggesting a linkage of the two inactivation modes. Therefore, a direct transition path between IC2 and IN was considered in the final model approach [39]. Forward transition rates α, λ, σ and backward transitions β, η, ε are voltagedependent and described by first-order exponential functions: where αi and βi represent specific gating parameters and V the applied voltage. c, d, m, k, x, and y denote rate constants without voltage-dependence. Defining ( ) as the probability of being in a specific state  As found in several optimization runs, a better fit of the data was obtained when assuming a transition possibility between the fast and slow inactivation, suggesting a linkage of the two inactivation modes. Therefore, a direct transition path between IC2 and IN was considered in the final model approach [39]. As found in several optimization runs, a better fit of the data was obtained when assuming a transition possibility between the fast and slow inactivation, suggesting a linkage of the two inactivation modes. Therefore, a direct transition path between IC2 and IN was considered in the final model approach [39]. Forward transition rates α, λ, σ and backward transitions β, η, ε are voltagedependent and described by first-order exponential functions: where αi and βi represent specific gating parameters and V the applied voltage.  Forward transition rates α, λ, σ and backward transitions β, η, ε are voltage-dependent and described by first-order exponential functions: where α i and β i represent specific gating parameters and V the applied voltage. c, d, m, k, x, and y denote rate constants without voltage-dependence. Defining P S i (t) as the probability of being in a specific state S i at time t leads to the equation for the time evolution of the channels' open probability P O (t) [2,40]: where the first three terms represent transitions entering the open state O and the term furthest to the right transitions leaving the open state O.
Since HMMs model the current through a single ion channel, optimization based on measured whole-cell currents requires estimating the number of ion channels in addition to the model parameters for simulating the macroscopic current. For sufficiently large numbers of the same channel, the fluctuations in the stochastic opening of individual ion channels average out and the quantities in Equation (10) can be replaced by their macroscopic interpretation. Moreover, the probability of being in state S i can be interpreted as the fraction of channels in S i . The transition probabilities become rate constants, r i,j , which describe the number of channels that change from S i to S j in a given time period [2,40]. The macroscopic current I Kv1.1 is given by the open probability P O , the ion channel number N c , the single channel conductance g Kv1. 1 , and the reversal potential E K : The rate constants were parameterized using a particle swarm optimization (PSO) algorithm from the Global Optimization Toolbox (MathWorks Inc., Natick, MA, USA) based on averaged activation and deactivation measurements. The number of sample cells considered for the activation and deactivation currents differed, and since the magnitude of the macroscopic current varied considerably from cell to cell, the magnitude of the macroscopic currents at similar voltage levels also showed considerable deviations between patch clamp experiments assessing the activation and deactivation characteristics. To account for these variations, the number of ion channels N c was individually optimized for each measurement protocol [41]. For the given dataset, the channel number was determined as N c_act = 3088 for the measured activation and N c_deact = 2588 for the deactivation currents. The final model parameters are summarized in Table 1.   Figure 7 shows the corresponding simulation results of model parametrization (RMSE act_HMM = 0.0714, RMSE deact_HMM = 0.1098). For detailed information on model parametrization and simulations, see Appendix A.

Rate Constants and Parameters
The basic idea of HMMs is to model the specific changes in the conformational states of the protein represented by the different states in the model. To determine whether the transitions and the occupancy of states in the HMM in response to a stimulus corresponded to the underlying kinetics of the channel, we simulated the model stochastically by generating a random sequence of states using the hmmgenerate function in MATLAB (MathWorks Inc.). The basic idea of HMMs is to model the specific changes in the conformational states of the protein represented by the different states in the model. To determine whether the transitions and the occupancy of states in the HMM in response to a stimulus corresponded to the underlying kinetics of the channel, we simulated the model stochastically by generating a random sequence of states using the hmmgenerate function in MATLAB (MathWorks Inc.).  To model the specific Kv1.1 channel conductance, Ranjan et al. [15] adapted the original model of Hodgkin and Huxley [42] (see Appendix D) describing the non-linear potassium conductance in the squid giant axon and added an additional gate h to account for channel inactivation, with single gates for activation and inactivation (see Equation (12)). with The process of model adaption and fitting of the Kv1.1 HH model can be briefly summarized as follows: the steady state variables m ∞ and h ∞ were fitted to single Boltzmann functions: where V 1/2 denotes the half activation and inactivation voltage, k the slope factor, and A the starting point. The time constant for activation τ m was fitted by two Boltzmann equations, and a single Boltzmann equation was again used for τ h :   Normalized conductivities of measured current traces from activation voltage step protocols between −40 mV and 50 mV were fitted for each cell and temperature level (15 • C, 25 • C, and 35 • C). Single-cell models that had a residual sum of squares (RSS) less than 0.36 were considered, and the median values for each gating parameter and temperature level were used for the final model [15].
To account for the temperature-dependent conductivity of the Kv1.1 channel, the median gating parameters of h ∞ , τ m and τ h obtained at each temperature were further fitted with Q10 functions. In comparison, m ∞ was considered to be temperature-independent, despite different values of gating parameters in the revised HH model. The model equations and gating values of the proposed model by Ranjan et al. [15] are given in Equations (17)- (23).
with h ∞ Q 10 = (0.032 · • C) − 0.365 (18) Table 2 summarizes the key features of the HH, HMM, and STB models, including the number of unknown parameters to be optimized, the extent of the mathematical description of the models, and the data used for model parameterization. Detailed information on model parameterization is provided in Appendix A. For the model evaluation and verification, different voltage protocols performed to determine the channel kinetics (see Figure 1) were simulated using each of the three model approaches and compared accordingly. Corresponding simulation results of the Hodgkin-Huxley formalism, the developed hidden Markov model, and the system theory-based approach for the activation, deactivation, inactivation, and ramp protocols are shown in Figure 9.  Since, in contrast to the developed approaches, the HH model was parametrized to fit the normalized currents and, thus, defined only for them, the simulated currents of the activation, deactivation, and inactivation protocols were each normalized to the maximum measured current at 70 mV for comparison. For the ramp currents, the maximum value of the entire trace was used for normalization. The goodness of fit of the simulated current curves was evaluated directly using the root mean square error (RMSE) and averaged over all voltage levels for both the normalized (RMSE norm ) and absolute (RMSE abs ) currents (Equation (24)).
The developed hidden Markov model and the system theory-based model outperformed the HH model in terms of data fitting and reproduced the specific activation, deactivation, and especially the recorded ramp currents very accurately. Remarkably, the activation currents simulated with the STB model were almost identical to the measured current, as shown by the obtained RMSE value, summarized in Table 3. The deviations of the HMM in the obtained deactivation curve at −30 mV, showing an increase in the current after the corresponding deactivation, could be explained by the high-voltage level, which naturally led to an activation of the channels (see Figure 9e). In turn, the disturbances in the STB model resulted from the capacitive spikes that were not filtered out and removed from the measured current traces. Because of these spikes, the model did not reproduce the raw output data more accurately (Figure 9f).
Similarly, the fast inactivation could be modelled with high precision by the newly developed STB compared to the HH approach, that revealed a too strong and prolonged inactivation. A slightly higher RMSE, in turn, was obtained for the simulation of the inactivation curves by the HMM due to moderate deviations of the absolute currents. In general, however, the kinetics correlated well with the measured current dynamics, which was also an acceptable modeling result for the developed HMM. Thus, both models, which were parametrized on a few single-current curves only (see Table 2), were suitable for different input functions and were able to simulate the specific Kv1.1 current, which serves as the first step in verifying and proving the validity of the model. Additional simulations using the HMM and STB approaches performed with action potential (AP) and recovery protocols can be found in Appendix B, Figure A1, which also shows useful simulation results and confirms the potential of the new STB-based modeling approach.
For a thorough evaluation of the accuracy of the models, basic electrophysiological parameters describing the activation, deactivation, and inactivation properties were extracted and compared.
Activation characteristics were evaluated using the conductance voltage relation and the time constant for activation, measured by the activation protocols performed. For this purpose, the normalized conductivities calculated from the peak currents at each voltage step (Equation (25)) were plotted as a function of the test pulse voltages and fitted to a Boltzmann function (Equation (26)): where G max is the maximal conductance measured at a voltage step of 70 mV, V 1/2_act the hemi-activation voltage, and k act the slope factor. The activation time constant was determined by fitting a single exponential function to each individual current curve from the start of the stimulus to the peak current: Tail currents evoked by hyperpolarization pulses following a depolarization step of 300 ms duration were measured and analyzed to determine the deactivation properties. Each individual tail current obtained by the deactivation protocol was fitted to a single exponential function to estimate the time constant of deactivation: with τ as the time constant of deactivation and A 1 and A 2 the initial and final values, respectively. Inactivation characteristics were determined based on the steady-state availability protocols performed that included conditioning pulses of longer duration at different voltage steps to establish a steady-state inactivation after channel activation, followed by a depolarizing voltage step (to activate channels still in an activatable state). The inactivation time constants were calculated based on the activation pulse by fitting a single exponential function from peak to steady state for each current trace according to Equation (28). The half-inactivation voltage V 1/2_inact and the slope factor of inactivation k inact were, again, calculated by fitting the normalized peak currents of the depolarizing voltage step to a Boltzmann function according to Equation (29): Slow voltage ramps were used to determine the voltage at which the channels had maximum conductance V max_cond . Following Ranjan et al. [15], the peak value during the rising phase of the first ramp was used as the parameter V max_cond .
The calculated and extracted electrophysiological parameters are shown in Figure 10 and summarized in Table 3.

Discussion
Single-channel modeling is a central component of computational electrophysiology. Today, extensive experimental investigations and a steadily growing body of knowledge about ion channels enable the development of highly detailed models that simulate the specific gating behavior and the bioelectric properties of ion channels. The increasing biophysical detail, however, also inevitably leads to high computational costs, which, to some extent, limit both the construction and the application of complex wholecell models, especially for simulations on the tissue and organ level. Hence, while detailed HMMs that map the protein structure and better address the processes behind channel gating are mainly considered in biomolecular and pharmacological research, HH models, for example, are still the golden standard in neuroscience, since they provide a low computational cost method and, thus, a high integrability into complex models to represent the electrophysiological activities of cells, tissue layers, and up to whole organs.
Beyond conventional methods, following the phenomenological approach of Hodgkin and Huxley, we proposed for the first time a new system theory-based concept of deterministic ion channel modeling and the simulation of ion currents, which provide an easy-to-use method with remarkable performance and accuracy, especially with respect to the structurally comparable HH models. Using the example of Kv1.1 (KCNA1) delayed rectifier channels, which are strongly expressed in the central and peripheral nervous system and "regulate" neuronal subthreshold excitability and spike initiation [20][21][22]24], the newly introduced method was compared with the concepts of the HH model and HHM, and evaluated on several parameters relevant in the computational modeling of cellular electrophysiology. The measured activation characteristics with a half activation voltage of V 1/2_act_measured = −22.45 mV and slope factor k act_measured = 10.81 mV were best reproduced by simulations with the HMM (V 1/2_act_HMM = −22.64, k act_HMM = 11.82 mV). For the STB model, the curve and, thus, the half activation voltage were slightly shifted towards a more depolarized value, but comparable results to the HH model could be obtained with V 1/2_act_STB = −18.39 and k act_STB = 14.97 mV relative to V 1/2_Act_HH = −14.94 and k act_HH = 9.913 mV. With respect to the activation time constant, both the HMM and the STB model better reflected the actual voltage-dependent dynamics of activation by showing a faster activation time at higher clamp voltages and a slower activation as the voltage decreased, compared to the HH model with the same time constant over the entire voltage range. However, the activation in the STB model was instantaneous and, thus, somewhat too fast, while the activation in the HMM, especially at lower voltages, was too slow compared to the measured values. The simulation results for the deactivation of the HMM and STB models revealed a slower deactivation, but they, again, better reflected the measured deactivation behavior compared to the HH model, as shown by the determined deactivation time constants (see Figure 10 and Table 3).
The model simulations of the inactivation curves showed a slightly better, but comparable half inactivation voltage compared to the HH model with respect to the measured parameters for the HMM model (Figure 10b). In contrast, the STB approach again outperformed the accuracy of the HMM and HH models, and showed a nearly perfect fit of the measured inactivation time constants; see Figure 10e and Table 3.
Taking all the results obtained into consideration, both the newly developed HMM and the STB approach provided an accurate modeling of the channel kinetics that better reflected the underlying dynamics of the channel in response to various input functions than the established HH model used here as benchmark or state-of-the-art model. In particular, the HMM and the STB models provided two valuable new approaches for ion channel modeling and the simulation of the Kv1.1 current at a physiological temperature.

Discussion
Single-channel modeling is a central component of computational electrophysiology. Today, extensive experimental investigations and a steadily growing body of knowledge about ion channels enable the development of highly detailed models that simulate the specific gating behavior and the bioelectric properties of ion channels. The increasing biophysical detail, however, also inevitably leads to high computational costs, which, to some extent, limit both the construction and the application of complex whole-cell models, especially for simulations on the tissue and organ level. Hence, while detailed HMMs that map the protein structure and better address the processes behind channel gating are mainly considered in biomolecular and pharmacological research, HH models, for example, are still the golden standard in neuroscience, since they provide a low computational cost method and, thus, a high integrability into complex models to represent the electrophysiological activities of cells, tissue layers, and up to whole organs.
Beyond conventional methods, following the phenomenological approach of Hodgkin and Huxley, we proposed for the first time a new system theory-based concept of deterministic ion channel modeling and the simulation of ion currents, which provide an easy-to-use method with remarkable performance and accuracy, especially with respect to the structurally comparable HH models. Using the example of Kv1.1 (KCNA1) delayed rectifier channels, which are strongly expressed in the central and peripheral nervous system and "regulate" neuronal subthreshold excitability and spike initiation [20][21][22]24], the newly introduced method was compared with the concepts of the HH model and HHM, and evaluated on several parameters relevant in the computational modeling of cellular electrophysiology.

Model Accuracy
The introduced STB model, parametrized on the ramp data only, allows the accurate simulations of the specific kinetics of the Kv1.1 channel and fits almost perfectly with the measured currents for the different voltage protocols performed (see Figure 9 and also Figure A1), even in a currently highly simplified and well-interpretable form where only two breakpoints were used to approximate the nonlinear input and output function. The accuracy could be further improved by considering additional breakpoints. Figure A2 in Appendix C shows an example simulation of the ramp data using 10 breakpoints in the STB model with an almost perfect fit. However, a higher number of breakpoints resulted in a more complex system description, represented by an even higher order and a less interpretable differential equation in the time domain.
As shown by a direct comparison with a recently published HH model of the Kv1.1 and the new HMM model developed here based on the same experimental data, the STB model outperformed the established models in accuracy and better reproduced the specific activation, deactivation, and inactivation properties of Kv1.1 channels at a physiological temperature. It is important to note that the HH model, used as a benchmark for comparison and model validation, also accounted for the temperature-dependent modulation of the channel kinetics and was parameterized based on the activation curves of different temperature levels, i.e., 15 • C, 25 • C, and 35 • C. For this reason, the HH model represented an average best model for simulating the Kv1.1 current within this temperature range, but did not perfectly match the measured currents at a single temperature. However, simulation results that were within the deviations of the HH model were considered sufficiently reasonable and valid.
Comparable results were obtained for the newly developed HMM in terms of fitting the experimental data to the HH model. The optimization of the HMM to the activation data only allowed an almost perfect simulation of the activation curves, while the deactivation and inactivation characteristics were not represented at all. Furthermore, parametrization based on the ramp curves, as performed for the STB approach, did not lead to a satisfying modeling result. The HMM model was, thus, finally parametrized based on the activation and deactivation curves, which also allowed the inactivation to be adequately represented by the model approach and acceptable model simulations of all voltage-protocols (see Figures 9 and A1). However, the model showed a lower accuracy with regard to the inactivation characteristics. Therefore, in a next step, more attention should be paid to the inactivation path, e.g., by considering additional inactivation protocols in the model parameterization, experimental investigations, and an appropriate redefinition of the number of states, representing the slow and fast inactivation, in order to improve the validity of this newly introduced hidden Markov-based Kv1.1 model.

Model Complexity, Explainability, and Adaptability
Compared to the HMM, but similar to Hodgkin and Huxley, the STB approach is entirely data-driven and does not take into account any electrophysiological knowledge, which, currently, does not allow for inference or insights into the inherent channel gating mechanisms by the model approach. By contrast, even at a highly simplified level, the kinetic schemes of HMM, which map the transitions between different conformational states, offer better explainability compared to the HH and STB models, and the study of specific modifications in the opening and closing behavior of channels, as particularly needed, for example, in pharmacological studies. Moreover, since HMMs describe the dynamics of single channels, they provide a high degree of flexibility and allow its application to different datasets with varying dynamics or current amplitudes by adjusting the rate constants or number of ion channels. HH models, as well as the newly introduced STB approach, always represent the measured macroscopic currents and are valid only for a specific dataset. Therefore, a direct adoption to other experimental data, sample populations, or cells with varying ion channel compositions, is usually not possible without an appropriate and comprehensive reparameterization.
However, the proposed HMM represents a simplified kinetic scheme derived solely on the basis of macroscopic currents and does not take into account further electrophysiological studies such as single-channel recordings or structural studies of protein conformation, which limits the degree of the explainability and adjustability of this first HMM of the Kv1.1 channel. Furthermore, with respect to the inactivation characteristics, no characterization of the slow and fast inactivation was performed, e.g., an assessment of the respective proportion using specific blockers. Additionally, the assumption of a possible transition between and, thus, an interaction of the slow and fast inactivation implemented by a cross-link between the IC2 and IN1 states was based only on achieving a better modeling result as shown in several optimization runs, but without experimental validation. Thus, the states in the model do not correspond to the actual protein conformational states and microscopic conformational changes of the protein, but can be viewed as aggregates of molecular configurations grouped into a set of distinct functional states separated by large energy barriers [1].
Despite the aforementioned simplifications, the HMM model allowed an accurate and reliable simulation of the different measured kinetics, as shown by the occupancy diagrams. The occupancy of states was consistent with the measured and known kinetics, which confirmed the validity of the proposed kinetic scheme and parameterization for modeling the kinetics of the Kv1.1 channel.

Computational Burden
Together with the complexity and level of detail, the high computational cost is one of the major drawbacks limiting the application of HMM. Even simplified kinetic schemes, such as the one developed in this work, include a great number of parameters and states that are implemented in the model by a set of first-order differential equations, implying the need for a very a high computational effort not only in terms of simulation runtime, but also for parametrization. In contrast to the HMM and HH methods, the system theory-based approach significantly reduced the typically huge set of differential equations in the HMM approach to one single higher-order differential equation that describes the current-voltage relation of the ion channels as a nonlinear system with regularly coupled subsystems. This enormously reduces the computational cost for parameterization and model simulation.
Together with the remarkable model accuracy, this represents the main advantage of the newly developed STB model compared to the traditional modeling approaches in electrophysiology.
For HMM in particular, the large number of parameters relative to the comparatively few data also increases the risk of overfitting and, thus, limits the predictive power and reliable simulation of additional data. Therefore, it makes sense to keep the HMM as simple as possible by involving different measurement protocols in model optimization. However, if more data were included in the model optimization, the time for parameterization would increase again. For the developed HMM parametrized on the activation and deactivation curves, each optimization run took about 30 h on a high-performance computer with 12 cores working in parallel for model parametrization. By contrast, the parameterization of the STB model function based solely on the ramp curves was, for example, performed in less than 10 min using the same computer infrastructure with MATLAB (System Identification toolbox, MathWorks Inc.).
Compared to the HH models, HMM also had, on average, a higher computation time, even with a smaller number of states, as shown, for example, in a study by Andreozzi et al. [1], which yielded a 5% higher runtime of a simplified HMM compared to the corresponding HH model. However, given the simulation results obtained, which showed excellent accuracy compared to the HH approach, the increased computation time was considered to be acceptable. For our Kv1.1 simulations, the runtime of an example cell with 3500 individual Kv1.1 channels was about 20 times higher for the HMM than for the HH and the STB model, with the latter requiring less than 1 s.

Experimental Data for Model Parameterization
It is important to note that electrophysiological studies are generally time-consuming, and obtaining representative, quality-assured results usually requires a high experimental effort. The experimental data used in this study are publicly available and include measured whole-cell currents from transfected cells, stably expressing Kv1.1 channels recorded with different voltage protocols. For phenomenological modeling, the data required for model parametrization were rather small and comparable for all modeling approaches examined in this work. They included measured macroscopic currents from patch clamp recordings with standard voltage step protocols to characterize the activation, deactivation, and inactivation characteristics. In order to fully characterize the kinetic properties and improve the validity of HMMs, however, extensive experimental investigations are required, such as single-channel patch clamp measurements, determination of fast and slow inactivation and possible cross links, or structural studies to gain a deeper understanding about the protein conformational states. All these together increase the experimental effort required for HMM development and validation in general enormously compared to the HH model or, in particular, to the newly proposed system theory-based modeling approach.

Which Method Should Now Be Chosen? When, How, and Why?
The three different modeling approaches presented in this work all have both strong advantages and disadvantages, and should always be selected with respect to the particular application. Table 4 summarizes the three modeling approaches by qualitatively comparing the key parameters in computational electrophysiology. The system theory has been an established tool for modeling physical or biological processes for decades, and it is used traditionally in the field of control engineering. In this work, we introduced the concept of a transfer function for the kinetic characterization of single ion channels for the first time. We investigated the extent to which its properties could be used to simulate the activation, deactivation, and inactivation of channels without knowing the intrinsic biological and physical mechanisms, but only using the data characteristics of the input and output function of the "system", which is presented by only the one third-order differential equation, taking the input and output nonlinearities into account. Today, available software tools, such as MATLAB, allow for an easy and automated characterization of the transfer function of the biological system, enabling simple and fast model parameterization compared to the conventional methods such as the HH model and HMM. With this easy to use parameterization strategy, this strongly data-driven modeling approach can be adapted simply to different datasets of sample populations with varying ion channel composition, and could make the system theory-based modeling approach the method of choice for high-performance simulations at the tissue and organ level. Further investigations could show whether and to what extent this concept can also be applied to other ion channel types with divergent kinetics, such as channels with a slow inactivation (e.g., Kv3.1) or constant activation (e.g., Kv7.1), after an appropriate system identification.
In contrast, by embedding knowledge from biophysical and structural studies, the HMM allows a detailed modeling of the specific functionality and structural changes underlying channel gating, representing possible dependencies of activation and inactivation, transitions from closed to inactivated states, or multistep activation processes. In particular, ligand-or second messenger-dependent changes as well as drug-induced effects on specific conformational states and, thus, on the functionality and kinetics of the channel can be investigated at the microscopic and macroscopic level by appropriate kinetic schemes, as ultimately required in pharmacological or molecular-biological investigations. For these applications, the Markov models, which take into account the inherent gating properties and better address the stochastic gating behavior, represent a perfectly suitable method despite the higher experimental effort and computational load [1,11]. Moreover, the HMM with a sufficient complexity and low computational cost used in whole-cell applications can overcome the limitations of the currently most widely used HH approaches, for example, by better accounting for the complex interplay of ion channels, calcium dynamics, or specific responses to changing environmental conditions such as the temperature, pH, or ionic composition. To this end, HMMs are increasingly considered for detailed modeling approaches to further improve the reliability and validity of complex single-cell applications. However, we can expect that if extensive experimental data-representing mechanisms such as druginduced effects, changes to environmental conditions, or intracellular ionic compositions are available, appropriate STB models, because of their simple parameterization, could also be introduced.
In summary, the system theory-based modeling approach combines the positive features and properties of both the HH model and HMM. The proposed concept outperformed the HH model and HMM in accuracy, although it strongly abstracted the underlying elec-trophysiological mechanisms, while overcoming the current computational limitations of the HMM. In particular, for applications requiring high computational power, this newly introduced modeling approach offers a promising new possibility that could be used alongside or even instead of HH-based ion channel models in computational electrophysiology, while further improving the simulation accuracy and runtime. Thus, beyond single-cell applications, STB models have high potential to further improve the simulation performance of complex cell and organ models and may represent a valuable tool in the field of next-generation in silico electrophysiology. The different modeling approaches were implemented in the simulation environment MATLAB (R2019b, MathWorks Inc.).
HH model. For the HH model, differential equations of activation and inactivation gates were solved numerically by the Forward Euler method according to Ranjan et al. [15], using a step size of ∆t = 1.10 −4 s. HMM model. The parametrization of the HMM was based on the averaged activation (n = 60) and deactivation (n = 37) data by a particle swarm optimization (PSO) algorithm (swarm size: 600; number of iterations: 10,000; function tolerance: 1 × 10 −6 ) using the Global Optimization Toolbox (MathWorks Inc.). Defining P S i ,k as the fraction of channels in a specific state S i at time-step k, the time evolution of the system could be described by the following set of autonomous difference equations: The system was simulated with the MATLAB lsim function (MathWorks Inc.) over the entire simulation protocol, with a step size of dt = 5 × 10 −7 . The output vector was defined as c T = 0 0 0 0 1 0 0 0 , to obtain the fraction of channels in the open state P O,k for each time-step k. The optimization defined the best choice for the voltage-dependent forward (α, λ, σ) and backward state transition rates (β, η, ε) and the constant state transition rates c, d, m, k, x, and y as well as the number of ion channels (N C Kv1.1 ) by fitting the resulting macroscopic (I model ) current to the measured whole-cell current (I data ): Measured activation curves between 10.4 ms and 601.7 ms were considered for parameterization. For the deactivation curves, again, only the test pulses starting from 401.2 ms to 598.2 ms were used, excluding the depolarization pulse. In addition, to account for the lower number of voltage-levels of the deactivation protocol used, the deactivation curves were weighted for adequate consideration by a factor of 3.5.
To stochastically model the opening and closing of the single Kv1.1 ion channels for model validation, the hmmgenerate function (MathWorks Inc.) was used to generate a random sequence of states. The transition probability matrix T was defined as follows: Occupancy was summarized for closed, open, and inactivated states at each time step for corresponding fractional occupancy plots.
STB model. The main advantage of this approach is that the model equations (see Equation (7)) do not have to be solved numerically. Instead, the identified model can be exported into the workspace of MATLAB where the obtained model can be further analyzed, linearized, or inserted into Simulink for a further application and simulation. The dynamic behavior of the ion channel is finally characterized by the transfer function and input/output nonlinearities.
In detail, we used the MATLAB Control System toolbox and the System Identification App. Note that a general system identification methodology contains key elements such as the experiment design, experiment, data preprocessing, fitting model to data, model validation, and audit [42].
For the STB model development, we assumed that a ramp stimulus protocol was a proper stimulus for the system identification. Using the System Identification App, we imported the time domain data and performed some data operations, including filtering, removing means, or transforming the data. Finally, we were able to define a mathematical model of the system represented as a nonlinear polynomial transfer function in state space. After analyzing different model concepts, we decided to use a nonlinear HW model because this model best fit the experimental data. After the model was validated, it could be exported to the MATLAB workspace and inserted into Simulink.
It should be noted that the HW model obtained with the System Identification App was represented in discrete-Z-domain-, which we then converted to the continuous domain since the opening/closing of ion channels is a continuous-time process. Finally, the HW model developed in MATLAB was fully functional for a further analysis, synthesis, and simulations of the Kv1.1 dynamics.
with g x = g x y z .
In the Hodgkin-Huxley formalism, the macroscopic ion conductance g x is described by various gates, controlling the flow of ions through the membrane. Each gate contains several independent gating particles z, which change between the open and closed positions, depending on the membrane potential. The gating variable y represents the probability of a single gating particle being in the open state. For several independent gating particles z, the probability of the entire gate being open, is given by y z .
The movement of gating particles between the closed and open state can be expressed as a reversible reaction with forward and backward rates α(V) and β(V): Cells 2022, 11, x FOR PEER REVIEW 26 of 28 with = ̅̅̅ . In the Hodgkin-Huxley formalism, the macroscopic ion conductance gx is described by various gates, controlling the flow of ions through the membrane. Each gate contains several independent gating particles z, which change between the open and closed positions, depending on the membrane potential. The gating variable y represents the probability of a single gating particle being in the open state. For several independent gating particles z, the probability of the entire gate being open, is given by y z .
The movement of gating particles between the closed and open state can be expressed as a reversible reaction with forward and backward rates α(V) and β(V): where y0 is the starting point at time zero, ∞ the steady state value, and the time constant. Both, ∞ and are related to the voltage dependent rate coefficients α(V) and β(V), which can further be modeled by fitting empirical functions of the membrane potential to experimental data: As proposed by Hodgkin and Huxley [42], the best fit for the non-linear potassium conductance in the squid giant axon was achieved by assuming four independent gating particles for the activation gate n, leading to the following model equations for the potassium current IK: Assuming a large number of ion channels, the probability y of individual gating particles being in the open state can be interpreted as the fraction of gating particles in the open position. Correspondingly, the fraction of gating particles in the closed state is 1 − y. Thus, the time evolution of the gating variable y can be described by a first-order differential equation: The general form of the time evolution for y(t) to a voltage step is: where y 0 is the starting point at time zero, y ∞ the steady state value, and τ y the time constant. Both, y ∞ and τ y are related to the voltage dependent rate coefficients α(V) and β(V), which can further be modeled by fitting empirical functions of the membrane potential to experimental data: y ∞ = α y α y + β y (A6) τ y = 1 α y + β y (A7) As proposed by Hodgkin and Huxley [42], the best fit for the non-linear potassium conductance in the squid giant axon was achieved by assuming four independent gating particles for the activation gate n, leading to the following model equations for the potassium current I K : with α n = 0.01