Introduction

Disorders in uterine excitation or contractility can lead to a range of complications for the mother and child, including preterm birth, ineffective and long labour, and post-partum haemorrhage. Premature births are associated with an increased chance of morbidity and mortality for the child and are the leading cause of death in children under 5 worldwide1. If they survive beyond five years these children tend to be small, lightweight and possibly suffering from mild to severe disabilities2. To reduce the risk of preterm birth, a better understanding of the mechanism(s) underlying the initiation and control of electrical and mechanical activities of the uterus is needed.

The electrophysiology of the human myometrium is complex, both in terms of how the cell and tissue properties change during gestation, and how electrical activity and mechanical contractions are initiated and maintained. In the past several decades, extensive experimental studies3,4,5,6,7,8,9,10 have been conducted to investigate in detail the cellular and tissue electrophysiology of the uterus. It has been shown that uterine mechanical contractions result from the integrated electrophysiology, biochemistry and biomechanics of uterine smooth muscle cells within the myometrial tissue, and their synchronization, with the mechanically passive supporting tissue and its architecture. Pace-making sites, which may play an important role in initiating uterine electrical activity triggering mechanical contraction, have been identified and mapped in the pregnant rat uterus3. However, the exact mechanism(s) underlying the initiation and control of uterine electrical and mechanical activities remains incompletely understood.

Biophysically detailed computational models of uterine cells provide an alternative approach to experimental studies to investigate possible mechanism(s) underlying the genesis of uterine electrical and mechanical activities. Over the past decade, several cell models for uterine smooth muscle cell electrophysiology have been developed. These include the historical Bursztyn et al.11 and Rihana et al.12 models that laid the foundation for current myometrial cell models. In 2011, Tong et al.13 developed the first comprehensive cell model for rat uterine cells that coupled cellular electrophysiology to the intracellular Ca2+ handling and the generation of active force of myofilament. The model was later updated to include some more potassium channels in 201414. In 2016, Atia et al.15 developed the most comprehensive human myometrial electrophysiology model available, based on combining transcriptomic and biophysical data.

The aim of this study was to further develop and update a computational model for simulating the membrane potential and currents, intracellular calcium dynamics and mechanical activities of rat isolated uterine myocytes based on the Tong et al. 2011 model13 of rat uterine myocytes. First, the model was modified to incorporate some newly available experimental data on the kinetics of some membrane ion channels that underlie the membrane potential. This included reformulated equations for the L-type calcium current (ICaL), a voltage-activated potassium current (IK2), the calcium-activated potassium current (IK(Ca)) and the calcium-activated chloride current (ICl(Ca)). Secondly, the model was modified to incorporate more detailed descriptions of the intracellular calcium handling and mechanical dynamics of the uterine cell, including a sarcoplasmic reticulum (SR). This enables the model to be related to cell and tissue recordings obtained by optical imaging of membrane potentials and intracellular calcium16,17. Finally, the developed model was validated by its ability to reproduce the functional impacts of hormone and drug actions, such as those of oestradiol, oxytocin and nifedipine on the electrical and mechanical behaviours of uterine cells. Table 1 summarises and compares the elements included in models previous to the one presented in this study (referred to as Testrow et al. 2018), highlighting the progress of cellular model development, especially in the areas of intracellular Ca2+ handling and cellular biomechanics.

Table 1 A summary of mammalian uterine smooth muscle cell models.

The paper is structured in the following way. In section 2, a brief description of the model development is presented (additional information can be found in the supplementary material). In section 3, simulation results are presented to demonstrate capability of the model in reproducing experimentally observed uterine cellular electrical and mechanical activities, including the spike train action potentials and staircase calcium transients associated with uterine cells. In addition, we demonstrated the model’s capability in reproducing the effects of some drugs and hormones on cellular electrical and mechanical activities. All of these simulation results qualitatively and quantitatively match to experimental data, validating the cell model development. With the validated cell model, we theoretically investigated the role of changing Na+ and Ca2+ ion channel current density in modulating cellular action potentials, the intracellular Ca2+ handling and the genesis of active force, exploring effects of possible ion channel remodelling during the gestational period on uterine cellular electrical and mechanical activities. In Section 4, we discuss the major contributions of this study.

Model Development

This section presents a brief overview of the model development; full details of the model components are provided in the Supplementary Material.

