**1. Introduction**

Nanoelectromechanical systems, nanomachines, and nanorobots represent a challenging and future-orientated research area in biomedical engineering. A nanoelectromechanical system (NEMS) is a device that combines electrical and mechanical system behavior at the nanoscale level. The NEMS has electromechanical parts that have been developed on the nanoscale level, such as sensors, actuators, controllers, and drives [1]. Nanorobots are nanoelectromechanical or nanomechatronic systems, which are innovative devices expected to have revolutionary applications in health care, cancer therapy, monitoring, and drug delivery [2]. These nanoelectromechanical systems and nanorobots represent miniaturizations of microelectromechanical systems and microrobots that travel in the human body and can be used in applications to monitor, interact, and control processes at a cellular level. In particular, these robotic applications involve sensing, control, actuation and propulsion, communications, interfacing, programming and coordination at macro, micro, and nano-levels [2].

In the field of biomedical engineering and medicine, microelectromechanical/ nanoelectromechanical systems (MEMS/NEMS) are also called bioMEMS/bioNEMS [3],

**Citation:** Šaji´c, J.L.; Langthaler, S.; Baumgartner, C. Creating a Novel Mathematical Model of the Kv10.1 Ion Channel and Controlling Channel Activity with Nanoelectromechanical Systems. *Appl. Sci.* **2022**, *12*, 3836. https://doi.org/10.3390/ app12083836

Academic Editors: Luis Gracia and Carlos Perez-Vidal

Received: 14 March 2022 Accepted: 8 April 2022 Published: 11 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

adapted based on material published in [3,9]. Atom Bacteria 1mm 100μm 1μm 100nm 1nm Hair Cell Red blood cell Chromosome Virus Antibody DNA Ribossome Molecule DNA Microarray Microelectromechanical systems (MEMS) Carbon nanotubes NEMS Nanorobots Quantum dots NEMS MEMS MACRO SYSTEMS Hearth Microscope

**Figure 1.** Scales in nature, ranging from micro-systems to nano-systems.

In this work, we chose a voltage-gated potassium ion channel as an example of a biological system at a nanoscale. Ion channels are membrane proteins the size of about 4 nm that enable the passive transport of ions through the cell membrane [10,11]. Based on the scale depicted in Figure 1, we assumed that the dynamic behavior of the ion channel could be measured and controlled with MEMS/NEMS and/or nanorobots. Depending on the gating mechanism, three main types of ion channels are classified as: voltage-gated, receptor/ligand-gated, and second messenger gated channels [12]. Cells usually express a variety of ion channel types, e.g., there are more than 300 different ion channels in an inner ear cell alone [13].

with almost unlimited practical applications, e.g., in surgery, drug delivery, or gene therapy [4–8]. Figure 1 shows a schematic nanoscale comparison of different items in nature,

We here selected a single channel from a voltage-gated ion channel family representing the biological system. Voltage-gated ion channels open (activate) and close depending on the cell membrane potential. Subtypes of voltage-gated ion channels (Nav, Cav, Kv, CLC, and Hv) are specifically selective for sodium (Na+), calcium (Ca2+), potassium (K+), chloride (Cl<sup>−</sup>), and proton (H+) channels. Voltage-gated potassium channels are referred to as the giant subtype and denoted as Kv. These Kv channels encode more than 100 human genes [14]. Potassium channels represent the most complex class of voltage-gated channels [15] and are found in basically all types of mammalian cells, such as cells in the nervous, muscular, and immune systems, among others [15–18].

Each ion channel exhibits a specific dynamic behavior with different physiological and pharmacological properties [15]. In particular, aberrant expression and functionality of various Kv channels in cancer cells have been associated with tumor development and progression. The voltage-gated Kv10.1 potassium channel, encoded by the gene KCNH1 (subfamily H member 1, known as EAG1 or Ether-à-go-go 1), is implicated in various cellular processes, including, for example, cell proliferation [15].

An analysis of Kv channel expression in human cancer cells [18] revealed an upregulation of Kv10.1 in a variety of tumors and can be found in blood, bone, brain, breast, stomach, colon, cervix cancer, and prostate cancers cells, thus providing a novel biomarker candidate and potential oncological target for cancer [5,19,20].

