For these simulation studies we used our previously published USMC model (Tong et al., 2011). A schematic of the ionic currents contributing to the model is shown in Figure 1A. The model source code, parameter values and the full description are provided in Tong et al. (2011). All the USMC AP simulations were started at the same resting state. The numerical values of all the dynamical variables at this state—the initial conditions—are provided in the Supplementary Materials.

There are numerous computational models of cardiac cells and many are listed in the CellML model repository (http://models.cellml.org/). They represent different parts of the heart from the central and peripheral sinoatrial nodal cells (Zhang et al., 2000; Garny et al., 2003), the atrial cardiomyocytes (Aslanidi et al., 2009), the atrial-ventricular node and His-bundle (Inada et al., 2009), Purkinje fiber cells (Stewart et al., 2009; Li and Rudy, 2011) and ventricular cardiomyocytes (Faber and Rudy, 2000). These cell types can be grouped, rather roughly, into those of a contractile or conducting function. In our preliminary simulations using these cardiac models, we selected the ventricular and sinoatrial nodal cell models that showed the greatest propensity to alter AP form following a reduction of Ca current thus forming the most sensitive situation for comparing the same maneuvre to that on the USMC model.

For the sinoatrial nodal cell, we used the Garny et al. (2003) rabbit sinoatrial nodal cell model (SAN) and its 0D-capable version configurations for the SAN central cell. A schematic of the ionic currents contributing to the model is shown in Figure 1B. The model equations and parameter values and are listed in the CellML model repository. The initial conditions for the SAN model are provided in the Supplementary Materials.

For the ventricular cardiomyocyte cell, we used as a platform the Faber and Rudy (2000) guinea pig cell model (LRd00) and its M cell configurations. A schematic of the ionic current components of the model is shown in Figure 1C. The model equations, parameter values and initial conditions are listed in the CellML model repository. A copy of the model source code can also be downloaded from http://rudylab.wustl.edu/research/cell/methodology/cellmodels/LRd/code.htm. However, as detailed below in the Results, we had to modify this model in order to incorporate a more robust description of I_CaT. The LRd00 modified model was paced for 20 min at 2 Hz to allow it to reach a stable state with its M cell configurations. Then, the values of all the dynamical variables in this stable state were saved: these stable state conditions of the modified LRd00 model are provided in the Supplementary Materials.

The effects of reduced I_CaL and I_CaT on the uterine and cardiac APs were assessed and quantified. The I_CaL and I_CaT were reduced by multiplying their maximal conductances with a scaling factor between 0 and 1 and we assumed the same proportional reduction would occur in all three cell types. As the APs between the three unique cell types have different forms and characteristics, it is difficult to use a common measure to assess their AP behaviors. Therefore, we chose to assess a characteristic that best reflects the main function of each of the cell types: the bursting activities in the uterine cell, the pacemaking ability of the SAN cell and the action potential duration (APD) of the ventricular cell. For the USMC model, after adjusting the levels of I_CaL and I_CaT, an AP was induced by a 5 s current clamp at −0.5 pA pF⁻¹ and the maximal peak membrane voltage (V_peak>) reached after the first initial spike was measured. This assessment gives an indication of the presence or absence of bursting spikes in a USMC AP. For the SAN model, after adjusting the levels of I_CaL and I_CaT, the model was simulated for 10 s and the pacemaking frequency was measured from the APs in the last 2 s of simulations. For the LRd00 model, an AP was induced by a 0.5 ms current pulse at −90 pA pF⁻¹ and the APD at 90% repolarization level (APD₉₀) was assessed.

All simulations were computed using XPPAUT Version 6 (Ermentrout, 2002) in a Dell Optiplex PC with an Intel(R) Core(TM) i7-3770 CPU @ 3.40 GHz. For the USMC model, all simulations were computed with a fixed time step of 0.02 ms and the Euler method. For the SAN model, all simulations were computed with a fixed time step of 0.1 ms and the fourth-order Runge–Kutta method. For the LRd00 model, all simulations were computed with a fixed time step of 0.002 ms and the fourth-order Runge–Kutta method.

We first verified each source code of the three individual cell models with their respective published results. Both the USMC and the SAN models were validated in this regard. Surprisingly, we found that the LRd00 model could not be validated so in respect of I_CaT: the simulated current tracings (Figure 2A), and the corresponding current-voltage (I-V) relationship (Figure 2B), of the LRd00 model did not match the experimental results that they were rated against (Balke et al., 1992). The simulations produced current responses at membrane voltages that were much farther to the right (positive) than observed in the experimental data or, indeed, to that anticipated from other I_CaT models. When we traced back the original formulation of this I_CaT to Zeng et al. (1995), and the experimental values from which its kinetics details were based to Droogmans and Nilius (1989), we did not find any typographical errors in the model equations or the experimental values to explain the discrepancy. As an aim of our study was to investigate the effects of altering I_CaT, it was essential that the models incorporated the correct formulations. We were left, therefore, with no option but to attempt to modify the I_CaT in the LRd00 model such that it did reflect more closely the experimental data. For those interested in the details of this process, our attempts to modify the model were undertaken as follows: we had first tried, but failed, to produce a reasonable I_CaT model by either adjusting the LRd00 I_CaT equations or reformulating new equations using data from Droogmans and Nilius (1989). We next tried substituting the LRd00 I_CaT with a validated I_CaT model from Li and Rudy (2011), which was developed for canine Purkinje fiber cells. The simulated current tracings with this canine I_CaT resemble the I_CaT current tracings in Balke et al. (1992) but the kinetics were too fast and the peak of the I-V relationship was still too positive compared to the experimental data from Balke et al. (1992). Using the I_CaT from Li and Rudy (2011) as our template, we corrected the differences between the simulations and the experimental data with the modifications that resulted in Equations (1–8) for the new ventricular I_CaT that are listed in the Supplementary Materials. As discussed by Droogmans and Nilius (1989), we also found that the kinetics of the ventricular I_CaT currents were best described by a gating product of b²g where b and g are gating variables for activation and inactivation (Equation 1 in the Supplementary Materials). The I-V data in Balke et al. (1992) was matched with a half-activation at −50 mV in the activation steady-state function (Equation 3 in the Supplementary Materials). Although this half-activation value deviated from the reported value for I_CaT (Droogmans and Nilius, 1989), it is similar to those reported values from clonal Cav3.1 expression data (Serrano et al., 1999; Hering et al., 2003). As the kinetics of the canine I_CaT model were too fast compared to the ventricular cell experimental data within the LRd00 model (i.e., the data from Balke et al., 1992), we slowed down the activation and inactivation by scaling the time constants (Equations 5, 6 in the Supplementary Materials). The resultant new I_CaT model satisfactorily described the ventricular I_CaT data (Balke et al., 1992) and we replaced the original LRd00 I_CaT with this new I_CaT model. The ratio between peak I_CaT and I_CaL in guinea pig ventricular cells are reported to be between 0.13 and 0.19 (Balke et al., 1992; Masumiya et al., 1998; Zorn-Pauly et al., 2004). Therefore, we chose a value for the maximal conductance of our new I_CaT model (g_CaT) so that the ratio of the maximum I_CaT: I_CaL in the modified LRd00 cell model was 0.15.

USMC APs can take a variety of forms of variable durations. Bursting type USMC AP, wherein an initial spike of an AP followed by rapid and repetitive fluctuations in membrane potential that persist for the whole AP duration, are often observed. These may assist in the maintenance of long-duration contractions of many tens of seconds that are a feature of the uterus during labor. Therefore, a promising feature of a tocolytic would be one that could dampen, or delay, the bursting in the USMC APs with minimum effects on the APs of the cardiac SAN and the ventricular cells.

Figure 3A shows the simulated APs from the three cell models under control conditions before any maneuvre. The USMC model showed a resting membrane potential (RMP) of −55.5 mV before stimulation. During the evoked USMC AP, burstings occurred throughout the AP with a frequency between 2 and 2.54 Hz and the amplitudes of the repetitive spikes were around 40 mV. The V_peak> after the initial spike was −7 mV. The SAN cell is autorhythmic where periodic APs occur without external stimulation. The frequency and amplitude of the pacemaking activities were 3.15 Hz and 75 mV (between −56 and 19 mV), respectively. These values were the same throughout the whole 10 s simulation. The LRd00 ventricular cell showed a comparatively hyperpolarized (negative) RMP at −87 mV. A single AP was evoked in respond to a stimulus and the ventricular APD₉₀ was 165 ms. These quantified characteristics are the reference values before any maneuvre.

I_CaL is a major depolarization current in USMCs and thus a logical target for tocolytics. We examined the effect of reducing I_CaL by 10% (Figure 3B). In the USMC model, this change dampened the initial repetitive spikes of the USMC AP and reduced the duration of the burstings by around 60%. The amplitudes of the remaining, late-onset, bursting spikes were reduced by 30% to about 28 mV and the V_peak> of the USMC AP was lowered to −14.8 mV. A less severe effect of 10% less I_CaL was evident in the APs from the SAN and the LRd00 models. The amplitude of the SAN pacemaking APs was modestly reduced by 12% and the frequency remained unchanged. The LRd00 APD₉₀ was slightly reduced by 5%. When I_CaL was reduced by 20% (Figure 3C), the burstings of the USMC AP were suppressed completely. However, the SAN cell also stopped pacemaking. The LRd00 APD₉₀ was slightly reduced by ~10% at 150 ms. These data indicated that although reducing some I_CaL suppressed uterine bursting APs, this maneuvre also affected both the cardiac pacemaker APs and the ventricular AP.

I_CaT may be one of the SAN pacemaking currents in the heart and may also be involved in the generation of APs in USMCs. The effects of our maneuvres to reduce I_CaT in each of the three cell types are shown in Figure 4. The control cases of the three models are displayed in Figure 4A. Reduction of I_CaT by 40% (Figure 4B) delayed the onset of the burstings in the USMC AP. However, both the V_peak> and the amplitude of the later-onset burstings remained the same as the control. With 40% less I_CaT, the pacemaking by the SAN cell slowed down to about 2.94 Hz with the peak potential at a slightly lower level of 15 mV. Reduction of I_CaT by 40% did not affect the LRd00 APD₉₀. When more I_CaT was reduced, by 80% (Figure 4C), the onset of the USMC burstings was further delayed with only one spike appearing at the end of the AP. However, the amplitude of this spike was similar in size to those of the control. Also, the RMP level became more hyperpolarized at ~−57 mV. With this large reduction of I_CaT (by 80%), the SAN pacemaking APs slowed down to 2.63–2.7 Hz. The pacemaking potential in the SAN model remained at the same level at −56 mV but the peak potential was lowered to 12 mV. At this level, the SAN model also showed a longer (~0.5 s) transient during simulation. With 80% less I_CaT, the LRd00 APD₉₀ remained unchanged. Without I_CaT (Figure 4D), the RMP level of the USMC became more hyperpolarized at ~−58 mV and the bursting of the USMC AP stopped. The SAN pacemaking also ceased without I_CaT but the LRd00 AP was unaffected. These data indicated that reduction of I_CaT suppressed the uterine bursting AP with no impact on the ventricular AP. However, depending on the magnitude of I_CaT reduction, there was also an effect on the SAN APs.

The above data suggested that there was some promise in reducing either I_CaL or I_CaT and effecting an inhibition of USMC APs with little deleterious action on ventricular APs. The manouevres that reduced I_CaL or I_CaT and had mild effects on SAN AP may be tolerable. However, reductions in either current that effected a larger inhibition of uterine APs, as one may desire the action of a tocolytic in our experimental setting to be, also had dramatic actions on SAN AP form. The outcomes of reducing I_CaL or I_CaT in the USMC models were slightly different: reducing I_CaT altered the RMP and the onset of spike bursting whereas changes in I_CaL had a predominant action on spike amplitude. This begat the question, could a combined reduction of I_CaL and I_CaT exert the desired action on USMC APs with little impact on cardiac cell APs? An example is illustrated in Figure 5. The control cases of the three models are shown in Figure 5A and the effects of reducing only I_CaL by 10%, or only I_CaT by 40%, are illustrated again in Figures 5B,C, respectively. This enables comparison to the effect arising when both I_CaL and I_CaT were reduced, by 10 and 40%, respectively, (Figure 5D). In this case, the burstings of the USMC AP were completely suppressed (Figure 5D). With this combination, the pacemaking frequency of the SAN APs remained the same at 3.125 Hz and the amplitude was reduced by 20% (from −52 to 7.7 mV). The LRd00 APD₉₀ was slightly reduced by ~4% to 158 ms. Compared to the cases when only either I_CaL or I_CaT were reduced alone to completely suppress the burstings of the USMC AP (Figures 3C, 4D), this combined approach clearly performed better while preserving many of the properties of the cardiac APs.

We plotted two-parameter maps of percentage current inhibitions—of I_CaL vs. I_CaT—in each cell model to examine the effective ranges of different combinations of I_CaL and I_CaT reductions on the properties of uterine and cardiac APs (Figure 6). After the initial spike of a USMC AP, the peak of the bursting spikes typically reach values between −20 and 0 mV whereas without bursting, the plateau AP voltage usually stayed below −30 mV. Therefore, we assessed the presence or absence of burstings in a USMC AP from the level of V_peak> when the proportions of I_CaL and I_CaT were changed and the results of the parameter plots are color-coded in Figure 6. In Figure 6A, the white region indicated the combinations of I_CaL and I_CaT that would produce a USMC AP with some burstings and the shaded region indicated the combinations of currents with which the USMC cell would not produce such activity. The USMC AP example shown previously in Figure 5D, with 10% reduction of I_CaL and 40% reduction of I_CaT, lay just within the non-bursting domain (as indicated by the green star). For the SAN cell, the pacemaker frequency was monitored as one changed the proportions of I_CaL and I_CaT and the results are shown in Figure 6B. The colored domain indicated pacemaking occurrence and the white region suggested the parameter combinations would result in no pacemaking. The frequencies of the pacemaking APs were similar throughout the pacemaking zone and the transition between the two regions was steep (not shown). The SAN example shown in Figure 5D with 10% reduction of I_CaL and 40% reduction of I_CaT was within the pacemaking region. For the LRd00 ventricular cell, the APD₉₀ was monitored as the proportions of I_CaL and I_CaT were changed and the results are shown in Figure 6C. Changes in the APD₉₀ occurred only with changes in I_CaL, indicating that I_CaT did not influence the APD of the ventricular cell model. From overlapping these parameter maps, we derived the boundaries of combinations of I_CaL and I_CaT for attaining specific sets of AP forms in the uterine and cardiac cells as shown in Figure 6D. The colored lines separate different properties of the uterine and cardiac cells obtained from the simulations in Figures 6A–C. The blue line traces the level of V_peak> at −30 mV in the parameter space of I_CaL and I_CaT of the USMC cell model and it separates the bursting and quiescent conditions. The red line separates the pacemaking conditions with I_CaL and I_CaT in the SAN cell model from the non-pacemaking. We allowed for 10% variations in the APD₉₀ of the LRd00 ventricular cells and this threshold at 150 ms was indicated by the gray line. These three criteria were superimposed on the same parameter map of I_CaL and I_CaT and revealed an overlapping area that satisfied all three conditions. In other words, any combination of reduced I_CaL and I_CaT within this area can produce our desired outcome for tocolytic treatment, namely, effective suppression of USMC APs with minimal effects on cardiac APs.

The above experiments have considered the outcomes on AP form when imposing the same proportional reduction of I_CaL and I_CaT in all three cell types. However, another interesting possibility would be the same absolute quantities of currents being reduced in all three cell types. This consideration arises from examining the range of current densities of I_CaL and I_CaT reported from the different cell types as noted in Table 1. Cells of the cardiac conduction system tend to have much more I_CaT than contractile cardiomyocytes and, when comparing between uterine and cardiac cell types, the total I_CaT found in USMCs amounts to around 42% of that reported from guinea pig SAN. We re-performed our simulations to test the effects of reducing the same quantities of currents in each of the three cell type models. This resulted in a transformation of the XY axes of the parameter maps in Figure 6 to current densities and the result is shown in Figure 7. Now, the simulations reveal a larger domain within which one may reduce I_CaL and I_CaT to produce effective inhibition of USMC model APs with little detriment to cardiac cell model APs. In particular, it is possible to completely block the USMC bursting by large reduction of the USMC I_CaT without affecting the SAN pacemaking. Therefore, our first experimental protocol of examining the effects of proportional current reductions likely represents the most stringent case.

Our work brings to the fore the notion that, instead of searching for a single tocolytic compound with high specificity, there may be merit in considering a cocktail of more than one compound. The present recommendations in the UK are to avoid combinations of tocolytic drugs for fear of increasing the risk of side effects (RCOG, 2011). However, this caution arises from a lack of depth to the background research data and the paucity of currently available tocolytic options. The progression to labor has been described as a modular accumulation of (patho)physiological systems—or MAPS—and it is such modularity that is a major challenge to overcome in seeking to prevent or inhibit labor (Mitchell and Taggart, 2009). The suspicion is that once module A of the labor process is inhibited, there is module B, or module C or more that eventually facilitates a similar function and outcome of labor and delivery. If this is so, then the possibilities offered by combination drug strategies must continue to be explored.

Our assessment only considered effects at the cell level. However, to fully evaluate the actions of any tocolytic compounds, we need to also consider their effects at the tissue and organ levels. Currently there is no validated computational organ model for the uterus that would serve well this purpose although it is recognized by many that efforts toward this are required (Aslanidi et al., 2011a; Sharp et al., 2013). We are under no illusion that there is a need for much more “wet” experimentation—particularly in USMCs—to furnish computational model improvement, and to test with increasing rigor and clarity, important hypotheses of relevance to tocolysis.

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.