To model the physiology of the single uterine cell, two distinct subsystems (electrochemical and chemomechanical) were coupled as shown in Fig. 1. The electrochemical subsystem is divided into the membrane (electrophysiology), intracellular ion handling and sarcoplasmic reticulum (SR) components. The electrochemical subsystem communicates with the chemo-mechanical subsystem using intracellular Ca2+ as a mediator. The mechanical force generation component handles the Ca2+ activation of contractile proteins and their regulation, and the subsequent production of tension.

Figure 1
figure 1

Schematic illustration of the myometrial cell model49. Left: The electrochemical and chemomechanical subsystems are coupled via the intracellular calcium concentration. The arrows indicate the direction of influence of the coupled systems. Right: Model components of 7 inward membrane channels (red arrows), 5 outward membrane channels (blue arrows), a calcium membrane pump, a Na-K-Cl symporter, 2 membrane ion exchangers, a SERCA pump (Iup), an SR transfer current (Itr) and an SR release current (Irel) controlled by the ryanodine receptors (RyR) were shown. The colours in the flow diagram match their respective systems in the cell diagram.

The cell membrane contains ion channels, pumps and exchangers that permit the passage of charged ions (Fig. 1). These ionic channels are voltage-gated, but some are dependent on ionic concentrations. Typically a channel has activation gates and inactivation gates; each is described by a set of mathematical equations for describing steady states and time constants of the channel (Supplement A).

The Tong et al. model was used as the basal model, which was modified to incorporate detailed models of Ca2+ release from the sarcoplasmic reticulum (SR) and four-component force models, as developed by Yang et al.18, and the kinetic cross-bridge model developed by Hai and Murphy19 (details are provided in Supplement A).

Cellular Electromechanical Activity

The myometrial cell is either at rest or exhibiting electrical and contractile activity. The rest state represents the quiescence exhibited by the cell between periods of activity; it is prevalent during most of gestation20. The active state represents the bursting activity that triggers contraction. Myocytes in the rat uterus are organised in two layers – the inner circular and outer longitudinal layer. Myometrial cells from the circular layer produce plateau-like action potentials (AP), and those from the longitudinal layer are spike-like21.

The action potential of the uterine single cell can also take various forms depending on the stage of pregnancy. Typically they are plateaux, spikes or a combination of the two. A successful labour is a sequence of contractions, and during each contraction there are rapid bursts of consecutive spikes that lead to the accumulation of intracellular Ca2+. Bengtsson et al.22 demonstrated AP can have variant morphologies in late pregnant rat myometrial cells from the circular layer. They begin as single spikes with plateaux. Over time the plateau potential begins to oscillate with increasing frequency and magnitude until it builds up to a rapid burst of spikes. Spike bursts increase in frequency and contractions increase in strength and duration until parturition23.

Gestational changes take place over the course of pregnancy. In non-pregnant rats uterine myocytes have a volume of ~2.5pL that increases 8-fold during pregnancy9. Spontaneous action potentials, typically last tens of seconds with a resting potential of ~−55 mV and a peak spike potential of ~10 mV24. During gestation the peak sodium current approximately doubles, while the peak calcium current halves25.

General Model Equations

The total membrane ionic current density (Itot) is a sum of the currents from the various ion channels, pumps and exchangers. The change in membrane voltage is given by,

$$\begin{array}{rcl}\frac{dV}{dt} & = & -\frac{{I}_{tot}}{{C}_{m}}\\ {I}_{tot} & = & {I}_{h}+{I}_{CaL}+{I}_{CaT}+{I}_{Na}+{I}_{K1}+{I}_{K2}+{I}_{KA}+{I}_{KCa}+{I}_{NSCC}+{I}_{Cl}\\ & & +\,{I}_{Capump}+{I}_{NaK}+0.5{I}_{NaCa}+{I}_{sus}+{I}_{soc}+{I}_{Kleak}+{I}_{NaKCl}\end{array}$$
(1)

where Ih is the hyperpolarisation-activated current, ICaL the L-type calcium current, ICaT the T-type calcium current, INa the fast sodium current, IK1 and Ik2 the voltage-dependent potassium currents, IKA the A-type transient potassium current, IKCa the calcium-activated potassium current, INSCC the non-selective cation current, ICl the calcium-activated chloride current, ICapump the calcium pump current, INaK the sodium/potassium exchanger current, INaCa the sodium/calcium exchanger current, Isus the sustained background current, Isoc the store-operated non-selective cation current, IKleak the background potassium leakage current and INaKCl the sodium-potassium-chloride co-transport current. Current densities are given in units of pA pF−1 and specific membrane capacitance (Cm) has a value of 1.4 µF cm−29.