For single-channel modeling, we finally selected the Kv10.1 ion channel, as it represents a significant ion channel in cancer development and progression. The Kv10.1 model is based on control system theory, which is similarly based on the modeling concept of the Kv1.1 voltage-gated channel recently presented in [21].

#### **2. Methods**

For mathematical modeling of voltage-gated ion channels, the method of system identification known from control engineering was used. We considered the voltage-gated ion channel Kv10.1 as a system, object, or process in accordance with the control system theory. Kv10.1 was treated as a separate object (or system/process) within a space that can interact or connect with other elements. If a system is interpreted as a combination of elements that act together and perform specific objectives [22], then a single ion channel can be considered as a system or, in more precise terms, as an object in a control loop. As such, the channel can achieve the desired dynamics under nominal conditions and display acceptable behaviors under random conditions that deviate from the maximum prescribed boundary conditions.

System identification is a scientific methodology that is used to develop mathematical models of a dynamical system. This method is based on observed or measured data of the system [23]. A block diagram of a general system (*G(s)*) represents a transfer function with a complex variable *s* in the Laplace domain, specified with an input *c*(*t*), output *r*(*t*), measured disturbances *dc*(*t*), and unmeasured disturbances *d*(*t*) (Figure 2a). An illustration of a voltage-gated ion channel as a system is shown in Figure 2b. The block diagram of the system shows its unilateral property, where system inputs and outputs, in general, are vectors. If there are multiple inputs and multiple outputs, then the system is called a MIMO system. In the case of a single input and a single output, however, the control system is called a SISO system. For example, if a voltage-gated ion channel has a voltage stimulus as input and a measured current as output, we assume that it can be described as a SISO system.

**Figure 2.** Block diagram of a system (**a**); illustrated voltage-gated ion channel as a system (**b**).

Furthermore, we assume that the voltage-gated ion channel Kv10.1 is a linear timeinvariant system based on the experimental results given in the following sections. A linear time-invariant system is defined by

$$\begin{array}{ll} a\_{n}r^{(n)}(t) + a\_{n-1}r^{(n-1)}(t) + \dots + a\_{1}\dot{r}(t) + a\_{0}r(t) = \\ b\_{0}c(t) + b\_{1}\dot{c}(t) + \dots + b\_{n-1}c^{(n-1)}(t) + b\_{n}c^{(n)}(t), \quad n \ge m. \end{array} \tag{1}$$

where *r*(*t*) is the system's output and *c*(*t*) is the system's input. The transfer function, *G*(*s*), is defined as the ratio of the left Laplace-transformed output and the left Laplace-transformed input under zero initial conditions, where *s* is a complex variable in the Laplace (*S*) domain.

$$G(s) = \frac{\mathcal{L}\{r(t)\}}{\mathcal{L}\{c(t)\}} = \frac{\mathcal{R}(s)}{\mathcal{C}(s)} = \frac{b\_0 + b\_1s + \dots + b\_{n-1}s^{n-1} + b\_ns^n}{a\_ns^n + a\_{n-1}s^{n-1} + \dots + a\_1y(t) + a\_0y(t)} \tag{2}$$

$$= \frac{\sum\_{0}^{m} b\_m s^m}{\sum\_{0}^{n} a\_ns^n}$$

The application of the transfer function concept makes it possible to represent the system's dynamics in the form of algebraic equations. If the highest power of *s* in the denominator of a transfer function is equal to *n*, the system is called an *n*th-order system [22]. The transfer function of a system is thus a mathematical model of this system. If the system's transfer function is unknown, we can estimate it based on known inputs and measured outputs for the system. This resulting transfer function provides a complete description of the system's dynamic behavior [22,24,25].

A transfer function can be determined experimentally based on measured input and output signals if a mathematical model cannot be developed using physical correlations and equations. Therefore, the process of developing a model from measured input and output data is known as system identification [23,26]. The main elements in the system identification process cycle are experimental design, the experiment itself, data preprocessing, fitting the model to the data, testing the model structure, validating and auditing the model [26].

The process of system identification has been greatly simplified due to the availability of modern software tools, such as MATLAB or LabVIEW. In this work, we used the System Identification toolbox provided by MATLAB [27,28]. This toolbox includes tools to identify a modeling approach and to define the model structure properties. Measured data can be imported as time-domain or frequency-domain data, and functions for preprocessing, representing, filtering, and estimating the model can be applied. Furthermore, a model can be estimated as a transfer function, a state-space, process, polynomial, or nonlinear model by the System Identification toolbox.

