Computational modeling reveals key contributions of KCNQ and hERG currents to the malleability of uterine action potentials underpinning labor.

The electrical excitability of uterine smooth muscle cells is a key determinant of the contraction of the organ during labor and is manifested by spontaneous, periodic action potentials (APs). Near the end of term, APs vary in shape and size reflecting an ability to change the frequency, duration and amplitude of uterine contractions. A recent mathematical model quantified several ionic features of the electrical excitability in uterine smooth muscle cells. It replicated many of the experimentally recorded uterine AP configurations but its limitations were evident when trying to simulate the long-duration bursting APs characteristic of labor. A computational parameter search suggested that delayed rectifying K(+) currents could be a key model component requiring improvement to produce the longer-lasting bursting APs. Of the delayed rectifying K(+) currents family it is of interest that KCNQ and hERG channels have been reported to be gestationally regulated in the uterus. These currents exhibit features similar to the broadly defined uterine IK1 of the original mathematical model. We thus formulated new quantitative descriptions for several I(KCNQ) and I(hERG). Incorporation of these currents into the uterine cell model enabled simulations of the long-lasting bursting APs. Moreover, we used this modified model to simulate the effects of different contributions of I(KCNQ) and I(hERG) on AP form. Our findings suggest that the alterations in expression of hERG and KCNQ channels can potentially provide a mechanism for fine tuning of AP forms that lends a malleability for changing between plateau-like and long-lasting bursting-type APs as uterine cells prepare for parturition.

Introduction

The timely onset and maintenance of regular contractions of the uterus are necessary features for ensuring successful parturition and safe delivery of a baby and placenta. These contractions are driven by episodic, spontaneous myogenic action potentials (APs) that exhibit a broad spectrum of form and duration, the variability of which is likely to be beneficial in facilitating the co-ordination of uterine contractile effort during labor. For example, as pregnancy progresses uterine APs become more frequent and regular, the duration of the APs progressively lengthens [1] and they are often manifested as bursting APs that can last for many seconds to minutes [2], thus contributing to a greater time-averaged force of uterine contractions at parturition. Although more frequently shown in the recordings of multi-cellular uterine tissues, these APs have also been reported in single uterine smooth muscle (myometrial) cells isolated from near-labor, and in-labor, animals and humans [2], [3]. The quest for a deeper understanding of the mechanisms that regulate these changes in uterine AP forms would benefit from computational biology approaches. A recently established mathematical model for uterine smooth muscle cell (USMC) electrical excitability is a useful starting point for these considerations [4]. The electrogenic components of this model were largely based on published data from uterine and other smooth muscles. These components were validated by experimental data and the model reproduced several published types of AP forms (spike, plateau and short bursts of spikes) of variable duration ≈0.5 – 10 secs. However, despite showing these positive simulation performances, this model was limited in its ability to computationally reproduce all uterine AP behaviors, in particular, the very long duration (tens of sec to mins) bursting-type APs that have been experimentally noted in uterine smooth muscle cells from late pregnancy [1]–[3]. Also, another related problem is that the membrane voltage failed to repolarize upon cessation of the stimulation when it had remained at a depolarized level for too long, such as during attempts to simulate a prolonged AP. These limitations restrict the use of the cell model to simulating uterine AP responses of shorter duration.

Our aims in this study, therefore, were two-fold. First, to identify membrane currents that are essential for forming the characteristics of long duration bursting-type APs with full repolarization capacity that are, presently, lacking in the uterine cell model. Second, to perform simulation experiments to explore the possible mechanistic roles of these components in the malleability of uterine cell excitability.

Results and Discussion

Insight from one-at-a-time sensitivity measures

The original USMC model [4] (the electrogenic components are depicted in Figure 1A) reliably simulated uterine bursting-type and plateau-type APs of short duration, examples of which are given in Figure 1B-C. However, to computationally reproducing APs of long duration bursting was a challenge (Figure 1D). We, therefore, performed a simple computational parameter sensitivity analysis to assess which kinetic properties would have the ability to improve the functionality of the model.

Figure 1

Electrogenic components of the initial USMC model.

A, Schematic representation of the 14 ion channels/exchangers that contributed to the initial USMC model of myometrial electrical excitability [4]. B-D, this USMC model [4] reproduced different AP forms of short duration but it could not simulate long-lasting APs. B, short bursting-type AP evoked by a 2 s current clamp stimulus (Ist) of -0.2 pA pF−1. C, short plateau-type AP evoked by a 2 s Ist of −0.5 pA pF−1. D, an attempt to simulate a long duration AP with a 50 s Ist of −0.25 pA pF−1 resulted in a plateau-like AP that failed to repolarize at the end of the stimulus. All three simulated APs were initiated from the same initial conditions.