Details of the model equation development for each of the ion channel currents in Equation (1) are provided in the Supplementary Materials. The data for each ion channel can be found in Supplement A, the equations and parameters in Supplement B and the initial conditions in Supplement C. In simulations, Equation (1) is numerically solved using a forward Euler scheme with a time-step of 0.01 ms that provides stable solutions.

Results

The developed model for the uterine myometrial cell is able to reproduce the typical electrical behaviours of the myometrial cell, including a stable resting potential in the absence of an external stimulus and bursting membrane potentials, which may be periodic and perhaps quasiperiodic and even chaotic in response to external stimuli. Figure 2 shows the simulated bursting action potentials in response to an external depolarising stimulus pulse, with a duration of 5 s and amplitude of −0.4 pA pF−1, which is <12% of the maximal inward ionic current density. It was shown that in response to such an external stimulus, a sequence of bursting action potentials were evoked (Fig. 2A) from a stable equilibrium resting potential, resulting in a stair-case increase of the intracellular Ca2+ concentration (Fig. 2B). Such bursting action potentials were generated by the integral action of opening and closing of many ion channels, some of which were shown in Fig. 2C. The simulated bursting action potentials well matched to those recorded experimentally. The contractile force of the membrane in response to the action potentials is shown in Fig. 2D. The simulated bursting action potentials, stair-case increase of the intracellular Ca2+ transient and the resultant cellular active force matched to examples of experimental data as shown by the insets in Fig. 2, which show the typical morphology of membrane potentials, calcium transients and force profiles in uterine cells from late-pregnant rat.

Figure 2
figure 2

Bursting electrical activity of the model and comparison with experimental data. (A) Simulated bursting action A Inset: experimental data of bursting AP from Bengtsson et al.22. (B) Simulated staircase-like calcium transient generated by bursting APs. B Inset (centre): the stimulus applied during this simulation, −0.4 pA pF−1 for a duration of 5 s. B Inset (right): experimental data from Shmygol et al.28. (C) Computed results of four principal currents: L-type calcium, sodium and voltage gated potassium currents (IK1 and IK2) during the time course of bursting APs. (D) Simulated total contractile force of the cell model. D Inset: experimental data from Zafrah and Alotaibi32.

The effects of AP morphology on [Ca2+]i

Experimentally variant morphologies of the AP have been recorded from myometrial cells22,26. It has also been shown that AP morphology affects the intra-cellular Ca2+ transient8,27 and therefore the resultant force. To test the capability of the model in reproducing the effects of AP morphology on the intracellular [Ca2+]i, in simulations, we implemented variant experimentally measured membrane potentials recorded from freshly isolated 21 day pregnant rat myometrial cells (digitized at 2.5 kHz) as clamp commands to drive the cell model. Results are shown in Fig. 3, in which two distinct forms of action potentials were used: a single spike followed by either a smooth plateau or a train of spikes superimposed on the plateau. In responding to the AP clamp, the model produced variant forms of Ca2+ transient (red lines).

Figure 3
figure 3

Simulated calcium transients produced by voltage clamps. The model voltage was clamped by using the voltage profiles obtained from experimental current clamp recordings, rather than constant values. The experimental data was recorded at 35 °C using the fluorescent indicator indo-1 in patch-clamped, single uterine myocytes from pregnant rats. A comprehensive methodology can be found in Shmygol et al.28. The parameters of each experimental current clamp are described in (AD). (A) A single spike followed by a long plateau, using a holding current of −30 pA and a pulse of 60 pA for 5 s produced a short-lived, low maximum linear rise in [Ca2+]i. (B) A spontaneous single spike followed by a short plateau, using a holding current of −4 pA and no pulse produced a very short-lived spike in [Ca2+]i followed by decay. (C) A spike train burst, using a holding current of −20 pA and a pulse of 60 pA for 20 s produced a regular stepped rise in [Ca2+]i. (D) A spontaneous spike train burst, using a holding current of −5 pA and no pulse produced a stepped, linear rise in [Ca2+]i.