#### **3. Results**

We considered experimental patch-clamp data [15] from CHO (Chinese hamster ovary) host cells, stably expressing rat Kv10.1 channels, which exhibit a highly comparable electrophysiological behavior to human Kv10.1 channels [29]. The data provided comprises voltage-clamp measurements at three different temperature levels of 15, 25, and 35 ◦C. The dynamic behavior of the voltage-gated ion channels displays three main characteristics: activation, deactivation, and inactivation (Figure 3). In the case of the depolarization of the cell membrane, the Kv channel transits from the resting (closed) to the active (open) state. During prolonged depolarization, Kv channels switch to an inactivated state [30]. Figure 3 shows the whole-cell current response to applied voltage-step protocols for the determination of the activation, deactivation, and inactivation characteristics of the channel. In developing the dynamic behavior of the system, the three different stimulus signals, as shown in Figure 3, were used to provoke these activities. In addition, these activation, deactivation, and inactivation characteristics could also be stimulated alternatively by ramp, action potential (AP), and recovery protocols [15].

For model development, the dynamic behavior of Kv10.1 was considered at a temperature of 25 ◦C, which corresponds to standard experimental conditions in in vitro electrophysiological studies. The activation or opening is stimulated with appropriate depolarizing input voltage step functions. Nevertheless, the Kv10.1 response corresponds to the transient response of the system. Based on the experimentally measured response to given stimuli, we concluded that Kv10.1 behaves as a first-order system. We assumed all initial conditions as equal to zero and at a nominal voltage of 70 mV. It was possible to obtain an average model based on averaged data for all cells, as given in [21] for the Kv1.1 voltage-gated ion channel. We also considered a randomly chosen active cell with highquality recording data from the database available via Channelpedia, a web-based freely accessible information management network and electrophysiology data repository [15].

The chosen cell ID is 9514 at 25 ◦C, and the host cell is CHO\_FT; the species is rat. Input and output experimental data from patch-clamp measurements are given in Figure 3.

**Figure 3.** Dynamic behavior of the voltage-gated ion channel Kv10.1. Created with BioRender.

Figure 4 shows the activation voltage step protocol to determine the opening or activation of the ion channel and the measured macroscopic ion current as input and output for system identification. We assumed a nominal operating point at an input voltage of 70 mV, based on the experimental design and input functions for inactivation and deactivation. Input-output data from the patch-clamp experiments for the nominal operating point are shown in Figure 5a. Time-domain experimental data were imported into the System Identification toolbox in MATLAB. The data were not preprocessed and, based on the measured output signal and input step signals, we adopted a first-order linear time-invariant system. A comparison between measured and simulated model output is shown in Figure 5b.

The transfer function was estimated with one pole and null zeros in the continuoustime domain, the i/o delay was fixed to zero, and we set the initial condition to zero. The initialization method implies that the algorithm is used to initialize the numerator and denominator. Algorithms applicable only for estimating continuous-time transfer functions using time-domain data are the instrument variable (IV), state variable filters (SVF), generalized Poisson moment functions (GPMF), subspace state-space estimation (N4SID), and a combination of all of the preceding (ALL) approaches [27,28]. We chose to use the initialization method ALL. This method of systems modeling is known as black-box modeling.

**Figure 4.** Input: activation voltage step protocol (**above**); output: measured macroscopic ion-current (**below**) for system identification.

**Figure 5.** Input-output data for nominal operating point at 70 mV (**a**); output: macroscopic ion current (**b**) used for system identification.

The estimated transfer function model for the assumed nominal step input of 70 mV and measured output macroscopic current in nA is given by

$$G(s) = \frac{3.687}{s + 21.22} \tag{3}$$

where units of the transfer function *G*(*s*) are ratio units of output (*s*)—the output macroscopic current (nA), and units of input (*s*)—the input voltage step (mV).

The model, Equation (3), is written in a general form of a first-order system in the Laplace domain by

$$G(s) = \frac{K}{Ts + 1} \tag{4}$$

where *K* is the gain with *K* = 0.1738 and *T* the time constant, *T* = 0.0471. The system's gain is the ratio between the input signal and the steady-state value of the output.