Parameters such as the maximal conductances (), the half activations and inactivations (V 0.5), the slope factors (k) of the steady-state activations and inactivations, and the time constant functions of activation or inactivation (τ), were changed one-at-a-time while holding the others at their fixed values. At each alteration, the influence of the parameter change on the model output was assessed by examining the form and duration of the initial burstings in a simulated AP evoked with a 50 s current clamp stimulus. The results were ranked and evaluated. We found that no single intrinsic parameter alteration of the ionic currents could produce a transition from short to long duration bursting APs. However, modification of some parameters could prolong the initial spiking period during an AP; in particular, changes in the maximal conductance ( K1) and the fast inactivation time constant (τ r1) of the uterine delayed rectifying K+ current, I K1 (not to be confused with the time-independent inward rectifying potassium current commonly designated also as I K1 in cardiac cells). Figure 2 shows the result of changing the kinetic properties of I K1 on the uterine AP. Compared to the control case (Figure 2A), the initial bursting duration was slightly increased when K1 was increased (Figure 2B). A similar result was obtained when τ r1 was scaled up but there were no apparent changes when either the activation time constant (τ q) or the slow inactivation time constant (τ r2) of I K1 was changed. Since all these parameters influence the same ionic current, I K1, their combined influences could exert more noticeable changes on the AP form than changing one parameter in isolation. Thus, we examined the effect of I K1 having overall slower dynamics (τ q, 2X; τ r1, 20X; τ r2, 20X) but this only slightly prolonged the initial bursting duration (Figure 2C). However, when an overall slower I K1 was combined with a larger conductance then the bursting continued until the end of the stimulus (Figure 2D). The resultant AP resembled the experimentally recorded long bursting behaviors [1], [2]. This simple parameter analysis provided an insight that, in order to simulate the long-lasting bursting APs, the uterine cell model required larger and slower delayed rectifying K+ currents. To make such changes we needed to either describe a known mechanism for causing such changes in uterine I K1 or replace the uterine I K1 with other molecularly identified K+ currents with such biophysical properties. There is no experimental information to support the former, and so we proceeded with the latter.

Figure 2

Altering the kinetic properties of I K1 in the initial USMC model was found to favor longer bursting-type APs.

A, the standard configuration of the initial USMC model [4] as depicted in Figure 1D. B, increasing the I K1 conductance ( K1 from 0.52 nS pF−1 to 0.64 nS pF−1) prolonged the initial bursting phase a little without enabling repolarization of the AP. C, alternatively, scaling the I K1 time constants (activation: τ, 2X; inactivation: τ 1, 20X and τ 2, 20X) also prolonged the bursting phase and facilitated the membrane voltage repolarization at the end of the stimulus. D, by changing both the conductance and time constants of I K1 together, the bursting persisted until the end of the stimulus and the membrane voltage repolarized on stimulus cessation.

Uterine whole-cell voltage-gated K+ currents have only been roughly categorized based on their inactivation properties and their sensitivity to different broad-spectrum potassium channel blockers. For example, the uterine I K1 has been described as a prominent voltage-gated, delayed rectifying K+ current rather loosely identified by its pharmacological profile and activation/inactivation dynamics [5]. Under voltage-clamp conditions, the activation threshold where I K1 first appears is between −60 mV to −40 mV, it has a half activation of ≈1.2 mV and a half-inactivation of ≈−65 mV. The dynamics of I K1 are slow and the decay of the current can take as much as 10 secs. Also, I K1 is sensitive to 10 mM TEA but insensitive to 5 mM 4-AP. Given these rather broad characterizations, there are contributions to this current that are likely to come from a number of molecularly distinct channels. When considering the possible identity of these channels we noted the following:

(i) In rat uteri at late pregnancy mRNA encoding a very slowly activating K+ channel increased in expression. This K+ current shares similarities to the cardiac delayed rectifying currents I Ks and I Kr [6], [7].

(ii) The delayed rectifying cardiac I Ks consists of the KCNQ voltage-gated pore-forming channel(s) and KCNE subunits and recent evidence suggests that uterine transcripts encoding KCNQ and KCNE subunits are gestationally regulated [8].

(iii) KCNQ activation time constants and conductances are subject to modulation by the pattern of KCNE subunit co-expression [9].

(iv) KCNQ potassium channels can be blocked by XE991 and XE991-sensitive currents in portal venous smooth muscle showed broadly similar characteristics as uterine smooth muscle I K1 with a half-activation of 7 mV and a half-inactivation of −54 mV; the activation threshold is −30 mV [10].

(v) KCNQ currents have pharmacological profiles similar to I K1 as many members of the KCNQ channels are sensitive to TEA and insensitive to 4-AP [11], [12].

(vi) A major contributor to the delayed rectifying cardiac I Kr is the human Ether-à-Go-Go-related (hERG) channel current. Recent evidence has shown that hERG channel subunits are expressed in uterine cells and the currents sensitive to hERG inhibitors, e.g. E4031 or dofetilide, are functionally suppressed at late pregnancy [13], [14].

(vii) The dynamics of hERG currents are comparable to I K1: the activation threshold is −50 mV [15], [16]; a half-inactivation is between -94 mV and −48 mV [15]–[18] and, depending on the recording conditions, the half-activation of hERG varies between −31 mV to 17 mV [15]–[21].