It was shown that the maximum [Ca2+]i predicted by the model varies depending on the morphology of the AP. Of the various forms shown, high frequency action potential spike trains are the most effective way to build up the intracellular calcium to levels sufficient to initiate contraction over these long periods (Fig. 3C,D). These bursts tend to lead to a step-like morphology in the calcium transient. The 6 second plateau (Fig. 3A) produced a maximum [Ca2+]i 2–3 times smaller than the 20 second bursts (Fig. 3C,D), and the 1 second plateau (Fig. 3B) produced half as much again. The simulated positive correlation between spike frequency and amplitude of [Ca2+]i as shown in Table 2 qualitatively resemble those observed experimentally in the uterus8,27, further validating the model development.

Table 2 Properties of simulated [Ca2+]i transient during voltage clamp.

Shmygol et al.28 demonstrated experimentally that repetitive depolarising current stimuli also produced a stair-case rise in intracellular Ca2+. This phenomenon is also re-produced by the model, and results are shown in Fig. 4, in which the simulated intracellular Ca2+ transient (Fig. 4A) and the underlying ICa (Fig. 4B) during the time course of repetitive stimulus pulses (inset in Fig. 4A) are shown. During the time course of repetitive current stimuli, the simulated stair-case rises in [Ca2+]i matched to experimental data (Fig. 4A), however, the time course of ICa shows notable difference in the amplitude of ICa between simulation and experimental data. In both model and experimental data, the total calcium current exhibits a spike-like behaviour switching between positive and negative spikes. However, the model’s current magnitudes were skewed into the negative compared with experimental data, having an inward current magnitude of approximately double. Note that ICa does not flip between an inward and outward direction. Experimentally the switching spikes are a result of capacitive currents that were not fully compensated for the cell’s capacitance exceeded 100 pF, whereas the amplifier’s range of compensation is limited to 100 pF. In the model the positive spikes are a result of the stimulus pulse causing large depolarising spikes that exceed the calcium channel reversal potentials. Such differences explain the discrepancy between the experimental results and the model in ICa morphology as shown in Fig. 4B. Table 3 shows a range of frequency and magnitude values for ICa, and our simulation data are qualitatively similar to experimental data.

Figure 4
figure 4

[Ca2+]i and ICa were observed under repetitive depolarising pulse train stimulations of −0.4 pA pF−1 at 2.5 Hz. Experimental results were obtained using a depolarising pulse of −80 to 0 mV at 2.5 Hz. The experimental data was recorded at 35 °C using the fluorescent indicator indo-1 in patch-clamped, single uterine myocytes from pregnant rats. A comprehensive methodology can be found in Shmygol et al.28. In order to reproduce the experimental data, the following model parameters from the standard model were multiplied by the factors shown: gCaL × 0.9 and gNa × 6. (A) Myometrial calcium transients exhibit a staircase-like increase. [Ca2+]i ranges from 0.1 to 0.3 μM. A inset: the pulse train protocol. (B) The whole-cell calcium current recorded under voltage clamp shows a negative stair-case phenomenon. Note incompletely compensated capacitive currents. B Inset: a detailed view of the highlighted area.

Table 3 Properties of simulated ICa under repetitive depolarising pulse train stimulations of −0.4 pA pF−1 at variant stimulation rates.

During the repetitive pulse stimulation, the characteristics of the simulated stair-case rise of [Ca2+]i matched to experimental data, which are summarised in Table 4.

Table 4 Comparison of the properties of [Ca2+]i transients between simulated and experimental data under repetitive depolarising pulse train stimulations of −0.4 pA pF−1 at 2.5 Hz.

In the model, the Ca2+ transient has a τrise of 2.9 s, a τdecay of 6.8 s and a magnitude of 0.23 µM (Fig. 4A), which is in close agreement with the data provided by A. Shymgol, which has respective values of 3.4 s, 6.9 s and 0.22 µM. The model’s ICa had a frequency of 2.5 Hz and a magnitude range of 4.3 pA pF−1 (Fig. 4B). This is also in good agreement with the data, which has respective values of 2.5 Hz and 2.75 pA pF−1. All of these show a close match between simulation and experimental data, which further validated the model development.

Modelling drug effects

In order to further validate the developed model, further simulations were conducted to examine if the model could reproduce the action of a variety of hormone and drugs on modulating the electrophysiological action potential, calcium handling and force production.

