A computational model of excitation and contraction in uterine myocytes from the pregnant rat

Aberrant uterine myometrial activities in humans are major health issues. However, the cellular and tissue mechanism(s) that maintain the uterine myometrium at rest during gestation, and that initiate and maintain long-lasting uterine contractions during delivery are incompletely understood. In this study we construct a computational model for describing the electrical activity (simple and complex action potentials), intracellular calcium dynamics and mechanical contractions of isolated uterine myocytes from the pregnant rat. The model reproduces variant types of action potentials – from spikes with a smooth plateau, to spikes with an oscillatory plateau, to bursts of spikes – that are seen during late gestation under different physiological conditions. The effects of the hormones oestradiol (via reductions in calcium and potassium selective channel conductance), oxytocin (via an increase in intracellular calcium release) and the tocolytic nifedipine (via a block of L-type calcium channels currents) on action potentials and contractions are also reproduced, which quantitatively match to experimental data. All of these results validated the cell model development. In conclusion, the developed model provides a computational platform for further investigations of the ionic mechanism underlying the genesis and control of electrical and mechanical activities in the rat uterine myocytes.

T-type calcium current -I CaT There have been reports of I CaT in human myometrial cells [4], [10]. Detailed electrophysiological data of cells expressing rat α1G/Cav3.1 are available [11], [12]. Spontaneous contractions in strips of pregnant rat myometrial tissue were inhibited by the reputed T-type calcium channel blockers mibefradil, NNC 55-0396 (a non-hydrolysable analogue of mibefradil) and Ni + [13], [14]. In addition Ohkubo et al. [15] showed that the expressions of mRNA encoding for the α1G and α1H protein subunits of the T-type calcium channel were gestationally regulated in murine myometrial cells. As a result, this model has been developed from the electrophysiological characteristics of the rat α1G/Cav3.1clona expression cell data, which was recorded at room temperature [11] and adjusted to the human myometrial cell current density [4], [10]. The activation and inactivation steady-state values and the I-V relationships of these different datasets are similar.
I CaT activates at V ≈ -60 mV, peaks between -20 to -30 mV and has fast activation, but inactivation varying between 7 to 100 ms [4], [10]- [12]. This might be due to different external divalent cation concentrations used between experiments.

Calcium pump current -I Capump
The calcium pump (plasma membrane Ca2+ ATP-ase) is an ionic transporter that maintains and regulates the calcium concentrations between the cytoplasm, sarcoplasmic and endoplasmic reticulum (SR and ER) in the cells [16]- [18]. There is a disparity between the concentrations of calcium ions in the cytosol (0.1 μM) and the ER and SR (0.1 to 1 mM) [19]; thus Ca 2+ ions must be pumped against the concentration gradient, using energy obtained from the ATP chemical energy store. We assume an instantaneous pump with quick bindings resulting in rapid equilibrium using an empirical Hill-type formula.

Sodium/potassium exchanger current -I NaK
This exchanger current plays a major role in regulating intra and extracellular sodium and potassium cations and has been found in late pregnant rat [20], [21] and human [1] myometrial cells. A net negative charge within the cell is achieved by using energy derived from hydrolysed ATP to expel three Na + ions for each two K + ions brought in. I NaK reportedly increases during pregnancy in the rat myometrium [21].
This model adopts and modifies the formulae from the Luo-Rudy mammalian ventricular model for I NaK , due to the absence of information relating to its role in the myometrium [22]. The inhibition of the exchanger changes dynamically, affecting the activities of K+ and its dependent currents, leading to a bursting AP.
Fast sodium current -I Na I Na has been found in pregnant human myometrial cells [23] and a fraction of late pregnant rat myometrial cells [4]. This channel is sensitive to tetrotodoxin and sagitoxin. It may be responsible for the fast upstroke of AP generation, the formation of the gap junction and quick excitation propagation during parturition [24]. Studies have found that the number density of sodium channels in the rat increases towards term [25].
Sodium/calcium exchanger current -I NaCa I NaCa expels a Ca 2+ ion and brings in three Na + ions against the concentration gradient. Abe [26] found that the prepotential and spontaneous spike discharge disappears at low [Na + ] o , but the tissue remains excitable; a depolarising current stimulus causes larger, faster spikes. This confirms the existence of the electrogenic I NaCa in the rat myometrium. We also borrow from cardiac models for this current, particularly the Hund-Rudy canine ventricular model [27], to make up for the lack of available data. Fig A2. T-type calcium current -I CaT A and B: the activation and inactivation steady states of I CaT , respectively. This is not an ideal fit to the steady state data, but produced a closer reproduction of the I CaT I-V relationship. C and D: the activation and inactivation time constants of I CaT , respectively. E: the normalised current-time traces resulting from a voltage-step protocol of -70 mV to 50 mV in steps of 10 mV with a holding potential of -80 mV, used to establish the I-V relationship of I CaT . F: the normalised peak I-V relationship of I CaT . Fig A3. Fast sodium current -I Na A and B: the activation and inactivation steady states of I Na , respectively. C and D: the activation and inactivation time constants of I Na , respectively. E: the normalised current-time traces resulting from a voltage-step protocol of -50 mV to 50 mV in steps of 10 mV with a holding potential of -90 mV, used to establish the I-V relationship of I Na . F: the normalised peak I-V relationship of I Na .
Voltage-dependent potassium currents -I K1 and I K2 I K1 and I K2 play a role in repolarising the cell. We have altered the activation steady state relationship of I K2 from the Tong et al. model in order to provide a better fit to the I-V relationship obtained by Knock and Aaronson ( Fig  A5A).
Knock et al. [28] separated I K1 and I K2 myometrial potassium currents and categorized them by their inactivation properties and sensitivity to pharmacological blockers of varying channel subtype specificity [28], [29]. I K1 and I K2 have rectifying properties and were found in myometrial cells of late pregnant rats [29] and humans [28]; they have very slow dynamics in comparison to other membrane currents and are unevenly distributed between cells. Their primary role appears to be to reset the action potential to its resting state by repolarising the cell membrane. I K1 and I K2 have similar voltage-dependent kinetics. I K1 activates between V ≈ -50 to -40 mV and half-inactivates between -62 to -68 mV. I K2 activates between V ≈ -40 to -30 mV and half-inactivates between -30 to -20 mV. I K2 inactivates faster than I K1 . Our model took advantage of the more abundant information on the electrophysiological characteristics of human myometrial I K1 and I K2 , complementing them with data from rat myometrium. I K1 and I K2 were separated by their activation thresholds, inactivation properties and current sensitivities to 4-aminopyridine (4-AP) and TEA. The Tong et al fit of the IK2 activation steady state was adjusted to allow a better fit to the unpublished I-V curve data from Knock and Aaronson.
A-type transient potassium current -I KA I KA has very fast activation and inactivation kinetics, is found in both rat and human myometrial cells (though it is rare in pregnant cells), is sensitive to 4-AP and insensitive to TEA [30]. It activates at V ≈ -40 mV, peaks within ≈ 10 ms and is almost completely inactivated within 50 ms [2], [30]. In human myometrium it has a half-inactivation ≈ -70 mV and a slope factor ≈ 5 mV [28], [30]. These are similar characteristics to the transient potassium current in myometrial cells isolated from immature rats.
Calcium activated potassium current -I K(Ca) I K(Ca) is associated with setting the resting membrane potential and regulating the shape and duration of the action potential. It may be responsible for up to 10% of the total K + current in late pregnant rate myometrial cells [29]. Quantitative data for the dynamics of this channel in the myometrium is lacking and the equations developed for rat mesenteric smooth muscle have been used [31]. Calcium activated potassium currents provide a means of burst termination in pancreatic cells (Chay Keizer model) and neurones (Chay Fan model).

Hyperpolarisation-activated current -I h
This current has reportedly been found in the circular and longitudinal cells of the pregnant murine myometrium [32], [33] and is otherwise known as the funny current (I f ) in cardiac cells. It is activated by hyper polarisation of the membrane potential, below -60 mV, has a half-maximal activation around -85 mV and does not inactivate. In voltage clamp experiments the activation is slow, taking > 1 s. K + ions permeate this channel more readily than Na + ions, with a ratio of P Na /P K = 0.35. The reversal potential E h ≈ -20 mV is dependent on the extracurricular concentration of K + and Na + ; and excess of K + increases the current amplitude and a low measure of Na + decreases it, while Cs + blocks the channel. This current is activated at the resting membrane potential and its role is to maintain this potential and spontaneous activity within the myometrium.
Data for this channel has been found experimentally by Okabe et al. [32] and Satoh [33]. Temperatures ranged from 24-30 o C; the time constant τ h was adjusted to compensate, with a Q 10 value of 3.5. The single cell model was developed for a longitudinal cell. The circular cell half-activation has been shifted by -15 mV to match the longitudinal data. At room temperature, with a holding potential of -50 mV, I h has a current density of -1.03 pApF -1 .   Calcium-activated chloride current -I Cl Spontaneous transient inward currents (STICs) are pro-contractile currents that exhibit the characteristics of calcium-activated channels (CaCCs) and have been observed in murine and human uterine smooth muscle cells, and in putative pacemaking interstitial cells of Cajal in the intestine [34]. They are blocked by Cl-inhibitors such as niflumic acid and flufenamic acid and are greatly diminished in calcium-free buffers, consistent with their being calcium-activated channels ( CaCCs) ], showing calcium activation, voltage dependence and a reversal potential of -15 mV [35], [36]. The anoctamin 1 and 2 (ANO1, ANO2) channels are responsible for these currents [36], and their kinetics were formulated from experimental data in Jones et al 2004 [9].
It is believed that CaCCs prolong the duration of contraction in pregnant tissue [9], [37]. They may play a role in depolarising the membrane, leading to increased Ca entry [38], [39]. In the interstitial cells of Cajal (ICC) of the intestines they play a role in pace making; ICC-like cells have been described in the myometrium [40] and Bernstein et al suggest there may be evidence of this role in the myometrium [36]. The contribution of I Cl to the total inward current is a matter for debate, with contrasting reports indicating it is responsible for up to 30% [9] or unimportant [41].
Jones et al [9] assessed the electrophysiology of I Cl in the single cell using two distinct voltage clamp experiments. The activation time relationship of I Cl is relatively complex. In order to provide a more representative activation model for I Cl , its time coefficient relationship has been split into two functions ( Fig A8C); the current is calculated differently for voltages above and below 0 mV. The Tong model used an alternative continuous function for the activation time constant and fitted I-V data from Jones 2004. I-V curves for I Cl for our model reproduce those seen in the more recent Bernstein 2014 [36], [42]- [44]. the normalised current-time traces resulting from a voltage-step protocol of -60 mV to 40 mV in steps of 10 mV with a holding potential of 0 mV, used to establish the I-V relationship of I Cl . E and F: the peak and tail I-V relationships of I Cl respectively. This channel was reformulated to account for new data describing the peak I-V values [36].

Ionic concentration dynamics
[Ca 2+ ] i increases as the calcium channels respond to a depolarising cell membrane. Na + , K + and Clion regulation are also regulated by the intracellular ion handling component. Additional Ca 2+ is also released by the sarcoplasmic reticulum store (SR), by the process of calcium-induced calcium-release (CICR).The [Ca 2+ ] i -dependent kinetics respond to this influx of Ca 2+ into the cytoplasm, leading to the development of contractile forces. Calsequestrin (CSQN) binds to Ca 2+ in the SR, storing it despite the higher concentration within the SR relative to the cytosol.
The Tong et al. model calculates the calcium current using a simple formula consisting of three components: contributions from the L-type calcium current, the T-type calcium current and the sodium-calcium exchanger. Notably absent from this description are contributions from the non-selective cation current (NSCC), the storeoperated non-selective cation current, the calcium pump and the SR currents. We included these factors for a more complete description of the calcium dynamics (B17.1).
Another advance from the previous Tong et al. model was the inclusion of formulae describing the sodium, potassium and chloride ion currents. Changes in the concentration of these ions contribute less significantly to contraction than calcium ions; however their effects are not negligible and should be accounted for in a comprehensive model. We modelled the changes to intracellular ionic concentrations of Ca 2+ , Na + , Kand Clwith a series of first order differential equations (B17.2, B17.4, B17.6, and B17.8, respectively). A calcium buffer was implemented to reduce sensitivity to [Ca 2+ ] i (B17.2).
The calcium currents include contributions from the L-type calcium current, the T-type calcium current, the sodium-calcium exchanger, the non-selective cation current, the store-operated non-selective cation current, the calcium pump and the SR currents in (B17.1). Changes in the concentration of Na + , Kand Cldo not directly contribute to contraction, but they do influence the behaviour of several ionic currents: the fast sodium current, the voltage-dependent potassium currents, the A-type transient potassium current, the calcium-activated potassium current, the background potassium leakage current, the hyperpolarisation-activated current, the sodium/potassium exchanger current, the sodium/calcium exchanger current, the sodium-potassium-chloride co-transport current, the store-operated non-selective cation current and the non-selective cation current. The feedback from changes to these currents can in turn affect the morphology of the AP, calcium transients and ultimately contractile force. The intracellular ionic concentration dynamics are measureable in in vitro experiments and provide additional constraints on the behaviour of the cell model [45].

Cell contraction
In order to model the contractile activity in the cell we must couple the electrical systems to mechanical systems (Fig 1), using the intracellular calcium concentration as a mediator.
Calcium plays a key role in the mechanical contraction and relaxation of myometrial tissue. Ca 2+ enters the cell through voltage-gated channels, primarily through the L-type Calcium channel. It can also be released from its store in the sarcoplasmic reticulum ( Fig A9) via IP 3 and the ryanodine receptors. Some experimental evidence suggests the SR makes a minimal contribution [46]- [48] to contractile activity; though other findings suggest it may be more significant [49], [50]. Optical measures of intracellular [Ca 2+ ] in tissue preparations have been used as a surrogate measure of myometrial cell excitation [51].

The sarcoplasmic reticulum (SR)
The SR is a calcium store within the cell. Noble et al state that the myometrial SR contributes to sustained depolarisation, bursts of APs and phasic contractions [50]. They remark that the model of the SR developed in cardiac muscle, whereby a small amount of L-type Ca 2+ entry triggers a substantial rise of cytoplasmic Ca 2+ via RyRmediated SR Ca 2+ release leading to contraction is not substantiated by their findings. Rather the Ca 2+ released from the SR stimulates store-operated calcium entry (SOCE), meaning that although the SR plays an important role, the greater contribution to contraction in myometrial cells comes from the membrane channels, unlike cardiac myocytes.
We modelled the SR as a two-compartment system; an uptake compartment and a release compartment, as described for vascular smooth muscle cells by Yang et al [52]. Calcium enters the SR via the SERCA pump (I up ); the release mechanism is governed by a Hai and Murphy four-state ryanodine receptor model (Fig A10) via the process of CICR (I rel ) [53]. The three currents I up , I tr and I rel (B18.1 -B18.3) describe the uptake into the SR, transfer between compartments and release of calcium ions into the cytoplasm.
The Michaelis constant of 1µM suggested by Kapela et al [31] was used rather than the 80 µM used by Yang et al, as it is closer to physiological values. To avoid excessive loading of the SR a leakage component was included, R leak . This and the release time constant τ rel were adjusted from the Kapela et al model to better fit myometrial [Ca 2+ ] i data. F rel is a scaling factor used to adjust the magnitude of the release current.

The ryanodine receptors (RyR)
In the process of CICR [Ca 2+ ] i binds to the RyR; they open, allowing calcium from inside the SR store to enter the intracellular space. In order to simulate CICR we have incorporated the model described by Yang et al [52]; the states of the ryanodine receptors are R 00 (free receptors), R 10 (activation sites with bound calcium), R 01 (inactivation sites with bound calcium) and R 11 (activation and inactivation sites with bound calcium) [52]. Transitions between these states were controlled by the activation and inactivation rate constants, K ±r1 and K ±r2 .

The kinetic cross-bridge cycle (KCBC)
We modelled the kinetic cross-bridge cycle using the four-state hypothesis described by Hai and Murphy [53] (Fig  A10). Unlike striated and cardiac muscles, smooth muscles do not contain troponin for the calcium to bind to; instead [Ca 2+ ] i binds with calmodulin and activates the myosin light-chain kinase (MLCK), forming cross-bridges between actin and myosin filaments, via the attached phosphorylated and attached de-phosphorylated non-cycling cross-bridges (AMp and AM) [55], [56], leading ultimately to muscle contraction. There are also free cross-bridges and phosphorylated cross-bridges (M and Mp). The flow between states was determined by the rate constants, K 1-6 , which are [Ca 2+ ] I dependent. K 7 was the slow detachment rate constant, which leads to relaxation.
Hai and Murphy established the values for the rate constants by first fitting phosphorylation kinetic data for K 1 , K 2 , K 5 and K 6 . These were then held constant while K 3 , K 4 and K 7 were fitted to stress data. The rate-limiting step was revealed to be K 7 , the latch-bridge detachment rate.  [53] and Yang et al [52]. The states are defined as being bound or unbound to Ca 2+ and if bound, whether the calcium is attached to activation sites, inactivation sites or both. The rate constants K r1 and K r2 represent the speed that Ca 2+ binds to the activation and inactivation sites, respectively. Unbinding rates are indicated with a minus sign. Right: Four-state kinetic cross-bridge cycling model in smooth muscle adapted from Yang et al [52]. the four-state kinetic cross-bridge cycling model in smooth muscle adapted from Yang et al [52]. The states are differentiated by their status as attached or detached, and phosphorylated or dephosphorylated. Cycling between the states is described by the various rate constants, K.

Force production
We modelled force production in the cell as described by Yang et al (Fig A11) [52], using a modification of the three-component Hill model. The cross-bridges in smooth muscles are arranged in parallel with equal spacing; hence force production by the muscles is proportional to the length of the myosin fibre. The total myometrial force is comprised of four components; they are the passive force F p , cross-bridge elastic force F x , active force F a and series visco-elastic force F s . [52]. The four force components F p , F a , F x and F s couple together to generate tension and stress in the cell. The k values are spring constants and μ s is the coefficient of viscoelasticity of the dashpot, tasked with absorbing shock from impacts caused by the spring. F p (B21.5) is dependent on the length of the cell l c , taken to be 123 μm. The myometrial cell length for the late pregnant rat can range from 80 to 700 μm [3], [25], [57], [58]. Other parameters necessary to calculate F p that were not available experimentally for the uterine cell have been substituted for values found in the Yang et al arterial smooth muscle model. F x (B21.6) is dependent on the length of the cross-bridge l x and the length of the active force component l a . l opt represents the optimal muscle length, i.e. when the maximum force is achieved from the overlapping of the myosin and actin filaments. F x and F a coupled to describe the sliding of the myosin and actin filament. F a (B21.7) is the component of the force generated by the activity of the cross-bridges. It is characterised by the friction coefficients for the attached phosphorylated (AMp) and attached dephosphorylated (AM) cross-bridges. F s (B21.8) is described by a spring and dashpot working in parallel. The force produced by the spring can deform the element and so the dashpot is included as a shock absorber to compensate for the effect. The displacement for the two components is the same and the total force produced by this element is a combination of the elastic and viscous forces. The total force produced by the cell is the sum of F a and F p (1.5).

Supplement B -Formulations and parameters
This supplement contains all of the equations and parameters used in our single cell model. Formulae that are not present in Tong et al or have been altered are highlighted in red.

mM
Half saturation of Na + in I NaK 2 Hill coefficient of Na + in I NaK Parameter Value Definition 3.58 x 10 -10 nmol -2 J -1 s -1 cm 2 Cotransport coefficient 1 Cotransport regulation    Cross-bridge cycling velocity 1.3 μN ms μm -1 Friction constant for phosphorylated cross-bridges 85.5 μN ms μm -1 Friction constant for latch bridges 0.01 μN ms μm -1 Viscosity coefficient of series element 4.5 Length modulation for series viscoelastic element 0.1 Length modulation for passive element 7.5 Length modulation constant for active and cross-bridge elements