These data, together with the results of the parameter analysis, led us to speculate that KCNQ and hERG currents may form part of the uterine I K1 current. Our next task was therefore to assess how best to incorporate this experimental information into the USMC model. However, there are no experimental measurements specific to uterine smooth muscle cells on the dynamics of KCNQ currents and there is only limited information about hERG-like currents in these cells. Thus, we next sought to (i) establish quantitative descriptions of the biophysical details from clonal data for KCNQ and hERG currents and (ii) modify the USMC model by substituting I K1 with KCNQ and hERG currents.

Formulation of mathematical descriptions of KCNQ1, KCNQ4, KCNQ5 and hERG currents

The KCNQ subunits form the pore of K+ channels. There are five members in the KCNQ family: Q1–Q5. In mouse uterine tissues, the mRNA expressions of all KCNQ subunits were suppressed from early to mid-gestation until late term when KCNQ1, KCNQ4 and KCNQ5 mRNA levels were increased [8]. Based on this information, we decided to formulate descriptions for currents determined by KCNQ1 (I KCNQ1), KCNQ4 (I KCNQ4) and KCNQ5 (I KCNQ5).

There are at least three main variants of hERG protein (hERG1-3), the expressions of which impose different kinetic characteristics on the channel current. The full-length hERG (also known as hERGA or hERG1a) is found in all types of muscles cells. hERG1a is known to form a N-terminus splice variant (known as hERGb or hERG1b), which shows faster activation and deactivation, and smaller tail currents compared to hERG1a [22]. A C-terminus splice variant of hERG (also known as hERGC [19], HERGUSO [20] or erg-sm [16]) has also been reported, including from smooth muscle. However, the descriptions of its characteristics vary from being non-functional channels to fully conducting channels. Greenwood et al. (2009) [13] showed that both hERG1a and hERG1b, but not hERG2 or hERG3, were expressed in mouse myometrium with hERG1a being the most abundant isoform. Therefore, we developed formulations for hERG current (I hERG) based upon published datasets for hERG1a.

Detailed information on the kinetics and steady-state voltage (V)-dependent parameters of I KCNQ or I hERG has been obtained from clonal data from heterologous expression cell systems for each channel subunit: KCNQ1 [10], [23]–[35], KCNQ4 [25], [30]–[32], [34], [36]–[39], KCNQ5 [31], [32], [34],[36],[38]–[42] and hERG [15], [16], [19], [20], [43]. Therefore, mathematical descriptions of the channel biophysical details were developed from these data. Mathematical descriptions of the biophysical characteristics of these currents are given in Appendix S1: for I KCNQ1 see equations 1–13; for I KCNQ4 see equations 14–20; for I KCNQ5 see equations 21–32; for I hERG see equations 33–41.

The equations for I KCNQ1 incorporate two activation gating variables (n Q1f and n Q1s) and two inactivation gating variables (w Q1 and s Q1). Steady-state values for activation and inactivation are shown in Figure 3A-B. The time constants of activation (τ Q1f and τ Q1s) and inactivation (τ Q1 and τ Q1) are depicted in Figure 3C-D and were each obtained by fitting the raw data current tracings from the literature [23]–[35]. τ Q1 is V-dependent and fitted to experimental data [23], [26], [27], [30], [33], [35]. τ Q1 is set as a constant at 50 s. This value was chosen to best reproduce the published time tracings from Jensen et al. (2007) [32]. Simulated traces of I KCNQ1 under voltage-clamp conditions are presented in Figure 3E–F and show dynamic profiles similar to the raw data [32]. The simulated I–V relationship of I KCNQ1 is shown in Figure 3G.

Figure 3

Biophysical characteristics of KCNQ1 current.

Properties of I KCNQ1 were developed using experimental data from human KCNQ1 clones expressed in Xenopus laevis oocytes (black) and Chinese hamster cell (gray) [23]–[35], as well as KCNQ-like currents in mouse smooth muscle cells [10]. A, V-dependent activation steady state (n Q1∞). B, V-dependent steady states for fast inactivation (w Q1∞) and slow inactivation (s Q1∞). C, V-dependent activation time constants (τ Q1f and τ Q1s). D, fast (τ Q1) and slow (τ Q1) inactivation time constants. τ Q1 is V-dependent and fitted to experimental data [23], [26], [27], [30], [33], [35]. τ Q1 is set as a constant at 50 s. This value is chosen to best reproduce the published time tracings from Jensen et al. (2007) [32]. E and F, simulated time tracings of I KCNQ1 with two different voltage-clamp protocols from [32] (insets) and the experimental time tracings are shown for comparison. (Experimental tracing adapted with permission from Jensen et al. (2007); copyright 2007, Biophysical Society.) G, simulated I-V relationship of I KCNQ1 using the same voltage-clamp protocol as in E. All I-V data are normalized to the maximal current value at V = 40 mV.

The equations for I KCNQ4 incorporate an activation (n Q4) and an inactivation (s Q4) gating variable. Steady-state values for activation and inactivation are shown in Figure 4A-B. Blom et al. (2009) [34] reported that KCNQ4 currents were best described with a single time constant for its activation within the voltage range of −40 mV to +40 mV, and two time constants for its deactivation. Thus we used a single time constant for activation (τ Q4) and its voltage-dependency is depicted in Figure 4C. We modeled the current with one inactivation time constant (τ Q4) (Figure 4D). Raw data concerning the dynamics of inactivation of I KCNQ4 are limited and therefore, for simplicity, a function was chosen to fit with the data available and which best reproduced the published raw data time tracings from Jensen et al. (2007) [32]. Simulations of raw data traces of I KCNQ4 under voltage-clamp conditions matched well to the experimental data [32] (Figure 4E–F) and the simulated I–V relationship is shown in Figure 4G.

Figure 4

Biophysical characteristics of KCNQ4 current.

Properties of I KCNQ4 are developed using experimental data from mammalian KCNQ4 clones expressed in Xenopus laevis oocytes (black) and other expression systems (gray) [25], [30]–[32], [34], [36]–[39]. A, V-dependent activation steady state (n Q4∞). B, V-dependent inactivation steady state (s Q4∞). C, V-dependent activation time constant (τ Q4). D, V-dependent inactivation time constant (τ Q4). E and F, simulated time tracings of I KCNQ4 with two different voltage-clamp protocols from Jensen et al. (2007) [32] (insets) and the experimental time tracings are shown for comparison. (Experimental tracing adapted with permission from Jensen et al. (2007); copyright 2007, Biophysical Society.) G, simulated I-V relationship of I KCNQ4 using the same voltage-clamp protocol as in E. All I-V data are normalized to the maximal current value at V = 40 mV.

The equations for I KCNQ5 incorporate two activation gating variables (n Q5f and n Q5s). Steady-state values for the activation are shown in Figure 5A. The dynamics of the activation are described by two time constants (τ Q5f and τ Q5s, Figure 5C). However, it was less straightforward to mathematically describe the steady-state V-dependent inactivation of I KCNQ5. For example, using functions that fitted well with the reported steady-state inactivation data of Jensen et al. (2007) [32], it was not possible to replicate the raw data current tracings from the same publication. Moreover, the hooked tail currents of I KCNQ5 shown in Jensen et al. (2007) [32] suggested that, in order to reproduce this feature, an inactivation faster than the activation may be required. Without sufficient other information concerning the inactivation dynamics of I KCNQ5, we have, therefore, used the fast inactivation of I KCNQ1 and the slow inactivation of I KCNQ4 to represent the fast and slow inactivation conditions for I KCNQ5, i.e., w Q5∞ =  w Q1∞, τw Q5 =  τw Q1, s Q5∞ =  s Q4∞ and τ Q5 = τ Q4. With these adjustments (Figure 5B,D), the simulated I KCNQ5 profiles were satisfactory compared to the published experimental raw data time tracings (Figure 5E–F). The simulated I–V relationship is shown in Figure 5G.

Figure 5

Biophysical characteristics of KCNQ5 current.

Properties of I KCNQ5 are developed using experimental data from cloned mammalian KCNQ5 expressed in Xenopus laevis oocytes (black) and other expression systems (gray) [31], [32], [34], [36], [38]–[42]. A, V-dependent activation steady state (n Q5∞). B, V-dependent steady states for fast inactivation (w Q5∞) and slow inactivation (s Q5∞). C, V-dependent activation time constants (τ Q5f and τ Q5s). τ Q5f is V -dependent but τ Q5s is set as a constant at 1 s. D, V-dependent time constants for fast inactivation (τ Q5) and slow inactivation (τ Q5). E and F, simulated time tracings of I KCNQ5 with two different voltage-clamp protocols (insets) used in Jensen et al. (2007) [32] and the experimental time tracings are shown for comparison. (Experimental tracing adapted with permission from Jensen et al. (2007); copyright 2007, Biophysical Society.) G, simulated I–V relationship of I KCNQ5 using the same voltage-clamp protocol as in E. All I–V data are normalized to the maximal current value at V = 40 mV.

The equations for hERG incorporate two activation gating variables (h and h) and an inactivation gating variable (h). Steady-state values for activation and inactivation are shown in Figure 6A–B. The time constants of activation (τ and τ) and inactivation (τ) are illustrated in Figure 6C–D. The kinetics of hERG currents have rather slow activation and deactivation (msec to sec) profiles but rapid inactivation (msec). As a result, the time tracings of raw data of I hERG under voltage-clamp conditions often show large ‘hook’ tail currents when the command voltage is stepped down from a depolarized level. Simulations of these current tracings are shown in Figure 6E. The corresponding I–V relationships, of either the end-of-pulse current or the peak tail current, are depicted in Figure 6F.

Figure 6

Biophysical characteristics of hERG current.

