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 I K1 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][2][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.

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.
Parameters such as the maximal conductances ( g), 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 (t), 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 ( g K1 ) and the fast inactivation time constant (t 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. 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 21 . C, short plateau-type AP evoked by a 2 s Ist of 20.5 pA pF 21 . D, an attempt to simulate a long duration AP with a 50 s Ist of 20.25 pA pF 21 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.
Compared to the control case (Figure 2A), the initial bursting duration was slightly increased when g K1 was increased ( Figure 2B). A similar result was obtained when t r1 was scaled up but there were no apparent changes when either the activation time constant (t q ) or the slow inactivation time constant (t 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 (t q , 2X; t r1 , 20X; t 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 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 ( g K1 from 0.52 nS pF 21 to 0.64 nS pF 21 ) prolonged the initial bursting phase a little without enabling repolarization of the AP. C, alternatively, scaling the I K1 time constants (activation: t q , 2X; inactivation: t r1 , 20X and t r2 , 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. 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.
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 260 mV to 240 mV, it has a half activation of <1.2 mV and a half-inactivation of <265 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 poreforming 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 254 mV; the activation threshold is 230 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].
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], HERG USO [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.
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 240 mV to +40 mV, and two time constants for its deactivation. Thus we used a single time constant for activation (t nQ4 ) and its voltage-dependency is depicted in Figure 4C. We modeled the current with one inactivation time constant (t sQ4 ) ( 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 voltageclamp conditions matched well to the experimental data [32] ( Figure 4E-F) and the simulated I-V relationship is shown in Figure 4G.
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 (t nQ5f and t nQ5 s, 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' 5 w Q1' , tw Q5 5 tw Q1 , s Q5' 5 s Q4' and t sQ5 5t sQ4 . 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.
The equations for hERG incorporate two activation gating variables (h n1 and h n2 ) and an inactivation gating variable (h s ). Steady-state values for activation and inactivation are shown in Figure 6A-B. The time constants of activation (t hn1 and t hn2 ) and inactivation (t hs ) 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.   [31,32,34,36,[38][39][40][41][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 (t nQ5f and t nQ5s ). t nQ5f is Vdependent but t nQ5s is set as a constant at 1 s. D, V-dependent time constants for fast inactivation (t wQ5 ) and slow inactivation (t sQ5 ).

Substitution of I K1 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 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 n' ). B, V -dependent inactivation steady state (h s' ). C, V-dependent activation time constants (t hn1 and t hn2 ). D, V-dependent inactivation time constant (t hs ). E, simulated time tracings of I hERG . Currents were evoked from a holding potential (V h ) of 270 mV to various voltage steps (V step ) for 2 s then stepped back to V h . 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 V560 mV.  [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 240 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 280 mV compared to the current evoked from a V h of 240 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, g K1 5 0.52 nS pF 21 ( Figure  S2A). We proceeded with the partial substitution in steps. The level of I K1 was first reduced by changing g K1 to 0.24 nS pF 21 ( 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 g KCNQ1 range of 0.026 -0.042 nS pF 21 , 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 g KCNQ4 range of 0.03216 2 0.04712 nS pF 21 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 g hERG range of 0.112 -0.208 nS pF 21 , 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 hERGcould 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 plateaulike ( 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 speciesdependent 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 ( g K1 5 0.24 nS pF 21 ), a combination of I KCNQ1 ( g KCNQ1 5 0.0032 nS pF 21 ), I KCNQ4 ( g KCNQ4 5 0.024 nS pF 21 ) and I hERG ( g hERG 5 0.08 nS pF 21 ) and with I b replaced by I KCNQ5 ( g KCNQ5 5 0.016 nS pF 21 ) -accomplished our objective of reproducing the experimentally reported long-duration bursting AP forms [1][2][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.

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.

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, hERG/KCNQ Modulation of Uterine Electrical Activity 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 hERG/KCNQ Modulation of Uterine Electrical Activity 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.

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 Ca 2+ 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 excitationcontraction coupling one also has to consider spatio-temporal patterns of excitation hERG/KCNQ Modulation of Uterine Electrical Activity [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 20.25 pA pF 21 at t51 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 240 mV or 280 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 5 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 5 280 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 5 70 mV from V h 5 280 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. Figure S1. I K1 cannot be completely replaced by KCNQ1, KCNQ4, KCNQ5 or hERG current. 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 ( g K1 5 0.52 nS pF 21 ), 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 h of -40 mV (blue) and another from a V h of -80 mV (red), and the corresponding experimental time tracings (broken lines) of Wang et al. (1998) [44]