Simulation Modeling of Reduced Glycosylation Effects on Potassium Channels of Mouse Cardiomyocytes

Dilated cardiomyopathy (DCM) is the third most common cause of heart failure and the primary reason for heart transplantation; upward of 70% of DCM cases are considered idiopathic. Our in-vitro experiments showed that reduced hybrid/complex N-glycosylation in mouse cardiomyocytes is linked with DCM. Further, we observed direct effects of reduced N-glycosylation on K_v gating. However, it is difficult to rigorously determine the effects of glycosylation on K_v activity, because there are multiple K_v isoforms in cardiomyocytes contributing to the cardiac excitation. Due to complex functions of K_v isoforms, only the sum of K⁺ currents (I_Ksum) can be recorded experimentally and decomposed later using exponential fitting to estimate component currents, such as I_Kto, I_Kslow, and I_Kss. However, such estimation cannot adequately describe glycosylation effects and K_v mechanisms. Here, we propose a framework of simulation modeling of K_v kinetics in mouse ventricular myocytes and model calibration using the in-vitro data under normal and reduced glycosylation conditions through ablation of the Mgat1 gene (i.e., Mgat1KO). Calibrated models facilitate the prediction of K_v characteristics at different voltages that are not directly observed in the in-vitro experiments. A model calibration procedure is developed based on the genetic algorithm. Experimental results show that, in the Mgat1KO group, both I_Kto and I_Kslow densities are shown to be significantly reduced and the rate of I_Kslow inactivation is much slower. The proposed approach has strong potential to couple simulation models with experimental data for gaining a better understanding of glycosylation effects on K_v kinetics.

1. Introduction

Heart disease is the leading cause of death globally, accounting for 23% of deaths in the U.S. in 2017 (Heron, 2019). Dilated cardiomyopathy (DCM) is the third most common cause of heart failure and the most frequent reason for heart transplantation (Weintraub et al., 2017). DCM is characterized by enlarged and weakened ventricular chambers, and it is associated with systolic and contractile dysfunction that has a high risk to heart failure, with approximately 70% of DCM cases regarded as idiopathic (Hershberger and Siegfried, 2011; Lakdawala et al., 2013; Weintraub et al., 2017). There has been consistent and increasing evidence of a link between aberrant glycosylation and heart failure (Gehrmann et al., 2003; Footitt et al., 2009; Marques-da Silva et al., 2017). Recently, we showed that reduction of hybrid/complex N-glycosylation in mouse cardiomyocytes, through ablation of the Mgat1¹ gene that encodes a critical glycosyltransferase (GlcNAcT1²) (Mgat1KO model), is sufficient to cause DCM (Ednie et al., 2019a,b). Mgat1KO mice develop DCM, heart failure, and 100% die early, likely from ventricular arrhythmias resulting in sudden cardiac death. Further, Mgat1KO ventricular myocytes demonstrated altered electromechanical functions, including excitation-contraction (EC) coupling, are consistent with observed changes in electrical signaling caused by acute and downstream (disease-related) effects on voltage-gated ion channel (VGIC) gating and activity (Ednie et al., 2019b).

We investigated the impact of reduced hybrid/complex N-glycosylation through the lens of electrophysiology. Electrical signaling in the heart has vital functions related to intra and extracellular communication, rhythmicity of heartbeats, and provides a driving force for contraction (Koenig et al., 2011). This electrical stimulation is called the action potential (AP), which is the net transmembrane potential varying over time.

The AP is the result of orchestrated activities of various underlying ionic currents through VGICs (Grant, 2009). Even small changes in VGIC function can contribute to aberrant AP waveform and/or conduction, leading to arrhythmias. As illustrated in Figure 1, Na⁺ channels (Na_v) are responsible for AP initiation, or the depolarization phase. Ca²⁺ channels (Ca_v) drive the prolonged depolarization phase particularly in larger species such as human (“plateau”) as in Figure 1A, while in small species, the AP re-polarizes rapidly without the long duration of intracellular Ca²⁺ flux as in Figure 1B because smaller contractile force is required. Several K⁺ channel (K_v) isoforms are collectively responsible for AP deactivation, i.e., re-polarization. The human ether-a-go-go gene (hERG) K_v is mainly responsible for the rapid re-polarization after the plateau phase. However, mouse ether-a-go-go gene (mERG) is almost absent in the mouse ventricular AP, where the re-polarization is generally delineated as kinetic components of a rapidly inactivating transient outward current and delayed rectifier-type currents. Like most transmembrane proteins, VGICs are heavily glycosylated membrane proteins (Ednie and Bennett, 2011). Studies have proven that glycosylation can affect VGIC function predominantly. For example, we reported that a saturating, electrostatic effect of negatively charged sialic acids, which are typically attached to the terminal of glycan branches, significantly altered cardiac electrical signaling (Ednie et al., 2013; Ednie and Bennett, 2015). Our recent studies observed aberrant electrical signaling, e.g., prolonged APs and abnormal early re-activations, in Mgat1KO ventricular apex myocytes (Ednie et al., 2019a,b). These observations strongly support direct and disease-related effects of reduced N-glycosylation on VGIC gating and activity.

FIGURE 1

Figure 1. Diagrams of major cardiac ion channels and their roles in the action potential. (A) Human ventricular myocytes. (B) Mouse ventricular myocytes.