Simulated effects of oestradiol

Previous studies have shown that in the presence of oestradiol, spontaneous AP bursts and muscle contractions were suppressed26,29. Further studies have shown that the suppressive action of oestradiol may be attributable to its inhibitive actions on the voltage-dependent calcium, calcium-dependent potassium and voltage-dependent potassium channel currents, contributing to the prevention of the genesis of spiking AP29,30. To test the ability of the model to reproduce the effect of oestradiol, the action of oestradiol on its targeted membrane ion channels was simulated by applying the conditions shown in Table 5 and comparing the response to the experimental data of Inoue et al.26.

Table 5 The effects of oestradiol are simulated by scaling down the conductance of the above channels.

Figure 5 shows the simulated spiking APs (Fig. 5A) and the underlying [Ca2+]i (Fig. 5B) in the control (left panels) and oestradiol (right panels) conditions in response to a stimulus current of −0.4 pA pF−1 applied to the model. In the control condition, such a stimulus generated a sequence of spiking APs that were quantitatively similar to the experimental data from Inoue et al. (Fig. 5A), which had a resting potential of ~−50 mV and peak of ~30 mV, with a frequency of ~4 Hz. In response to such repetitive spiking APs, the simulated intracellular Ca2+ concentration accumulated, generated a stair-case trace of [Ca2+]i (Fig. 5B). In the presence of oestradiol, the spiked APs were supressed, resulting in an AP with a long plateau phase. In such a case, the intracellular Ca2+ concentration accumulated to a lesser extent, resulting in a smaller [Ca2+]i amplitude (about 43% reduction in peak calcium) as compared to the control condition. Note that in the presence of oestradiol, the usual staircase-like increase in calcium instead became a smooth, linear increase.

Figure 5
figure 5

Simulated effects of oestradiol. (A) Experimental (black) and simulated (red) membrane potentials under control (left) and oestradiol (right) conditions. Experimentally a 1 s stimulus of 20 pA was applied, with or without 30 µM of oestradiol. In simulations a stimulus current of −0.4 pA/pF was applied to the cell for 1 s. Oestradiol was applied by altering channel conductances (Table 5). (B) Influence of oestradiol on simulated intracellular calcium levels (blue), under control (left) and oestradiol (right) conditions. (C) Effects of oestradiol on total calcium and potassium currents.

The underlying ion channel currents of ICa and IK in control and oestradiol conditions are shown in Fig. 5C. Okabe et al. demonstrated that oestradiol has a similar effect on both the total calcium and total potassium currents29; they reported a difference in reduction of ICa,max and IK,max of 4.4% to 8.4%, with ICa,max undergoing the greater drop. Simulated currents show a difference of 8%, with a greater fall in ICa,max than IK,max (Fig. 5C).

The simulated depressive effects of oestradiol matched to experimental observations. This in one hand validated the model by demonstrating the capability of the model to reproduce experimental observations of oestradiol. On the other hand the simulation results also provided mechanistic insights into understanding the depressive actions of oestradiol, which are attributable to the integral action of oestradiol on ICaL, IK1, IK2 and IK(Ca).

Simulated effects of oxytocin

The myometrium possesses specific receptors for oxytocin, which is a hormone used to induce labour. It has been shown that oxytocin promotes contraction by increasing [Ca2+]i and is up-regulated on the approach to labour31. Studies have also shown it increases both the frequency and amplitude of contractions in rat uteri32,33,34,35.

The mechanism underlying the promoting contraction of oxytocin is incompletely understood. In their study Inoue et al. demonstrated that oxytocin increases the calcium channel currents in rat myometrium, but has little to no effect on the sodium and delayed-rectifier potassium channels, unlike in rabbit where it’s effect on sodium channels is important36. It is unclear if the observed increases in the calcium channel currents can produce the promoting action of oxytocin.

In simulations, we investigated the role of an increased ICaL or ICaT on the electrical and mechanical behaviours of the myometrial cell model. Results are shown in Fig. 6 for control (Fig. 6A) and increased ICaL (Fig. 6B) and ICaT (Fig. 6C). It was shown that up-regulating L-type calcium increases intracellular calcium concentration by 38%, and consequently contractile force by 30%. However it reduces the frequency of the driving action potential by 26%, resulting in less frequent contractions (Table 6).

Figure 6
figure 6