Properties of I hERG for the myometrial cell model are developed using experimental data from full length hERG clones expressed in different expression systems [15], [16], [19], [20], [43]. A, V-dependent activation steady state (h ∞). B, V -dependent inactivation steady state (h ∞). C, V-dependent activation time constants (τ and τ). D, V-dependent inactivation time constant (τ). E, simulated time tracings of I hERG. Currents were evoked from a holding potential (V h) of −70 mV to various voltage steps (V step) for 2 s then stepped back to V. F, simulated I-V relationships of I hERG using the same voltage-clamp protocol as in E. Both the currents at the end of V step (empty points) and the following maximal tail currents after the end of V step (solid points) were plotted against V step. In both E and F, all data are normalized to the peak tail current value at V = 60 mV.

Substitution of IK1 in the USMC model by individual KCNQ or hERG currents

As rationalized above, in points (i)–(vii), I KCNQ or I hERG or both could be constituent components of the global, rather ill-defined uterine K+ current measured as I K1. We proceeded, therefore, to assess if I K1 in the USMC model could be replaced by I KCNQ or I hERG or a combination of each. As there is little information about the current densities for I KCNQ or I hERG in uterine cells, and the expression of uterine KCNQ and hERG subunits changes during pregnancy, we fitted the values of their maximal conductances so that the simulated maximal whole cell K+ currents matched the experimental value from Wang et al. (1998) [44].

First, we attempted to completely replace I K1 in the USMC model with, in turn, I KCNQ1, I KCNQ4, I KCNQ5 or I hERG. In particular, we were interested in seeing if this could change the AP forms. Each of these maneuvers resulted in the simulation of interesting long-duration AP forms, including a long bursting-type AP form (Figure S1). However, in each case, the corresponding simulations of whole cell K+ current tracings and I–V plots did not match the experimental data (middle and right panels of Figure S1B-E). In native USMC cells [5], [44], and in the original USMC model [4], the majority of myometrial I K1 would be inactivated by a holding potential (V h) of −40 mV due to the negative half inactivation of I K1. So, under voltage-clamped conditions, a larger outward K+ current would be evoked from a V h of −80 mV compared to the current evoked from a V h of −40 mV (Figure S1A middle panel). However, with complete replacement of I K1 by individual I KCNQ1, I KCNQ4, I KCNQ5 or I hERG, the resultant whole cell steady-state K+ current amplitude was the same when evoked from either V h. This was a deviation from the previously matched experimental results. Nonetheless, this exercise showed that I KCNQ or I hERG, could contribute to the development of long duration burstings APs in uterine cells.

Next, we partially replaced I K1 in the cell model with, in turn, I KCNQ1, I KCNQ4, I KCNQ5 or I hERG, while ensuring the simulated whole cell K+ current matched to the experimental data from Wang et al. (1998) [44]. With the standard model configuration, before any changes were implemented, K1 =  0.52 nS pF−1 (Figure S2A). We proceeded with the partial substitution in steps. The level of I K1 was first reduced by changing K1 to 0.24 nS pF−1 (Figure S2B). At this reduced level of I K1, the long-duration AP configuration shows only a plateau-type form with no initial bursting. However, the voltage-clamp experimental results for the whole cell K+ current from different holding potentials can still be reproduced and so this level of I K1 was chosen as a reasonable modification with which to investigate the effects of partially replacing I K1 by I KCNQ or I hERG. Within the constrain of matching the experimental whole cell K+ currents (middle panels) and the I-V relationships (right panels) from Wang et al. (1998) [44], we identified the conductance ranges for I KCNQ or I hERG that favor long-duration bursting APs. With reduced I K1 and incremental additions of I KCNQ1, within a KCNQ1 range of 0.026 – 0.042 nS pF−1, complex AP forms could be simulated whereby the AP plateau was interspersed with periodic bursting (Figure S2B-E). Alternatively, with reduced I K1 and incremental additions of I KCNQ4, within a KCNQ4 range of 0.03216 − 0.04712 nS pF−1 there was an increasing tendency towards bursting behavior persisting throughout the AP albeit with diminishing amplitude (Figure S3C–E). With reduced I K1 and incremental additions of I hERG, within a hERG range of 0.112 – 0.208 nS pF−1, there was an increasing tendency towards improving the bursting behavior of long-duration AP (Figure S5C-E). In contrast, if the reduced I K1 was replaced with incremental additions of I KCNQ5, there was no observable change in the AP form from that of a plateau-type AP (Figure S4C-E, left panel).