Although these in-vitro experiments showed changes in activities of ion channels, there are certain limitations that include: (1) Detailed ion channel kinetics are difficult to determine rigorously using in-vitro experiments alone. From the electrophysiological experiment, pathological electrical signaling can be observed through measurements of APs at the whole-cell level and currents at the ion-channel level. However, it is difficult to relate, rigorously, the impacts of VGIC gating changes to the altered AP waveform/conduction and vice versa. (2) Segregating different K⁺ currents is difficult using whole-cell recording experimental techniques only. I_Ksum is the result of joint activity of multiple isoforms with each isoform producing a different, but slightly overlapping (in activation and inactivation voltages) component of the total K⁺ current (I_Ksum) that contributes to various portions of AP re-polarization (Du et al., 2018). However, in-vitro experiments can only measure the sum of these component currents (I_Ksum) using whole-cell voltage-clamp protocols. Even pharmacologic separation of K_v current types is difficult, as the specificity of drugs for a single current type is not ideal. Hence, component K⁺ currents are usually estimated through multiple exponential fits of the decaying portion (due to channel inactivation) of the total K⁺ current (Brunet et al., 2004). Although the exponential function illustrates the shape of a component-current trace well, it cannot adequately describe kinetic dynamics of the current, nor fully and rigorously distinguish among currents produced by different K_v isoforms because of their slightly overlapping voltage-dependence of gating. As K_v have a critical function in re-polarizing the cell, it is imperative to model and examine K_v isoforms thoroughly.