Simulated effects of oxytocin on the membrane potential (top panels), calcium transient (middle panels) and force profile (bottom panels) of a rat uterine cell. A 10 s stimulus of −0.4 pA pF−1 was applied to the membrane, corresponding to a control dose of oxytocin at ~1 nM. (A) Control condition using the standard model. (B) L-type calcium was upregulated by adjusting its conductance by a factor of 1.8. (C) T-type calcium was upregulated by a factor of 10.

Table 6 Simulated effects of oxytocin.

Alternatively, up-regulating T-type calcium increased both the frequency of contractions by 25% and [Ca2+]i,max by 62%. This more closely represents the previous experimental observations32,33,34,35 and suggests that ICaT may play a more significant role than ICaL in the action of oxytocin on rat myometrium. This resembles the behaviour invoked by sodium channel up-regulation in other species and may be explained by the T-type calcium activation (Fig. A2C) being approximately twice as fast as the L-type (Fig. A1C) in the model. The sodium activation is ~8 times as fast as the L-type calcium (Fig. A3C).

Simulated effects of nifedipine

Nifedipine is an L-type calcium blocking agent used to treat pre-term labour, as well as a variety of other conditions relating to muscle contraction or spasms. It is a commonly used tocolytic (anti-contraction medication) favoured for its few side-effects37. It is cheap, given orally and significantly decreases the risk of delivery within 7 days of commencing treatment38. Reports state that the application of nifedipine can reduce the level of intracellular calcium or even eliminate it entirely28,39,40,41.

In order to further validate the model, effects of blocking ICaL on the intracellular Ca2+ concentration were simulated in the condition of a voltage clamp, from −80 to 0 mV for 10 s. Results are shown in Fig. 7. In simulations, blocking ICaL decreases the amplitude of total ICa by approximately 98.7% (Fig. 7A,B), and almost abolished the intracellular Ca2+ transient (Fig. 7C), which closely matched the experimental observation of Wray et al.42 (Fig. 7D).

Figure 7
figure 7

The influence of 10 µM of nifedipine on the total calcium current and Ca2+ transients. A voltage clamp from −80 to 0 mV was applied to the cell membrane for 10 seconds. (A,B) The membrane voltage response to the control and nifedipine conditions respectively, showing a drop of approximately 100 times in the maximum current value. B (inset): Experimental data from Wray et al. demonstrating the drop in current due to the action of nifedipine42. Both simulations (C) and experiments (Wray et al.42) (D) suggest that in the presence of nifedipine, applying a voltage across the membrane prevents an influx of calcium into the cell.

Further simulations were also conducted to investigate the effect of ICaL blocking on the membrane APs and contractile force of the cell model by applying a current pulse of −0.7 pA pF−1 for 0.2 s every 1.5 seconds for 40 seconds. The results are shown in Fig. 8. It was shown that blocking ICaL abolished cell APs (Fig. 8A), resulting in a small contractile force (Fig. 8B), which closely matched the experimental observation of Burdyga et al., who demonstrated that blocking the main entrance for calcium into the cell with nifedipine has the effect of minimising the contractile stress across the membrane (Fig. 8C,D).

Figure 8
figure 8

The effect of 10 µM of Nifedipine on [Ca2+]i and the contractile force of the single uterine cell. The model results were obtained by applying a current pulse of −0.7 pA pF−1 for 0.2 s every 1.5 seconds for 40 seconds. The experimental data was collected from tissue samples by Burdyga et al.7. Each dataset was normalised between its zero and one to allow a comparison of the morphology. (A,B) Experimental levels of intracellular calcium and normalised force measurements, respectively. (C,D) Simulated levels of intracellular calcium and normalised force measurements, respectively. In the control case a rapid increase in Ca2+ is followed by spike-like oscillations, ending in a smooth drop. Calcium accumulation was suppressed by the presence of nifedipine. Nifedipine was simulated by applying a 95% block to the L-type calcium channels.

With a well validated model, we theoretically explored the possible roles of altered Na+ and Ca2+ channel currents in modulating cellular actions potentials.

Experimentally variant action potentials have been recorded from uterine myocytes during different stages of gestation as shown in Fig. 9. Such varied profile characteristics in membrane action potentials may be due to different levels of ion channel expression, which may go through a process of remodelling during the pregnancy22. To test this hypothesis, in Fig. 10, the effect of varied maximal channel conductance of the L-type calcium current gCaL on membrane potential (top), intracellular Ca2+ concentration (middle) and the resultant force (bottom) were simulated. It was shown that under the condition of a current pulse stimulus–reducing gCaL from 0.12 nS pF−1 to 0 (i.e. mimicking block of the channels or down-regulation of the channel) changed the response of the cell model from a burst of spikes to a single spike and plateau, in a sequence that reflects the changes in action potential shape observed during gestation (Fig. 9). Correspondingly, the time trace of the intracellular Ca2+ concentration and the resultant force also showed different profiles and amplitudes.

Figure 9
figure 9

Characteristic spontaneous action potential waveforms recorded at various stages of gestation by Bengtsson et al.22, taken from rat circular muscle tissue at body temperature. Gestational age in days: (A) late 21, (B) early 21, (C) 20, (D) 18–19, and (E) 16–17. The chronology was reversed to more easily compare the change of morphology due to gestation with the effects of channel blocking in Figs 10 and 11.

Figure 10
figure 10

Effects of blocking L-type calcium channels on the AP and excitation contraction coupling. Action potentials (top panels), calcium transients (middle panels) and force profiles (bottom panels) produced by altering the maximum conductance of the L-type calcium current in the model. From left to right the conductance of gCaL (A) 100%, (B) 62%, (C) 57%, (D) 33% and (E) 0% with a 10 second stimulus of −0.4 pA/pF applied.

Similarly, the profile of membrane action potentials is also dependent on the maximal conductance of the sodium current gNa. Figure 11 shows that a reduction in gNa reduced the number of spikes evoked by the injected current pulse. As gNa was reduced the discharge rate reduced from 5/s to a solitary action potential for reductions of >77%, with a minimum maintained discharge rate of 1.4/s at a reduction of 37%. Each spike was accompanied by an increase in [Ca2+]i, and during repetitive spiking these summate. At the end of a burst of spikes the [Ca2+]i, fell with a time constant of ~5.6 s.

Figure 11
figure 11

Effects of blocking fast sodium channels on the AP and excitation contraction coupling. Action potentials, (top panels), calcium transients (middle panels) and force profiles (bottom panels) produced by altering the maximum conductance of the fast sodium current in the model. From left to right the conductance of gNa (A) 100%, (B) 74%, (C) 61%, (D) 57% and (E) 0% with a 10 second stimulus of −0.4 pA/pF applied.

Through the mechanism of electrical and mechanical coupling, changes in [Ca2+]i were expected to be reflected in the resultant contraction force. If L-type calcium channel expression increases as pregnancy progresses then we would expect to see a gradual increase in the levels of intracellular calcium and consequentially a rise in contractile force. Lowering the conductance of the sodium channels had a greater effect on the frequency and amplitude of individual [Ca2+]i and force spikes than the overall amplitude of their plateaus (Fig. 11), whereas lowering calcium channel conductance to 33% was enough to half [Ca2+]i (Fig. 10).

Figure 12 maps out the effect of simultaneous alterations of gCaL and gNa on the intracellular Ca2+ concentration. In the gCaL-gNa parameter space, changing the gCaL and gNa parameters from the bottom left hand corner to the top right resulted in an increase in [Ca2+]i, which is analogous to the increased force shown over the course of gestation43,44. A variant combination of changes in the gCaL and gNa parameters affects the rate of change of cell activity. Increasing sodium channel density had a smaller effect on [Ca2+]i than L-type calcium channel density until both reached ~50%. Changes in gCaL produced a smooth increase throughout the map, whereas increases to gNa in this top-right region resulted in a rapid change. For example, when gCaL is set to 100% increasing gNa from 0% to 50% results in a 10% increase in [Ca2+]i. Increasing it from 50% to 100% results in a larger increase in [Ca2+]i of 30%. This may indicate that an increase in sodium channel density in late pregnancy could produce a phase change in behaviour leading to contractions.

Figure 12
figure 12

The gCaL – gNa parameter-space. The L-type calcium and fast sodium conductances were varied simultaneously in increments of 10% from 0% to 100% to assess the dependence of the [Ca2+]i on the respective channels. The colour map indicates the concentration of peak calcium, relative to the maximum value where both conductances were set to 100%.

Discussion

In this study we present an updated model of the electrical, electro-chemical and mechanical activities of uterine smooth muscle cells based on modifications of the Tong et al. model13. The developed model was validated by its ability to reproduce variant action potential profiles, intracellular Ca2+ transients and generation of the active force, as well as the actions of some pharmacological agents, including oestradiol, oxytocin, nifedipine and wortmannin.

As compared to the previous Tong et al. model, the major advances of the present model are: (i) updated models for some of the membrane currents. The model considers membrane currents of ICaL, ICaT, INa, IK1, IK2, IKA, Ih, INaCa, INaK, INaKCl, ICl, ISOC, INSCC, IK(Ca) ICa-pump and IK-leak. Among them, model equations and parameters for ionic channel currents of ICaL, IK2, IK(Ca) and ICl(Ca) were updated to fit more closely to experimental data; (ii) a more detailed description of the calcium handling system that consists of an SR model to describe calcium sequestration, uptake and release with three currents (SERCA, RyR and a transfer current) and a calcium buffer; and (iii) a more detailed description of cellular force generation. The mechanical model consists of a 4-state kinetic cross-bridge model that describes the cycling between phosphorylated/non-phosphorylated and attached/unattached states, and a 4-component Hill model with passive, active, cross-bridge elastic and visco-elastic forces. As summarised in Table 1, these advances form the most comprehensive rat uterine smooth muscle model, representing a significant improvement as compared to other models for simulating the electrical and mechanical activities of the uterine smooth muscle cells. In addition, the presented model has been comprehensively validated by comparison of simulation data to experimental findings (Figs 25, 7 and 8) of the profiles of action potentials, intracellular Ca2+ transient and force generation in control and in drug actions.

Using the model, we investigated possible mechanisms underlying the transition from a stable steady state to genesis of spontaneous electrical bursting activities during the gestational period. In the absence of external stimuli, the model has a robust and quiescent stable steady state. However, with increased sodium and calcium channel density, the model may show spontaneous bursting activity. In simulations, it was shown that an increase in the sodium channel density appears to be a major factor in changing the model from a quiescent state to a spontaneous bursting state, a main feature of the electrical activity of the uterus in labour (Fig. 12). There is a threshold over which a myometrial cell becomes more sensitive to up-regulation in sodium channel activity, so small changes in sodium channel density would cause significant qualitative changes in cell electrical behaviour. It was also shown that both sodium and calcium channel density play a role in determining the morphology of the action potential and the calcium transient (Figs 10 and 11). While the calcium channel current is more important in determining the amplitude of the calcium transient, the sodium channel current, however, is more important in determining the frequency of AP spikes. This indicates that a combination of changes in the two channels is required to produce the rapid bursting activity associated with the large accumulation of calcium needed for strong, long lasting contractions. These simulation results imply increases in sodium and calcium channel density may happen as a natural process of gestation over time.

The developed model was also used to investigate electrical and mechanical activity of the uterus in various physiological and pharmacological conditions. By altering selected parameters it can reproduce a variety of physiological and pathological behaviours, such as drug effects via pore block or environmental conditions by altering the extracellular ion concentrations. In simulations, the effects of four drugs and hormones on the electromechanical activities of the model were examined and compared to experimental findings for validation purposes. The simulation results for oestradiol (Fig. 5), oxytocin (Fig. 6) and nifedipine (Figs 7 and 8) all closely matched to experimental data. Our results indicate that that action of oestradiol on myometrial activity is attributable to its integral action on the reduction of conductance in multiple channels, including ICaL, IK1, IK2 and IK(Ca); and the action of oxytocin may be attributable to the increased conductance of ICaT. Our results compare favourably to the effects of these drugs seen in the literature, suggesting the developed model may be used as a platform supplementary to experimentation to evaluate the action of a drug or hormone on myometrial electromechanical activities.

Similar to other models, the present model also has some limitations. While it inherits some limitations from the Tong et al. model that has been documented elsewhere, the primary limitation of this model is that it has been constructed based on data available in the literature, which were obtained over many years by different groups under variant conditions and from different species. The data used to build this model comes primarily from voltage-clamp experiments on late-pregnant rat cells. In some cases, where rat data was not available human data was substituted in its place; for example the current density of the T-type calcium channel was fitted to human values using data from α1G/Cav3.1 rat genes45,46. Ideally in future data obtained from the same species and consistent experimental conditions will be helpful to further improve the model development.