The tendency of I KCNQ1, I KCNQ4 or I hERG, to improve the bursting potential of the AP form, suggested that each current individually - I KCNQ1, I KCNQ4 or I hERG - could improve the capabilities of the USMC model in the manner we set out to achieve. For I KCNQ5 the situation appeared different to the other three currents in that after its partial replacement of I K1 the AP forms remained tonic and plateau-like (Figure S4). The effect of incremental inclusion of I KCNQ5 was to reduce the plateau potential and the resting membrane potential. Thus I KCNQ5 does not seem to play a role in bursting APs generation. In addition, there may be species-dependent differences in KCNQ5 expression patterns close to term and/or during labor with human myometrium expressing less KCNQ5 than KCNQ1 or KCNQ4 [8]. This, coupled to the rather negative membrane potentials for activation of the current (Figure 5A), led us to speculate that if I KCNQ5 is to have a role in USMC, it may be for setting the level of resting membrane potential or plateau potential or both rather than contributing to the macroscopic action potential generated K+ current. In our formulation of the original USMC model [4], a background K+ current, I b, was included to represent a collection of minor K+ currents with unknown kinetics that were necessary to hold the membrane potential at physiologically relevant resting levels. The above information suggested that the biophysically detailed I KCNQ5 may ably replace the rather vaguely defined I b. Indeed, when this was done, the cell model was still capable of reproducing all the previous validation results (data not shown). Moreover, with a small increase in I KCNQ5, the model repolarized at the end of the stimulation, while in the original cell model, a small increase of I b did not cause the same behavior (data not shown). Furthermore, replacing I b with I KCNQ5 did not affect the results of partially substituting I K1 by I KCNQ1, I KCNQ4 or I hERG (data not shown). Thus, our first confirmed modification of the initial USMC cell model [4] was to replace I b with I KCNQ5.

Given that uterine smooth muscle cells express components of KCNQ and hERG channels in a gestational dependent manner, we sought to test the proposition that the long duration bursting APs would be facilitated by the incorporation of a combination of I KCNQ1, I KCNQ4 and I hERG in place of reduced I K1 (Figure 7A). Implementation of this scenario - with reduced I K1 ( K1 =  0.24 nS pF−1), a combination of I KCNQ1 ( KCNQ1 =  0.0032 nS pF−1), I KCNQ4 ( KCNQ4 =  0.024 nS pF−1) and I hERG ( hERG =  0.08 nS pF−1) and with I b replaced by I KCNQ5 ( KCNQ5 =  0.016 nS pF−1) - accomplished our objective of reproducing the experimentally reported long-duration bursting AP forms [1]–[3]. Thus by incorporating biophysical details of I KCNQ1, I KCNQ4, I KCNQ5 and I hERG, these modifications corrected a main limitation of the original USMC model. It is likely that this was accomplished by strengthening the repolarization reserve [45] of the uterine cell thereby providing more ways to shape the USMC APs.

Figure 7

USMC model modification with incorporation of KCNQ1, KCNQ4, KCNQ5 and hERG currents.

A, Effect of I K1 being partially replaced by a combination of I KCNQ1 ( KCNQ1 =  0.0032 nS pF−1), I KCNQ4 ( KCNQ4 =  0.024 nS pF−1) and I hERG ( hERG =  0.08 nS pF−1). Note that I b has been replaced with I KCNQ5 ( KCNQ5 =  0.016 nS pF−1). With this configuration, long-duration bursting APs were induced (left panel) while the whole cell K+ currents (middle panel) and the I-V relationships (right panel) matched the experimental data [44] (© Wang et al., 1998). B. Summary of the electrogenic components of the modified USMC model.

Effects of changing hERG current conductance on USMC action potentials

Each of I KCNQ1, I KCNQ4 or I hERG contributed to the ability of the modified USMC model to produce longer bursting APs. Of these, I hERG appeared to exert the most potent influence in this regard (Figure S5). Indeed, in other electrically excitable cells, I hERG contributes to AP spike frequency adaptation [46]. It is noteworthy that hERG-like currents have been reported to be decreased in murine and human myometrial tissues at late pregnancy [13], [14] . Our simulations also showed that progressive reductions in I hERG conductance changed the simulated AP form from that of the long duration bursting AP (Figure 8A) to a more tonic plateau-type AP (Figure 8B-E). This result is similar to the effect of dofetilide on human myometrial APs [14] and provides a validation for the modified USMC model.

Figure 8

Effects of changing hERG current conductance on the modified USMC model AP form.

A, The same conditions and AP form of the modified USMC model as depicted in Figure 7 wherein hERG =  0.08 nS pF−1. B-E, the effects of progressively reducing I hERG conductance. B, hERG =  0.0464 nS pF−1; C, hERG =  0.0224 nS pF−1; D, hERG =  0.0064 nS pF−1; E, hERG =  0.0008 nS pF−1.

Potential effects of changing both KCNQ and hERG current conductances on USMC action potentials

Both KCNQ and hERG currents are present in cardiac muscles and each contributes to stable and timely repolarization [45]. In uterine tissues, the regulations of KCNQ and hERG currents, and/or their molecular components, appear different from each other. Currents sensitive to hERG inhibitors have been reported in uterine tissues in early to mid gestation but these are decreased close to term [13], [14]. In contrast, almost all KCNQ alpha subunits expressions were suppressed during early to mid-gestation and increased at late term [8]. These data suggest KCNQ and hERG currents/components are gestationally regulated in a different manner and raises the question what may the effect be of changing the contributions of each of these particular K+ channel conductances? We utilized simulations of the modified USMC model to examine this.

We began our simulations with no KCNQ currents but with full hERG current (Figure 9A). The result was a plateau-type AP with no initial bursting, qualitatively similar to the plateau-type APs recorded from uterine smooth muscle cells of mid-gestation pregnant rats [1]. We then examined the effects of either increasing I KCNQ or reducing I hERG. With a little I KCNQ, but no change in I hERG, complex oscillations of membrane potentials began to develop around the plateau level (Figure 9B). With more I KCNQ, but still no change in I hERG, a bursting-type AP could be formed (Figure 9C), qualitatively similar to those bursting-type APs recorded from the uterine smooth muscle cells at end-of-term in rats [1]. Note that this is the configuration for our modified USMC model described in previous sections. Further increase of I KCNQ, and still no change in I hERG, suppressed the AP (Figure 9D). (As we do not know the maximum quantities of I KCNQ in uterine smooth muscle cells, it is entirely possible for the I KCNQ to exceed our reference (1.0) values). With I KCNQ staying at the same level, but with reduced I hERG, the bursting-type AP can be formed again (Figure 9E). Further reduction of I hERG suppressed the size of the bursting spikes and changed the AP into a complex form (Figure 9F). These simulations show that although I KCNQ and I hERG may be inversely regulated at late pregnancy, the alteration of conductances for each current potentially yields a fine tuning mechanism for uterine cells to develop and maintain long-duration bursting AP forms that are required for parturition.

Figure 9

Effects of changing KCNQ and hERG current conductances on the modified USMC model AP form.

Simulating the effects of increased expression of KCNQ subunits and reduction of HERG current in uterine smooth muscle cells during late term. The conductance values for I KCNQ and I hERG in the modified USMC model, as in Figure 7, are used as the reference values. Both I KCNQ1 and I KCNQ4 were changed by the same proportions. A, no I KCNQ1 or I KCNQ4; I hERG at 1 fold. B, I KCNQ1 and I KCNQ4 at 0.25 fold; I hERG at 1 fold. C, I KCNQ1, I KCNQ4 and I hERG at 1 fold. D, I KCNQ1 and I KCNQ4 at 1.25 fold; I hERG at 1 fold. E, I KCNQ1 and I KCNQ4 at 1.25 fold; I hERG at 0.6 fold. F, I KCNQ1 and I KCNQ4 at 1.25 fold; I hERG at 0.05 fold.

Conclusions and future considerations

We have identified I KCNQ and I hERG as important components for forming the characteristics long duration bursting-type APs of uterine smooth muscle cells. Also, the above computational approaches have enabled the establishment of a biophysically detailed mathematical model of uterine smooth muscle cell electrical excitability that reproduces the repertoire of experimentally-derived AP forms. In addition, this computational work alerts one to important biological features that would otherwise not be foreseen and can serve as an important tool for hypotheses generations and subsequent experimental testing. A key example is the model prediction that alterations of I KCNQ and I hERG have the potential to radically affect uterine AP forms and, as such, suggests that the physiological regulations of I KCNQ and I hERG may be crucial factors in determining labor onset and maintenance. Of note, there is increasing appreciation of the important roles that computational biology can make to improving our understanding of the complexities of uterine excitability and labor onset [47]–[48], to which the open-access modified USMC model herein looks to contribute by providing the research community with the opportunity to iteratively improve upon the biophysical detail of the components. This will include considerations of intracellular (e.g. metabolic-excitation-contraction coupling characteristics [49] and Ca2+ binding protein affinities) and intercellular (e.g. interstitial K+ accumulation [50]) mechanisms that may act on ionic currents to modulate AP form. To provide tissue-level models that encapsulate the whole process of excitation-contraction coupling one also has to consider spatio-temporal patterns of excitation [48]. Notwithstanding these longer-term research aspirations, our modified model of uterine cell excitability adds to our wider knowledge of the molecular and electrophysiological basis of cellular rhythms [51].

Methods

The original USMC model from [4] describing 14 electrogenic membrane currents, intracellular calcium changes and associated force was used as the basis of this study. Detailed descriptions of the model equations and parameters values are given in [4].

This theoretical study was conducted in two parts. First, we performed an extensive one-at-a-time sensitivity analysis [52] on the USMC model kinetic properties to identify potential target parameters that can prolong or improve the initial bursting in an AP. Then, based on the results, we focused on including more kinetic details of the delay rectifying K+ currents. The USMC model is then modified by including four new voltage-gated K+ currents from the KCNQ and HERG families: I KCNQ1, I KCNQ4, I KCNQ5 and I hERG. The kinetics of these four K+ currents were described by Hodgkin-Huxley type formulations using published clonal data in various expression systems as there are no publish data from myometrial cells or tissues. The details of these currents are described in the Results and Discussion. The definitions of the symbols, the initial conditions and parameters values, the equations of the new potassium currents and the resulting modified USMC model are given in Table S1–3 and Appendix S1. The source code of the modified USMC model written in the C programming language is included in Appendix S2.

Unless otherwise stated, all simulated APs were evoked with a 50 s current clamp stimulus of −0.25 pA pF−1 at t = 1 s starting from the same initial conditions. For each maneuver of incorporating new currents into the cell model, the whole cell K+ currents and the corresponding I-V relationships were also constructed and validated against existing experimental data. Both the whole cell K+ currents and the I-V relationships were simulated by the same one-step voltage-clamp protocol: the membrane voltage was stepped up from a holding potential (V h) of either −40 mV or −80 mV to a step potential (V step) between -40 mV to 70 mV for 250 ms before stepping back to the V h. The model was started at the steady state conditions at V h for each voltage-clamp step. For clarity, and to aid the readers comparison with the corresponding experimental time tracings from Wang et al. (1998) [44], only one simulated time tracing at V step =  0 mV from each of the two V h conditions were shown in our figures and Supplementary figures, and both the simulated and the extracted experimental time tracings were normalized to the maximal value from V h =  −80 mV. The simulated I-V relationships were constructed from the maximum currents at each of the V step normalized to the peak value at V step =  70 mV from V h =  −80 mV.

Simulations were computed with a fixed time step of 0.02 ms, using XPPAUT [53] with the fourth-order Runge-Kutta numerical integration method, in a IBM laptop PC with a Intel(R) Pentium(R) M 1.5 GHz single processor.

Effects of I K1 being completely replaced by either I KCNQ1, I KCNQ4, I KCNQ5 or I hERG. A, left hand panel: the standard configuration of the initial USMC model ( K1 =  0.52 nS pF−1), as depicted in Figure 1C, before I K1 replacement; middle panel: two simulated USMC whole cell K+ currents (solid lines) at a V step of 0 mV, one tracing was from a V of -40 mV (blue) and another from a V of -80 mV (red), and the corresponding experimental time tracings (broken lines) of Wang et al. (1998) [44] (© Wang et al., 1998); right hand panel: the simulated USMC K+ I-V relationships (lines) from both V and the corresponding experimental data (points) of Wang et al., (1998) [44] (© Wang et al., 1998). B-E, the effects of I K1 being replaced by: B, I KCNQ1, C, I KCNQ4, D, I KCNQ5 and E, I hERG.

(TIF)

Click here for additional data file.

USMC model modification with reduced Effects of I K1 being partially replaced by I KCNQ1. A, the standard model configuration ( K1 =  0.52 nS pF−1) before I K1 replacement. Left hand panel: the standard AP simulation as depicted in Figure 1C; middle panel and right hand panel: the whole cell USMC I-V data. B-E, reduced I K1 ( K1 =  0.24 nS pF−1) with different levels of KCNQ1 added: B, KCNQ1 =  0 nS pF−1; C, KCNQ1 =  0.028 nS pF−1; D, KCNQ1 =  0.032 nS pF−1; E, KCNQ1 =  0.036 nS pF−1.

(TIF)

Click here for additional data file.

USMC model modification with reduced Effects of I K1 being partially replaced by I KCNQ4. A, the standard model configuration ( K1 =  0.52 nS pF−1) before I K1 replacement. Left hand panel: the standard AP simulation as depicted in Figure 1C; middle panel and right hand panel: the whole cell USMC I-V data. B-E, reduced I K1 ( K1 =  0.24 nS pF−1) with different levels of KCNQ4 added: B, KCNQ4 =  0 nS pF−1; C, KCNQ4 =  0.03568 nS pF−1; D, KCNQ4 =  0.0392 nS pF−1; E, KCNQ4 =  0.04272 nS pF−1.

(TIF)

Click here for additional data file.

USMC model modification with reduced Effects of I K1 being partially replaced by I KCNQ5. A, the standard model configuration ( K1 =  0.52 nS pF−1) before I K1 replacement. Left hand panel: the standard AP simulation as depicted in Figure 1C; middle panel and right hand panel: the whole cell USMC I-V data. B-E, reduced I K1 ( K1 =  0.24 nS pF−1) with different levels of KCNQ5 added: B, KCNQ5 =  0 nS pF−1; C, KCNQ5 =  0.0192 nS pF−1; D, KCNQ5 =  0.02 nS pF−1; E, KCNQ5 =  0.0208 nS pF−1.

(TIF)

Click here for additional data file.

USMC model modification with reduced Effects of I K1 being partially replaced by I hERG. A, the standard model configuration ( K1 =  0.52 nS pF−1) before I K1 replacement. Left hand panel: the standard AP simulation as depicted in Figure 1C; middle panel and right hand panel: the whole cell USMC I-V data. B-E, reduced I K1 ( K1 =  0.24 nS pF−1) with different levels of hERG added: B, hERG = 0 nS pF−1; C, hERG =  0.112 nS pF−1; D, hERG =  0.144 nS pF−1; E, hERG =  0.176 nS pF−1.

(TIF)

Click here for additional data file.

Definitions of the equation symbols.

(PDF)

Click here for additional data file.

Initial values of the dynamical variables used for the model.

(PDF)

Click here for additional data file.

Constant parameter values used for the model.

(PDF)

Click here for additional data file.

Equations used for the model.

(PDF)

Click here for additional data file.

Source codes in C for the model.

(BZ2)

Click here for additional data file.