Therefore, this article presents simulation modeling of dominant K⁺ currents in mouse ventricular myocytes and further calibrate their parameters using the in-vitro data under normal and reduced glycosylation conditions through ablation of the Mgat1 gene. Specifically, we propose an optimization procedure to calibrate simulation models of potassium channel isoform activity, I_Kto. Key statistics of component currents, amplitudes, and inactivation rates (τ's) are measured from in-vitro experimental data with a 4.5 s voltage-clamp protocol (see details in in-vitro Experimental Data Section in Material and Methods). The major contributions of this article are summarized as follows.

• Simulation models are integrated with in-vitro experiments to investigate the effects of reduced N-glycosylation on K_v activity of mouse myocytes. Computer simulation helps demonstrate the voltage-gating mechanism and conductance kinetics that cannot be readily available through traditional exponential fitting.

• Calibrated models facilitate the prediction of K_v characteristics at different voltages that are not observed in the in-vitro experiments. In the Mgat1KO group, both I_Kto and I_Kslow densities are shown to be significantly reduced and the rate of I_Kslow inactivation is much slower.

• Simulation modeling of cardiac myocytes is conducive to gain a better understanding of detailed K_v kinetics, as well as how reduced glycosylation through ablation of Mgat1 gene impacts K_v kinetics.

2. Materials and Methods

2.1. Glycosylation and Dilated Cardiomyopathy (DCM)

Protein glycosylation is an essential cellular process that impacts many cell functions (Marques-da Silva et al., 2017). Briefly, protein glycosylation is the sequential co-/post-translational process of attaching sugar residues (glycans) to proteins. Cardiac VGICs are heavily glycosylated proteins with upwards of 30% of their mass consisting of N- and O-linked glycans. A growing number of cardiac diseases, including DCM and hypertrophic cardiomyopathy, can present with concurrent, albeit, modest changes in glycosylation (Gehrmann et al., 2003; Footitt et al., 2009; Marques-da Silva et al., 2017). Mgat1 expression was implicated in cardiac function and shown to be down-regulated in human end-stage idiopathic DCM. Genome-wide searches identified changes in glycosylation-related gene expression in human idiopathic DCM, including glycosyltransferases (Barrans et al., 2002; Hwang et al., 2002; Yung et al., 2004); and proteomic/glycomic studies show changes in serum N-glycosylation in heart disease models and in humans with DCM risk factors (Nishio et al., 1995; Knezevic et al., 2009; Yang et al., 2015; Miura and Endo, 2016; Nagai-Okatani and Minamino, 2016). Models of DCM/heart failure were associated with subtle changes in glycosylation of proteins involved in electromechanical processes, and 20% of patients with congenital disorders of glycosylation (CDG) present with cardiac deficits, including idiopathic DCM (Gehrmann et al., 2003; Marques-da Silva et al., 2017). The data suggest a correlation between modest changes in extracellularly facing glycosylation and DCM/heart disease.

2.2. Potassium Currents Analysis and in-vitro Experimental Data

The orchestrated activity of different K_v isoforms produces unique but slightly overlapping currents, as illustrated in Figure 2A. However, it is difficult to rigorously delineate the multiple component K⁺ currents through voltage-clamp experiments alone because only the sum of K⁺ currents is measured and recorded. I_Ksum is usually decomposed mathematically by fitting the peak and decaying portion of each component current trace with an exponential function (Brunet et al., 2004). A standard exponential function, I = Ae^−t/τ + C has three parameters: amplitude A, time constant τ, and constant offset value C. t is time in millisecond. Amplitude refers to the value of the peak and τ to the time constant when current reduces by e⁻¹ ( 63%) of the peak.

FIGURE 2

Figure 2. (A) The action potential of mouse ventricular myocytes and underlying ioninc currents. (B) Predominant K⁺ currents in mouse cardiomyocytes and their shapes of whole-cell voltage clamp recordings.

As illustrated in Figure 2B, I_Ksum is often decomposed into three dominant currents: a rapidly inactivating transient outward current (I_Kto), a delayed rectifier-type current (I_Kslow), and a non-inactivating steady-state current (I_Kss) (Nerbonne and Kass, 2005). I_Kto has a high peak at the very beginning of activation and rapid inactivation. I_Kslow has a low peak and longer inactivation phase. I_Kss remains constant during the course of depolarization. A bi-exponential function with a constant component, as in Equation (1), is a mathematical form for decomposing I_Ksum into the three major K⁺ currents. In Equation (1), A_Kto and A_Kslow are the amplitudes; τ_Kto and τ_Kslow are time constants of I_Kto and I_Kslow, respectively. A_Kss is the constant current I_Kss. Although bi-exponential fitting captures essential characteristics of three major component K⁺ currents and describes the I_Ksum waveform, it cannot provide detailed kinetic dynamics of K_v isoforms (Plumlee et al., 2016).

$\begin{array}{ll}{I}_{\text{Ksum}}={\text{A}}_{Kto}{e}^{-t/{\tau }_{Kto}}+{\text{A}}_{Kslow}{e}^{-t/{\tau }_{Kslow}}+{\text{A}}_{Kss}& (1)\end{array}$

Recently, we reported electrophysiological experiment data to investigate the impact of reduced hybrid/complex N-glycosylation on left ventricular cardiomyocyte activity, through deletion of the Mgat1 gene, which encodes a critical glycosyltransferase (GlcNAcT1) (Ednie et al., 2019a,b). The detailed process of creation of Mgat1KO (Mgat1 Knock Out) strain, features of the cardiomyocyte-specific Mgat1KO strain, breeding, genotyping, and selection of wild type (WT) animals were previously described in Ednie et al. (2019a). In our previous study of K_v activities in the Mgat1KO (Ednie et al., 2019b), cells of 12 to 20-week-old mice were used and numbers of observation for WT and Mgat1KO are 35 and 38, respectively.

I_Ksum was measured using whole-cell voltage clamp protocols. Cells were held at -70 mV and then depolarized by 10 mV voltage steps from -50 to 50 mV for 4.5 s. Figure 3 shows our in-vitro experimental data about key statistics on I_Kto and I_Kslow. As shown in Figure 3A, peaks of both I_Kto and I_Kslow are reduced in the Mgat1KO group; and steady-state current I_Kss, remaining constant during the depolarization, is also reduced in the Mgat1KO group. A_Kslow is most reduced compared to WT, by 77%. Figure 3B compares the time constants τ. τ_Kslow of the Mgat1KO group is 30% larger than the WT group, which means I_Kslow is prolonged with reduced N-glycosylation. However, τ_Kto does not show significant difference.

FIGURE 3

Figure 3. In-vitro experimental results, adopted from Ednie et al. (2019b). (A) Averages and error bars of the peaks of I_Kto and I_Kslow, and steady-state amplitude of I_Kss. (B) Inactivation time constant of I_Kto and I_Kslow. Significant differences between WT and Mgat1KO at p ≤ 0.05 are indicated by an *. (n = 35 for WT and n = 38 for Mgat1KO).

2.3. In-silico Simulation Modeling

The AP is modeled as an orchestrated activity of various ionic currents as given in Equation (2). In this equation, C_m is the membrane capacitance, and there are several ionic currents: L-type calcium current (I_CaL), calcium pump current (I_p(Ca)), Na⁺/Ca²⁺ exchange current (I_NaCa), calcium background current (I_Cab), fast Na⁺ current (I_Na), background Na⁺ current (I_Nab), Na⁺/K⁺ pump current (I_NaK), fast transient outward K⁺ current (I_{Kto, f}), slower transient outward K⁺ current (I_{Kto, s}), which is essentially missing in apex myocytes, time independent K⁺ current (I_K1), slow delayed rectifier K⁺ current (I_Ks), ultrarapidly activating delayed rectifier K⁺ current (I_Kur), non-inactivating steady-state K⁺ current (I_Kss), rapidly delayed rectifier K⁺ current (I_Kr), Ca²⁺-activated Cl^- current (I_{Cl, Ca}), and stimulus current (I_stim). Although there are seven currents pertinent to K_v dynamics in this model, as mentioned earlier, three current types, I_Kto, I_Kslow, and I_Kss, are dominant K⁺ currents in a ventricular apex myocyte (also see Figure 2A). Thus, in this study, we focus on development of computer models of I_Kto and I_Kslow. I_Kss is assumed as a constant as in Equation (1). This model can be flexibly modified to include I_{Kto, s} for the septum. The details of ionic currents models and their kinetics can be found in Bondarenko et al. (2004).

$\begin{array}{l}-C_m\frac{dV}{dt}={\text{I}}_{\text{CaL}}+{\text{I}}_{\text{p}(Ca)}+{\text{I}}_{\text{NaCa}}+{\text{I}}_{\text{Cab}}\\ +{\text{I}}_{\text{Na}}+{\text{I}}_{\text{Nab}}+{\text{I}}_{\text{NaK}}+{\text{I}}_{\text{Kto},f}\\ +{\text{I}}_{\text{Kto},s}+{\text{I}}_{\text{K}1}+{\text{I}}_{\text{Ks}}+{\text{I}}_{\text{Kur}}+{\text{I}}_{\text{Kss}}+{\text{I}}_{\text{Kr}}\\ +{\text{I}}_{\text{Cl},Ca}+{\text{I}}_{\text{stim}}& (2)\end{array}$

There are two major modeling schemes for VGICs: Markov models and Hodgkin-Huxley models. While Markov modeling of ionic currents is popular such as O'Hara models (O'Hara et al., 2011; O'Hara and Rudy, 2012), Hodgkin-Huxley modeling is commonly used to formulate K⁺ currents for various species, for example, not only for mouse ventricular myocytes (Bondarenko et al., 2004) but also for rabbit (Mahajan et al., 2008) and human ventricular myocytes (ten Tusscher et al., 2004). I_Kto and I_Kslow are designed using Hodgkin-Huxley modeling as given in Equations (3) and (4), respectively, based on Bondarenko et al. (2004). Hodgkin-Huxley models have two state variables for describing activation and inactivation. In Equation (3), a_to and i_to are two state variables of activation and inactivation, respectively. G_to is the maximum conductance of K_v4.2, and (V−E_K) is the difference of voltage V and the K⁺ Nernst potential E_K. Similarly, the I_Kslow model is given as Equation (4), which also include two state variables for activation and inactivation gates. Relevant parametric equations of transition rates and time constants for the gating mechanisms are given below for each current equation.

$\text{}{\text{I}}_{\text{Kto}}=G_{to}{a}_{to}^3{i}_{to}(V-E_K)(3){\text{I}}_{\text{Kslow}}=G_r{a}_r{i}_r(V-E_K)(4)$

$\begin{array}{l}\frac{d{a}_{to}}{dt}={\alpha }_a(1-a_{to})-{\beta }_a{a}_{to}\frac{a_r}{dt}=\frac{a_{ss}-a_r}{{\tau }_a}\\ \frac{d{i}_{to}}{dt}={\alpha }_i(1-i_{to})-{\beta }_i{i}_{to}\frac{i_r}{dt}=\frac{i_{ss}-i_r}{{\tau }_i}\\ {\alpha }_a=0.18064e^{0.03577(V+x_1)}{a}_{ss}=\frac{1}{1+e^{-(V+x_1)/x_2}}\\ {\beta }_a=0.395e^{-0.06237(V+x_2)}{i}_{ss}=\frac{1}{1+e^{(V+x_3)/x_4}}\\ {\alpha }_i=\frac{0.000152e^{-(V+x_3)/x_4}}{0.067083e^{-(V+x_5)/x_4}+1}{\tau }_a=0.493e^{-0.0629V}+x_5\\ {\beta }_i=\frac{0.00095e^{(V+x_6)/x_7}}{0.051335e^{(V+x_6)/x_7}+1}{\tau }_i=x_6-\frac{170}{1+e^{(V+x_7)/x_8}}\end{array}$

2.4. Model Calibration

In the literature, optimization procedures have been proposed to calibrate simulation models with in-vitro experimental data. Du et al. (2015) integrates statistical metamodeling with in-vitro data to investigate the effects of reduced glycosylation on Na⁺ channels in mouse cardiomyocytes. Also, non-linear optimization methods show promising results for calibrating in-silico models of K_v channels, e.g., hERG channels (Du et al., 2013) and K_v isoforms in mouse cardiomyocytes (Du et al., 2018). Note that characteristic curves such as steady-state activation (SSA), steady-state inactivation (SSI), and time constant of inactivation at different voltage were used in these studies to evaluate goodness-of-fit as an objective function of optimization to minimize discrepancies between simulation results and in-vitro experimental data. However, these characteristic curves are not always available. In our Mgat1KO data, currents were too small at early-activating voltages to reliably calculate the characteristic curves.

This article presents a new method that calibrates the simulation models with sparse data when rich information of the characteristic curves are not available. To be specific, the proposed method directly fits models to amplitudes and τ's, which are the most critical characteristics to describe component K⁺ current traces. Equation (5) shows the objective function of the suggested optimization procedure, where $\widehat{A_i}$ and $\widehat{{\tau }_i}$ are estimations of amplitude and time constant from the simulation model. Potassium channels from mouse ventricular myocytes are modeled and parameterized based on partial differential equations. Then, full factorial designs are used to discover sensitive variables from the potential set of parameters. Further, a genetic algorithm-based heuristic optimization method is developed to calibrate the models to the data.

$\begin{array}{ll}min(|{\text{A}}_i-\widehat{{\text{A}}_i}|+|{\tau }_i-\widehat{{\tau }_i}|),i\in \left\{\text{Kto},\text{ Kslow}\right\}& (5)\end{array}$

In in-silico experiments, it is critical to identify a sparse subset of sensitive variables (or model parameters to be calibrated) to reduce computational burden and improve modeling accuracy. The curse of dimensionality causes the dramatic surge of required computing resources when the number of variables increases and counter-intuitive geometric properties, making the learning procedure difficult. While a sparse subset of variables restricts the model flexibility to fit data, the opposite also has a detrimental effect on modeling accuracy. We used two-level factorial designs to perform sensitivity analysis and screen the variables that impact the model outputs of our interest. Variables were adjusted at two different levels to assess their effects on model outputs and K_v characteristics. For I_Kto, x₁, x₂, x₃, x₆, x₇, and G_to are selected. For I_Kslow, x₁, x₂, x₃, x₄, x₆, and G_r are selected.

A metaheuristics optimization method is developed based on the genetic algorithm (GA) to calibrate the computer models. GA mimics the evolution of genes inspired by natural selection, and its procedure is described in Algorithm 1 (Reali et al., 2017). In the standard GA, a new population is constructed from superior solutions of the current population and new solutions are reproduced in Breeding step by so-called crossover. The crossover is a process of mixing up elements of superior solutions in the current population so that this process is also called reproduction. This procedure is expected to discover better candidates for combinatorial optimization problems by searching combinations of elite solutions. Although crossover is an intuitive searching method for combinatorial or discrete optimization problems, exploring combinations is not appropriate for continuous variables. Hence, in this article, we presented a new method named self-breeding genetic algorithm that reproduces new generations without crossover but directly breeds them from each superior solution as described in Algorithm 2.

ALGORITHM 1

Algorithm 1. Standard Genetic Algorithm.

ALGORITHM 2

Algorithm 2. Self-breeding Genetic Algorithm.

3. Results

The self-breeding GA was applied to a pair of mean values of A_i and τ_i in Figure 3, i∈{to, Kslow}, for each group, WT and Mgat1KO. This calibration process was repeated 30 times. The algorithm stopped when discrepancies of amplitude and time constant were less than tolerances. Tolerances for amplitude and time constant were (0.1, 1.0) for I_Kto and (0.1, 5.0) for I_Kslow, which are smaller than the standard error of mean values. Table 1 shows means and standard errors of calibrated parameters found by the suggested optimization process. The simulation models and optimization process were implemented in MATLAB 2020a, and the codes are available at https://github.com/haedong31/20_Kv_simulation.

TABLE 1

Table 1. Calibrated model parameters.

Figure 4 shows the current traces generated by the calibrated simulation models and exponential fitting with the mean amplitudes and time constants from in-vitro experiment data. Simulation models are shown to be compatible with results of exponential fitting in the electrophysiological data. The solid gray lines represent simulated current traces with 30 repetitions. Except for I_Kslow of the Mgat1KO, it is difficult to distinguish the two traces through a visual inspection. For this one, there was a slight dissimilarity during the early decay phase. Also, before applying the clamp voltage, it showed instability that currents had small positive values under the holding potential. Due to the small peak but large time constant, the trace of the Mgat1KO is flattened compared to a standard exponential function, which sharply decreases after the peak. This change is expected to show modeling differences between WT and Mgat1KO.

FIGURE 4

Figure 4. Comparison of current traces between calibrated simulations models (gray line) and benchmark exponential fitting. (A) I_Kto for WT. (B) I_Kslow for WT. (C) I_Kto for Mgat1KO. (D) I_Kslow for Mgat1KO. Each plot includes 30 replications of simulations results.

To further investigate differences between the WT and Mgat1KO groups, current traces were simulated at voltages between -60 and 50 mV in 10-mV increments as shown in Figure 5. Simulated current traces show different patterns as voltage changes between two groups. In Figures 5C,D, the amplitude declines rapidly as voltage decreases for Mgat1KO, while reducing evenly in WT as in Figures 5A,B. Because of the rapid reduction, I_Kslow traces for Mgat1KO below 0 mV depolarizations are negligible. In addition, from the perspective of decay phase, the current decreases slowly for I_Kslow in the Mgat1KO group. It is predicted that there are differences in conductance values between WT and Mgat1KO.

FIGURE 5

Figure 5. Predicted current traces. (A) I_Kto for WT. (B) I_Kslow for WT. (C) I_Kto for Mgat1KO. (D) I_Kslow for Mgat1KO. Clamp voltages of 4.5 s were applied to from -50 to +50 mV by 10-mV increments from the holding potential -70 mV.

Figure 6 further demonstrates the differences in predicted current density to voltage relationship between WT and Mgat1KO. For both I_Kto and I_Kslow, channels in Mgat1KO are less active at a given depolarized activation voltage compared to WT. The difference in the current peaks between the two groups is more prominent at depolarizations greater than -10 mV. This gap is more significant in I_Kslow than I_Kto. At all voltages, the Mgat1KO I_Kslow is smaller than the WT I_Kslow. This is consistent with in-vitro experimental results in which the average amplitude of I_Kto and I_Kslow were significantly reduced with N-glycosylation perturbation (in the Mgat1KO), and the reduction in Mgat1KO I_Kslow was bigger than for Mgat1KO I_Kto.

FIGURE 6

Figure 6. Predicted current density to voltage relationship under the WT and Mgat1KO conditions. (A) I_Kto. (B) I_Kslow. Voltage steps range from -60 to 50 mV by 10-mV increments.

K_v has unique inactivation kinetics that makes it difficult to distinguish component K⁺ currents. To investigate its inactivation gating kinetics, inactivation time constants at various voltages were simulated as shown in Figure 7. τ_Kto does not show a significant difference between Mgat1KO and WT when the depolarization is greater than -20 mV. There are gaps in the case where voltages are less than -20 mV, but the uncertainty represented by the error bars are too wide to make reliable predictions. τ_Kslow shows consistent differences, with the Mgat1KO inactivating significantly more slowly than WT at all voltages. Figure 8 provides further information of the steady-state inactivation (SSI) rate. The SSI relationships are not significantly different between Mgat1KO and WT for either current type. These data, the longer transition times from open to inactivated state and the lack of a shift in the voltage dependence of K_v distribution between open and inactivated states for less glycosylated K_v (Mgat1KO), are consistent with our previous studies(Schwetz et al., 2010; Ednie and Bennett, 2015; Du et al., 2018).

FIGURE 7

Figure 7. Predicted inactivation time constants to voltage relationship under the WT and Mgat1KO conditions. (A) I_Kto. (B) I_Kslow. Voltage steps range from -60 to 50 mV by 10-mV increments.

FIGURE 8

Figure 8. Predicted steady-state inactivation to voltage relationship under the WT and Mgat1KO conditions. (A) I_Kto. (B) I_Kslow. Voltage steps range from -60 to 50 mV by 10-mV increments.

K_v plays a vital role in the repolarization phase. In-silico modeling results show a reduction in the magnitudes of both I_Kto and I_Kslow currents. In particular, I_Kslow shows a significant decrease in the current magnitude and a slower inactivation rate. These results are consistent with our experimental data as summarized in Figure 3. Figure 9A shows in-vitro experimental data of AP traces from WT and Mgat1KO (Ednie et al., 2019b). Mgat1KO shows prolonged APs in the repolarization phase compared to WT. The magnitude reduction and slow inactivation of I_Kslow in Mgat1KO contribute to such AP prolongations. Figure 9B shows the prediction results from our in-silico model that are calibrated for I_Kto and I_Kslow under WT and Mgat1KO conditions. Figure 9B shows how Mgat1KO alters the K⁺ currents and causes the AP to be prolonged. The in-silico modeling results corroborate our in-vitro experimental data.

FIGURE 9

Figure 9. AP traces from WT and Mgat1KO. (A) Experimental data (Ednie et al., 2019b). (B) Model prediction. APs from the Mgat1KO group are significantly prolonged compared to WT in the repolarization phase.

4. Discussion

This article has developed a self-breeding GA method to calibrate simulation models of K⁺ channel of mouse ventricular apex myocytes based on key statistics and raw data obtained from in-vitro voltage-clamp experiments (Ednie et al., 2019b). In-silico simulation models allow for the investigation of underlying dynamics of observed current and to make inferences about different experimental conditions that were not yet or cannot be conducted.

Notably, ventricular myocytes consist of several K_v isoforms that lead to multiple component K⁺ currents. These component currents collectively produce the total I_Ksum that can be recorded during in-vitro experiments. Typically, the sum of potassium currents is approximated by fitting multiple exponential functions, each of which is used to estimate the component current. However, traditional exponential functions cannot adequately describe detailed kinetics.

Simulation models were calibrated by comparing their simulated current traces with those derived from the exponential decomposition of I_Ksum. The simulation models enable the prediction of K_v characteristics at different voltages that may not be observed in the electrophysiological experiments. In addition to K_v characteristics at 50 mV, simulation models show that changes of amplitudes, time constants, and SSI with different voltages are consistent with previous studies (Ednie et al., 2019b).

Experimental results show that the proposed optimization method effectively calibrate simulation models of K_v's with in-vitro experimental data. These simulations models allow for the prediction of component K⁺ currents at different experimental conditions. In the Mgat1KO group, both I_Kto and I_Kslow densities are significantly reduced and the rate of I_Kslow inactivation is much slower, which are consistent with experimental observations.

This study couples simulation models with in-vitro experiments to investigate the effects of reduced N-glycosylation on K_v activity. A new method is used to delineate and examine multiple potassium currents. The computational simulation helps demonstrate the voltage-gating mechanism and conductance kinetics, thereby gaining a better understanding of glycosylation effects on K_v kinetics. This article does not include the 25s depolarization, which often helps distinguish I_Kslow1 (K_v1.5 activity) and I_Kslow2 (K_v2.1 activity). Our ongoing experiments focus on voltage-clamp experiments on 25-s periods. In the future work, new sets of in-vitro data will be analyzed via the in-silico modeling approach for the Mgat1KO effects on a variety of ionic currents.

Data Availability Statement

The datasets generated and analyzed for this study are included in the article material and can be found in the GitHub repository https://github.com/haedong31/20_Kv_simulation.

Ethics Statement

Animals were used and cared for as outlined by the NIH's Guide for the Care and Use of Laboratory Animals. All animal protocols were reviewed and approved by the Wright State University Institutional Animal Use and Care Committee.

Author Contributions

HK: methodology, experimental analysis, and manuscript preparation. HY: methodology, experimental analysis, conceptualization, supervision, project administration, and manuscript preparation. AE and EB: data curation, clinical interpretation, and manuscript preparation. All authors contributed to the article and approved the submitted version.

Funding

This work was supported in part by the National Science Foundation (MCB-1856132 to HY) and (MCB-1856199, IOS-1146882, and IOS-1660926 to EB), an American Heart Association, Greater Southeast Affiliate Grant-In-Aid (14GRNT20450148 to EB), and Postdoctoral Fellowship (15POST25710010 to AE). HK and HY would also like to thank the NSF I/UCRC Center for Health Organization Transformation (CHOT) award NSF IIP-1624727 for the support of their research work.

Author Disclaimer

Any opinions, findings, or conclusions expressed in this article are those of the authors and do not necessarily reflect the views of the sponsor.

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.

Footnotes

1. ^Mannosyl (α-1,3-)-Glycoprotein β-1,2-N-acetylglucosaminyltransferase.

2. ^UDP-GlcNAc:α-3-D-mannoside-β1,2-N-acetylglucosaminyltransferase.

References

Barrans, J. D., Allen, P. D., Stamatiou, D., Dzau, V. J., and Liew, C.-C. (2002). Global gene expression profiling of end-stage dilated cardiomyopathy using a human cardiovascular-based cdna microarray. Am. J. Pathol. 160, 2035–2043. doi: 10.1016/S0002-9440(10)61153-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bondarenko, V. E., Szigeti, G. P., Bett, G. C., Kim, S.-J., and Rasmusson, R. L. (2004). Computer model of action potential of mouse ventricular myocytes. Am. J. Physiol. Heart Circulatory Physiol. 287, H1378–H1403. doi: 10.1152/ajpheart.00185.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunet, S., Aimond, F., Li, H., Guo, W., Eldstrom, J., Fedida, D., Yamada, K. A., and Nerbonne, J. M. (2004). Heterogeneous expression of repolarizing, voltage-gated k+ currents in adult mouse ventricles. J. Physiol. 559, 103–120. doi: 10.1113/jphysiol.2004.063347

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, D., Yang, H., Ednie, A. R., and Bennett, E. S. (2015). Statistical metamodeling and sequential design of computer experiments to model glyco-altered gating of sodium channels in cardiac myocytes. IEEE J. Biomed. Health Inf. 20, 1439–1452. doi: 10.1109/JBHI.2015.2458791

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, D., Yang, H., Ednie, A. R., and Bennett, E. S. (2018). In-silico modeling of the functional role of reduced sialylation in sodium and potassium channel gating of mouse ventricular myocytes. IEEE J. Biomed. Health Inf. 22, 631–639. doi: 10.1109/JBHI.2017.2664579

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, D., Yang, H., Norring, S. A., and Bennett, E. S. (2013). In-silico modeling of glycosylation modulation dynamics in herg ion channels and cardiac electrical signals. IEEE J. Biomed. Health Inf. 18, 205–214. doi: 10.1109/JBHI.2013.2260864

PubMed Abstract | CrossRef Full Text | Google Scholar

Ednie, A. R., and Bennett, E. S. (2011). Modulation of voltage-gated ion channels by sialylation. Comprehensive Physiol. 2, 1269–1301. doi: 10.1002/cphy.c110044

PubMed Abstract | CrossRef Full Text | Google Scholar

Ednie, A. R., and Bennett, E. S. (2015). Reduced sialylation impacts ventricular repolarization by modulating specific k+ channel isoforms distinctly. J. Biol. Chem. 290, 2769–2783. doi: 10.1074/jbc.M114.605139

PubMed Abstract | CrossRef Full Text | Google Scholar

Ednie, A. R., Deng, W., Yip, K.-P., and Bennett, E. S. (2019a). Reduced myocyte complex n-glycosylation causes dilated cardiomyopathy. FASEB J. 33, 1248–1261. doi: 10.1096/fj.201801057R

PubMed Abstract | CrossRef Full Text | Google Scholar

Ednie, A. R., Horton, K.-K., Wu, J., and Bennett, E. S. (2013). Expression of the sialyltransferase, st3gal4, impacts cardiac voltage-gated sodium channel activity, refractory period and ventricular conduction. J. Mol. Cell. Cardiol. 59, 117–127. doi: 10.1016/j.yjmcc.2013.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Ednie, A. R., Parrish, A. R., Sonner, M. J., and Bennett, E. S. (2019b). Reduced hybrid/complex n-glycosylation disrupts cardiac electrical signaling and calcium handling in a model of dilated cardiomyopathy. J. Mol. Cell. Cardiol. 132, 13–23. doi: 10.1016/j.yjmcc.2019.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Footitt, E., Karimova, A., Burch, M., Yayeh, T., Dupre, T., Vuillaumier-Barrot, S., et al. (2009). Cardiomyopathy in the congenital disorders of glycosylation (cdg): a case of late presentation and literature review. J. Inherited Metabolic Dis. Official J. Soc. Study Inborn Errors Metab. 32, 313–319. doi: 10.1007/s10545-009-1262-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gehrmann, J., Sohlbach, K., Linnebank, M., Böhles, H.-J., Buderus, S., Kehl, H. G., et al. (2003). Cardiomyopathy in congenital disorders of glycosylation. Cardiol. Young 13, 345–351. doi: 10.1017/S1047951103000702

CrossRef Full Text | Google Scholar

Grant, A. O. (2009). Cardiac ion channels. Circulat. Arrhythmia Electrophysiol. 2, 185–194. doi: 10.1161/CIRCEP.108.789081

PubMed Abstract | CrossRef Full Text | Google Scholar

Heron, M. (2019). Deaths: leading causes for 2017. Nat. Vital Stat. Rep. 68, 1–77. Available online at: https://www.cdc.gov/nchs/data/nvsr/nvsr68/nvsr68_06-508.pdf

Google Scholar

Hershberger, R. E., and Siegfried, J. D. (2011). Update 2011: clinical and genetic issues in familial dilated cardiomyopathy. J. Am. Coll. Cardiol. 57, 1641–1649. doi: 10.1016/j.jacc.2011.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, J.-J., Allen, P. D., Tseng, G. C., Lam, C.-W., Fananapazir, L., Dzau, V. J., et al. (2002). Microarray gene expression profiles in dilated and hypertrophic cardiomyopathic end-stage heart failure. Physiol. Gen. 10, 31–44. doi: 10.1152/physiolgenomics.00122.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Knezevic, A., Polasek, O., Gornik, O., Rudan, I., Campbell, H., Hayward, C., et al. (2009). Variability, heritability and environmental determinants of human plasma n-glycome. J. Proteome Res. 8, 694–701. doi: 10.1021/pr800737u

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenig, X., Dysek, S., Kimbacher, S., Mike, A. K., Cervenka, R., Lukacs, P., et al. (2011). Voltage-gated ion channel dysfunction precedes cardiomyopathy development in the dystrophic heart. PLoS ONE 6, e20300. doi: 10.1371/journal.pone.0020300

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakdawala, N. K., Winterfield, J. R., and Funke, B. H. (2013). Dilated cardiomyopathy. Circul. Arrhythmia Electrophysiol. 6, 228–237. doi: 10.1161/CIRCEP.111.962050

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L.-H., et al. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410. doi: 10.1529/biophysj.106.98160

PubMed Abstract | CrossRef Full Text | Google Scholar

Marques-da Silva, D., Francisco, R., Webster, D., dos Reis Ferreira, V., Jaeken, J., and Pulinilkunnil, T. (2017). Cardiac complications of congenital disorders of glycosylation (cdg): a systematic review of the literature. J. Inherited Metabolic Dis. 40, 657–672. doi: 10.1007/s10545-017-0066-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Miura, Y., and Endo, T. (2016). Glycomics and glycoproteomics focused on aging and age-related diseases–glycans as a potential biomarker for physiological alterations. Biochimica et Biophysica Acta (BBA)-Gen. Subjects 1860, 1608–1614. doi: 10.1016/j.bbagen.2016.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagai-Okatani, C., and Minamino, N. (2016). Aberrant glycosylation in the left ventricle and plasma of rats with cardiac hypertrophy and heart failure. PLoS ONE 11, e0150210. doi: 10.1371/journal.pone.0150210

PubMed Abstract | CrossRef Full Text | Google Scholar

Nerbonne, J. M., and Kass, R. S. (2005). Molecular physiology of cardiac repolarization. Physiol. Rev. 85, 1205–1253. doi: 10.1152/physrev.00002.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishio, Y., Warren, C. E., Buczek-Thomas, J. A., Rulfs, J., Koya, D., Aiello, L. P., et al. (1995). Identification and characterization of a gene regulating enzymatic glycosylation which is induced by diabetes and hyperglycemia specifically in rat cardiac tissue. J. Clin. Investigat. 96, 1759–1767. doi: 10.1172/JCI118221

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Hara, T., and Rudy, Y. (2012). Arrhythmia formation in subclinical (“silent”) long qt syndrome requires multiple insults: quantitative mechanistic study using the kcnq1 mutation q357r as example. Heart Rhythm 9, 275–282. doi: 10.1016/j.hrthm.2011.09.066

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput. Biol. 7, e1002061. doi: 10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Plumlee, M., Joseph, V. R., and Yang, H. (2016). Calibrating functional parameters in the ion channel models of cardiac cells. J. Amer. Stat. Assoc. 111, 500–509. doi: 10.1080/01621459.2015.1119695

CrossRef Full Text | Google Scholar

Reali, F., Priami, C., and Marchetti, L. (2017). Optimization algorithms for computational systems biology. Front. Appl. Math. Stat. 3, 6. doi: 10.3389/fams.2017.00006

CrossRef Full Text | Google Scholar

Schwetz, T. A., Norring, S. A., and Bennett, E. S. (2010). N-glycans modulate kv1. 5 gating but have no effect on kv1. 4 gating. Biochimica et Biophysica Acta (BBA)-Biomembranes 1798, 367–375. doi: 10.1016/j.bbamem.2009.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H., Noble, D., Noble, P.-J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. Heart Circulatory Physiol. 286, H1573–H1589. doi: 10.1152/ajpheart.00794.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Weintraub, R. G., Semsarian, C., and Macdonald, P. (2017). Dilated cardiomyopathy. Lancet 390, 400–414. doi: 10.1016/S0140-6736(16)31713-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Chen, L., Sun, S., Shah, P., Yang, W., Zhang, B., et al. (2015). Glycoproteins identified from heart failure and treatment models. Proteomics 15, 567–579. doi: 10.1002/pmic.201400151

PubMed Abstract | CrossRef Full Text | Google Scholar

Yung, C. K., Halperin, V. L., Tomaselli, G. F., and Winslow, R. L. (2004). Gene expression profiles in end-stage human idiopathic dilated cardiomyopathy: altered expression of apoptotic and cytoskeletal genes. Genomics 83, 281–297. doi: 10.1016/j.ygeno.2003.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar