uterine contractility is generated by contractions of the myometrial smooth muscle cells (SMCs) that compose most of the myometrial layer of the uterine wall. In the nonpregnant uterus, synchronous contractions of these SMCs produce changes in the geometry of the uterine fluid-wall interface. These changes induce intrauterine fluid motions that are essential during early phases of reproduction (3, 11, 28). During parturition, the synchronized contraction of these myocytes generates the forces required to deliver the baby out of the uterus. Depolarization of the cell membrane initiates calcium ion (Ca²⁺) entry into the cells through voltage-operated Ca²⁺ channels (VOCCs) and thereby a rise in the intracellular Ca²⁺ concentration (C_Ca,i). The elevated level of C_Ca,i allows binding of Ca²⁺ and calmodulin, thus activating myosin light-chain kinase (MLCK), which phosphorylates a regulatory myosin light chain (MLC) (29, 32). This subsequently allows the formation of cross bridges between actin and myosin filaments and the generation of muscle contraction.

The excitation-contraction process was studied in both rat and human myometria using the voltage-clamp technique. Stimulation of isolated myocytes using voltage pulses revealed the current-voltage relationships of the pregnant myometrial SMCs in rats (18, 34) and in humans (7, 37). Application of single voltage pulses demonstrated that the Ca²⁺ current (I_Ca) through L-type VOCCs significantly increased C_Ca,i, whereas repetitive stimulation with pulse trains revealed that both VOCC opening and Ca²⁺-induced Ca²⁺ release (CICR) from the sarcoplasmic reticulum (SR) are responsible for the increased C_Ca,i (18). The entry current following depolarization of rat myometrial cells comprises two components, which were identified based on differences in their activation and inactivation properties as well as in their kinetics. The first component is a fast Na⁺ current (I_Na), while the second is a slow I_Ca (34). Studies of the mechanisms responsible for the decay of C_Ca,i in the pregnant rat myometrium, which is a critical process for SMC relaxation, showed that Ca²⁺ pumps in the plasma membrane are responsible for 30% of the Ca²⁺ extraction from the cell, and Na⁺/Ca²⁺ exchangers are responsible for up to 60%. The remaining Ca²⁺ is probably handled by the intracellular stores (19). Accordingly, the sarcolemmal mechanisms of Ca²⁺ extrusion are crucial for C_Ca,i decay, whereas the significance of Ca²⁺ intake into the intracellular stores is lower.

A different group of studies was concerned with the relationship between C_Ca,i, MLC phosphorylation, and the contractile force of myometrial SMCs. The first study, performed in cultured human myometrial cells, revealed that an increase in C_Ca,i initiated an increase in MLC phosphorylation (12). Simultaneous measurements during spontaneous and electrically induced contractions of nonpregnant human myometria showed that the force develops at a slower rate than the C_Ca,i increases and the MLCs are phosphorylated (30). Maximal and steady-state values of MLC phosphorylation were reached prior to these values of C_Ca,i. It was suggested that desensitization of MLCK by phosphorylation after the first second following initiation of contraction causes a decrease in the rate of MLC phosphorylation and force production. Similar findings were obtained when the steady-state and transient C_Ca,i vs. force relationships were investigated in strips of human pregnant myometria during spontaneous and agonist-induced contractions (22). Experiments with strips of pregnant and nonpregnant human myometria also revealed that MLC phosphorylation during contraction was lower in pregnant compared with nonpregnant tissue, while the amount of generated stress per MLC phosphorylation level was higher in the pregnant myometrial tissue (29).

Several mathematical models have been developed to describe the control of C_Ca,i level, force production, and length changes in various types of SMCs. The Ca²⁺ transport system and membrane potential in a single arterial myocyte were modeled by two coupled oscillators that simulated the interaction between intracellular C_Ca,i and membrane potential due to cyclic release of Ca²⁺ from internal stores and cyclic influx of extracellular Ca²⁺ (15). The stress produced by the cell was calculated using the four-state cross-bridge model of Hai and Murphy (5), which relates C_Ca,i to cross-bridge formation and stress development. The process of excitation-contraction in a cerebrovascular SMC was also described using a mathematical model (33). The electrochemical behavior as well as the regulation of C_Ca,i were described using a Hodgkin-Huxley-type membrane model combined with a fluid compartment model. Material balance equations were used to calculate the concentration of various ions within the cytosol. For calculations of C_Ca,i the effects of buffering and Ca²⁺ fluxes through the membrane of the SR were also taken into account. In both models, equations describing various ionic membrane currents were established from their characteristic activation curves. The model parameters were set according to measurements performed in isolated cells, using voltage-clamp techniques and additional experimental data.

The existing models of the myometrium described uterine behavior at the tissue and organ level and predicted contractile forces that closely resembled clinical measurements of normal intrauterine pressure during contractions in labor (2, 27, 35). However, the relations between membrane depolarization, Ca²⁺ control, and force production at the level of a single myometrial myocyte were not described in detail. Accordingly, the present model was developed to simulate the complete process of a single myometrial smooth muscle contraction, which is initiated by depolarization. The model is based on the electrophysiological properties of the cell and the cellular mechanisms that relate the rise in C_Ca,i to stress production and is used to study the operation and properties of these mechanisms. The predicted variations in C_Ca,i, MLC phosphorylation, and stress produced by the contracting myocytes are compared with experimental data from human and rat myometrial cells.

GLOSSARY

CaM kinase 2	Ca2+/calmodulin-dependent protein kinase II
CICR	Calcium-induced calcium release
HH	Hodgkin-Huxley
MLC	Myosin light chain
MLCK	Myosin light chain kinase
MLCP	Myosin light chain phosphatase
MSE	Mean square error
SMC	Smooth muscle cell
SR	Sarcoplasmic reticulum
VOCCs	Voltage-operated calcium channels

Process 1

ECa	Ca2+ reversal potential (mV)
R	Ideal gas constant (mJ·mol−1·K−1)
T	Temperature (K)
F	Faraday constant (C/mol)
CCa,ex	Extracellular calcium concentration (mM)
CCa,i	Intracellular calcium concentration (mM)
IVOCC	Ca2+ current through L-type VOCCs (pA)
gCa	Ca2+ conductivity of L-type VOCCs (nS)
Vm	Membrane voltage (mV)
JVOCC	Ca2+ influx through L-type VOCCs (mM/s)
Vcell	Cell volume (liters)
ZCa	Valence of Ca2+
gm,Ca	Maximal Ca2+ conductance (nS)
ρVCa	Relative amount of open L-type VOCCs
VCa,1/2	Half-activation potential of L-type VOCCs (mV)
KCa,1/2	Slope of activation function at VCa,1/2 (mV)

Process 2

JCa,pump	Ca2+ efflux through Ca2+-ATPase pumps (mM/s)
Vp,max	Maximal flux through Ca2+-ATPase pumps (mM/s)
Kph	Ca2+ concentration in half-activation of Ca2+-ATPase pumps (mM)
n	Hill coefficient of Ca2+-ATPase pumps
JNa/Ca	Ca2+ flux through Na+/Ca2+ exchangers (mM/s)
GNa/Ca	Maximal flux through Na+/Ca2+ exchangers (mM·s−1·mV−1)
KNa/Ca	Ca2+ concentration in half-activation of Na+/Ca2+ exchangers (mM)
CNa,ex	Extracellular Na+ concentration (mM)
CNa,i	Intracellular Na+ concentration (mM)
ENa	Na+ reversal potential (mV)
Vm,Na/Ca	Reversal potential of Na+/Ca2+ exchangers (mV)

Process 3

FM	Fractional amount of free unphosphorylated cross bridges
FMp	Fractional amount of free phosphorylated cross bridges
FAMp	Fractional amount of attached phosphorylated cross bridges
FAM	Fractional amount of attached dephosphorylated cross bridges
K1	Rate constant of phosphorylation of M to Mp (1/s)
K6	Rate constant of phosphorylation of AM to AMp (1/s)
K2	Rate constant of dephosphorylation of Mp to M (1/s)
K5	Rate constant of dephosphorylation of AMp to AM (1/s)
K3	Rate constant of attachment of fast cycling phosphorylated cross bridges (1/s)
K4	Rate constant of detachment of fast cycling phosphorylated cross bridges (1/s)
K7	Rate constant of latch bridge detachment (1/s)
CCa,1/2MLCK	CCa,i at half-activation of MLCK (mM)
nM	Hill coefficient of MLCK activation by Ca2+/calmodulin

THE MODEL

Description of the Model

The biophysics of excitation-contraction in a uterine myocyte consists of three simultaneous processes: Ca²⁺ entry into the cell in response to membrane depolarization, Ca²⁺ extraction from the cell, and myocyte contraction (Fig. 1).

Process 1 (Eqs. 1–3).

It is assumed that Ca²⁺ entry following membrane depolarization occurs almost entirely via L-type VOCCs (31). This process is controlled by the membrane voltage and is stopped when membrane repolarization inactivates the channels. The contribution of T-type Ca²⁺ channels to Ca²⁺ entry into the myocyte can be neglected for a resting membrane voltage in the range of −50 to −35 mV since these channels are inactivated in a membrane voltage higher than −60 mV. In addition, these channels are thought to account for the recruitment of L-type channels, whereas the L-type channels contribute most of the inward Ca²⁺ current (13). The time-dependent properties of L-type VOCCs are not taken into account, and the Ca²⁺ entry through the channels is assumed to be solely dependent on the membrane voltage, or, equivalently, the time-dependent changes in the conductivity of the channels are assumed to be small relative to the voltage-dependent changes. This assumption is reasonable for brief depolarizations of the cell membrane, when the quantity of interest is the amount of Ca²⁺ entering the cell following the depolarization and not the time-dependent changes in the Ca²⁺ current through the channels. Since inactivation of the channels is not taken into account, it is possible that the model will overestimate Ca²⁺ entry into the cell at late stages of the response to depolarizations or action potentials.

Ca²⁺ release from the SR significantly changes the shape of C_Ca,i variations only for complex types of stimulation, and, accordingly, it is not explicitly taken into account at this stage. The underestimation of C_Ca,i caused by this assumption may be compensated for by an overestimation of Ca²⁺ entry through the L-type VOCCs due to neglect of time-dependent inactivation. The effects of these assumptions on the resulting C_Ca,i rise are examined in the results (Control of Intracellular Ca²⁺Concentration), where it is shown that the overall prediction of C_Ca,i is fairly accurate.

Process 2 (Eqs. 4 and 5).

Two mechanisms of Ca²⁺ extraction out of the cell are assumed to account for the C_Ca,i decay from elevated levels: Ca²⁺ pumps and Na⁺/Ca²⁺ exchangers in the plasma membrane. Since mechanisms responsible for Ca²⁺ intake into the SR have been shown to play only a minor role in C_Ca,i reduction (19), these are not taken into account. The Ca²⁺ pumps extract Ca²⁺ from the cell when C_Ca,i is high. The direction of the Ca²⁺ flux through the Na⁺/Ca²⁺ exchangers is set by the difference between the membrane potential and the reversal potential of the exchangers. The reversal potential, in turn, is set by C_Ca,i and the intracellular concentration of Na⁺ (C_Na,i). Below the reversal potential, the exchangers remove Ca²⁺ from the cell; above the reversal potential, Ca²⁺ is introduced into the cell, as long as the affinity of this mechanism to C_Na,i allows its operation in reverse direction. Since Na⁺ control mechanisms are not included in the model, the calculation of changes in C_Na,i as a function of time is beyond the scope of the model. Accordingly, C_Na,i is assumed to be constant. During voltage-clamp experiments, the cells are superfused with solutions containing 2 mM Ca²⁺ (18, 19). This extracellular Ca²⁺ concentration (C_Ca,ex) is much higher than the increase in C_Ca,i during contraction, which is <1 μM (18). Accordingly, C_Ca,ex is assumed to be constant throughout the simulations and unaffected by Ca²⁺ entry into the cell. Similarly, the extracellular concentration of Na⁺ during simulation is assumed to be 140 mM and constant, as used for superfusion of the cells (18, 19).

Processes 1 and 2 control the level of C_Ca,i, as described by Eq. 6.

Process 3 (Eqs. 7–9).

Stress production by the contracting cell is described using the four-state cross-bridge model of Hai and Murphy (5), which was developed for arterial and tracheal SMCs. The model describes the fractional amount of myosin phosphorylation and cross-bridge formation as well as myosin dephosphorylation and the formation of latch bridges. Thus, it allows calculation of the stress produced by the cell. These stresses are assumed to be proportional to the forces produced by the cell during experiments with either single cells or strips of myometrial tissue.

Governing Equations

Process 1: opening of L-type VOCCs.

The Ca²⁺ current through all the L-type VOCCs of the cell may be computed from the transmembrane voltage (V_m) and the channel conductivity. This is described by the following equation:

where g_Ca is the conductivity of L-type VOCCs in the cell membrane described by Eq. 3 and E_Ca is the Nernst potential for calcium. This current is directly related to the Ca²⁺ influx (J_VOCC). A positive current indicates ion flow from the cytoplasm to the extracellular space; accordingly, the flux of ions into the cell is given by (15): where Z_Ca = 2 is the valence of Ca²⁺, V_cell is the cell volume, and F is the Faraday constant. The Nernst potential for Ca²⁺ is given by (8): where R is the ideal gas constant and T is the absolute temperature (K).

The peak inward current through all the L-type VOCCs of the cell following depolarization from a given holding potential shows a U-shaped relation to the depolarization voltage. Accordingly, the voltage dependence of the Ca²⁺ conductivity can be described using a Boltzmann-type activation curve (15, 19):

where g_m,Ca is the maximal Ca²⁺ conductance and ρ_VCa(t) is a function of V_m(t) as described by the following activation function: where V_Ca,1/2 is the half-activation potential of L-type VOCCs (i.e., the potential at which half of the channels are open) and K_Ca,1/2 is the slope of the activation function at V_Ca,1/2.

Process 2: extraction of Ca²⁺ from the cell.

The efflux of Ca²⁺ through Ca²⁺ pumps in the plasma membrane can be described using Hill functions as follows:

where n is the Hill coefficient, K_ph is the Ca²⁺ concentration for half-activation of the pump, and V_pmax is the maximal velocity (in mM/s) of Ca²⁺ extraction by the pump (8, 15).

The Na⁺/Ca²⁺ exchangers can either extract Ca²⁺ from or insert Ca²⁺ into the cell. In normal physiological conditions, the exchangers mostly extract Ca²⁺ from the cell. The flux of Ca²⁺ through all the Na⁺/Ca²⁺ exchangers in the cell can be described by the following equation (15):

where G_Na/Ca is the conductance of all Na⁺/Ca²⁺ exchangers in the cell and K_Na/Ca is the Ca²⁺ concentration in half-activation of the Na⁺/Ca²⁺ exchangers by Ca²⁺. The reversal potential of the Na⁺/Ca²⁺ exchanger (V_m,Na/Ca) is given by The Nernst potential for sodium (E_Na) is given by where C_Na,ex is the extracellular Na⁺ concentration and C_Na,i is the intracellular Na⁺ concentration. The extracellular Na⁺ concentration is set according to the conditions of the simulated experiment. The internal Na⁺ concentration was optimized to fit experimental results.

The net concentration of C_Ca,i is controlled by Ca²⁺ influx and efflux from the cell; thus,

In this equation, the transient effects of Ca²⁺ buffering by the superficial SR are not taken into account. The steady-state effect of buffering is included in the equation, but not explicitly. It can be described by a multiplication factor of the right side of the equation, in the range of (0, 1) (21). Accordingly, only a fraction of the Ca²⁺ entering the cell increases C_Ca,i, whereas the rest undergoes uptake by the buffering proteins (20, 36). Since each of the fluxes in the right side of the equation are scaled by a characteristic conductivity that was set by fitting to experimental results, an explicit multiplication by a buffering factor is redundant. The effect of buffering is accounted for by low conductivities that can be interpreted as higher conductivities multiplied by a small buffering factor.

Process 3: stress production by the contracting SMC.

The four-state cross-bridge model of Hai and Murphy (5) for description of the kinetics of myosin phosphorylation and stress development consists of the following four species of cross bridges representing functional states: M = free unphosphorylated; M_p = free phosphorylated; AM_p = attached phosphorylated; and AM = attached dephosphorylated (latch).

The governing equations describe the changes in the fractional amount (F) of each of the cross-bridge species as a function of time (5):

These equations are solved under the constraint where K₁–K₇ are the seven rate constants of the chemical reactions involved in muscle contraction. K₂ and K₅ describe the rates of dephosphorylation of M_p to M and AM_p to AM by MLCP. K₃ and K₄ represent the rates of attachment and detachment of fast cycling phosphorylated cross bridges. K₇ is a rate constant for detachment of latch bridges and thus must be lower than K₄. K₁ (t) and K₆ (t) are rate constants of phosphorylation of M to M_p and AM to AM_p by the active MLCK-Ca²⁺/calmodulin complex. These constants are regulated by Ca²⁺ and therefore vary with time. MLCK is activated through binding of Ca²⁺ to calmodulin followed by binding of the Ca²⁺/calmodulin complex (Ca²⁺/calmodulin) to MLCK. However, sensitization of MLCK to further increase in C_Ca,i occurs because of phosphorylation of this enzyme by Ca²⁺/calmodulin-dependent protein kinase II (CaM kinase 2) (30). As a result, MLC phosphorylation is rapid at the initial stages of the contraction, when C_Ca,i is about twice as high as its resting level, and becomes more moderate when C_Ca,i continues to increase to a level four times higher than its resting level. Accordingly, the rate constant of myosin phosphorylation, K₁, is higher at the initial stages of the contraction and lower at advanced stages. The desensitization process of MLCK is not considered explicitly in the present model. Activation of MLCK is set by C_Ca,i, and, accordingly, MLC phosphorylation initially increases rapidly and then reaches a plateau later during the contraction. The dependence between the rate constant of MLC phosphorylation (K₁) and C_Ca,i is described by the following equation (15): where C_Ca,1/2MLCK is the C_Ca,i required for half-activation of MLCK by Ca²⁺/calmodulin and n_M is the Hill coefficient of activation. The fractional amounts of cross-bridge species are used to calculate myosin phosphorylation and stress production. The amount of phosphorylated myosin relative to the total amount of myosin within the cell is given by The stress produced by the muscle, relative to the maximal possible stress, is given by

Model Parameters

The parameters of the model can be divided into three groups. The general parameters relevant for all types of myometria are listed in Table 1. The parameters that were determined in experiments with tissue from a specific mammal, either pregnant or nonpregnant, are detailed in Table 2. These parameters were extracted from the relevant literature. The remaining parameters were obtained by minimizing the mean square error (MSE) between the model output and the measured reaction of the cells under similar experimental conditions using the optimization toolbox provided in Matlab (version 7.1.4). An initial guess for the parameter values was made based on similar models, physiological constraints, or a low-resolution grid search within the boundaries of the parameter values. A local minimization process was then used to find parameter values that bring the MSE to a local minimum. The computed parameters are listed in Table 3. For process 3, the original relations between the rate constants that were assumed by Hai and Murphy (5) have been preserved. Accordingly, K₄ = K₃/4, K₁ = K₆, and K₅ = K₂. The stress is normalized by 0.8, since the result of the first constraint is that the maximal number of attached cross bridges equals 80% of the total amount of myosin.

Solution of the Model Equations

The model equations were numerically solved using an explicit Runge-Kutta (4, 5) formula, the Dormand-Prince pair, which is implemented in the Matlab (version 7.1.4) environment. For simulations of changes in C_Ca,i in response to various known membrane voltages, Eq. 6 is iteratively solved. At each iteration, Ca²⁺ fluxes described by Eqs. 1b, 4, and 5a are calculated according to the instantaneous C_Ca,i and membrane voltage. These values are used in Eq. 6 to find the instantaneous derivative of C_Ca,i. For complex types of simulations, such as a train of pulses, the time step was limited to a lower bound of 10⁻² s.

For simulations of the stress produced by muscle contraction in response to a known C_Ca,i, which can be either measured or simulated, Eqs. 7–9 are iteratively solved as follows. At each iteration, K₁ is calculated using Eq. 8 and the instantaneous C_Ca,i. The calculated values of K₁ and K₆ are then used in Eqs. 7, which are simultaneously solved to find the instantaneous derivative of the cross-bridge populations. Finally, the relative myosin phosphorylation and stress are calculated using Eqs. 9.

RESULTS

Control of Intracellular Ca²⁺ Concentration

Simulations of C_Ca,i decay were performed to evaluate the parameters characterizing the operation of Ca²⁺ extraction mechanisms, Ca²⁺ pumps, and Na⁺/Ca²⁺ exchangers. The contribution of various mechanisms for extraction of Ca²⁺ from the cell to the decay of C_Ca,i has been studied in pregnant rat myometria by comparing the decay rate in control conditions to the decay rate reached when one of the mechanisms was inhibited (19). The rise of Ca²⁺ was elicited by a train of voltage pulses from a holding potential of −80 mV to a pulse potential of 0 mV. Simulations under the same conditions were performed, allowing one to determine the appropriate parameters for each mechanism independently from the other. The parameters characterizing Ca²⁺ extraction by Na⁺/Ca²⁺ exchangers (C_Na,i, G_Na/Ca, and n) were optimized using the measured data of C_Ca,i decay during inhibition of the Ca²⁺ pumps. The result of this simulation is shown in Fig. 2A. Then, the parameters characterizing the Ca²⁺ pumps (G_Na/Ca and V_pmax) were optimized using the measured data of C_Ca,i decay in control conditions, when both mechanisms were active. The result of this simulation is shown in Fig. 2B. The parameters used for these simulations (which yielded the best fit to the experimental data) are presented in Table 3. C_Ca,i values calculated using the model fit very well with the experimental data.

Simulations of the C_Ca,i rise and decay following voltage pulses having a duration of 200 ms, from a holding potential of −50 mV to various potentials, were performed to verify that the model correctly describes the cell's response to this basic type of stimulation. Two parameters were set according to the results of these simulations: the intracellular Na⁺ concentration (C_Na,i) and the maximal calcium conductance (g_m,Ca). The transmembrane potential of the cell affects the concentrations of all ions in the cell. Since the holding potential in these experiments is different from the holding potential used in the previous simulations, C_Na,i and hence the reversal potential of the Na⁺/Ca²⁺ exchangers are also expected to change. Therefore, the previously set value of this parameter was changed for these simulations. In contrast, the maximal fluxes through the Ca²⁺ extraction mechanisms (G_Na/Ca and V_pmax) are determined by the amount of exchangers and pumps in the cell membrane. The fluxes may vary in different species and hormonal states, but they are expected to be similar for cells of the same specie at the same stage of pregnancy, which are at the same hormonal state. Furthermore, the concentrations for half-activation of the Ca²⁺ extraction mechanisms (K_ph and K_Na/Ca) and the Hill coefficient of the pumps (n) are expected to show even less variation, since these are intrinsic properties of the mechanisms. Accordingly, the values of these parameters were not changed. Finally, by setting the value of the maximal Ca²⁺ conductance (g_m,Ca), it was possible to calculate the rise in C_Ca,i and match it to the experimental data reported by Shmigol et al. (18).

The parameters C_Na,i and g_m,Ca were set using the measured C_Ca,i following voltage pulses to 0, +10, and −10 mV as shown in Fig. 3. The simulated changes in C_Ca,i following the voltage pulses are very similar to the experimental data. The fit is better for the voltage pulse to +10 mV than for the pulse to −10 mV, where the simulated decay rate is lower, and the final C_Ca,i level higher, then the measured values of these quantities. The model's ability to predict the C_Ca,i following voltage pulses to +20 and −20 mV was tested by using the same parameters as in the previous simulation. These predictions are depicted in Fig. 4 compared with the experimental data of Shmigol et al. (18). The prediction for a voltage pulse to +20 mV is similar to the measured behavior, whereas the simulation for a pulse to −20 mV underestimates level of C_Ca,i at the early stage of the decay, and later on it overestimates the experimental data. The maximal difference in this simulation between measured and simulated C_Ca,i is ∼16%. Taking into account that the confidence intervals of experimentally measured C_Ca,i may reach 20% of the measured values (18), this error is reasonable and within what can be interpreted as an expected intersubject variability.

The model was further used to estimate the relative contribution to C_Ca,i decay of the different mechanisms responsible for Ca²⁺ extraction. The results for Ca²⁺ decay at a holding potential of −80 mV are depicted in Fig. 5A. At all stages of this simulation, the Ca²⁺ efflux through Ca²⁺ pumps is higher than the efflux through Na⁺/Ca²⁺ exchangers. The contribution of all Ca²⁺ control mechanisms during Ca²⁺ rise and decay initiated by a voltage pulse to 0 mV from a holding potential of −50 mV is shown in Fig. 5B. At all stages of this simulation, the Ca²⁺ efflux through the Na⁺/Ca²⁺ exchangers is higher than the Ca²⁺ efflux through the pumps. The contribution of the Na⁺/Ca²⁺ exchangers relative to that of the membranal Ca²⁺ pumps to Ca²⁺ extraction from the cell is shown in Table 4. For a holding potential of −50 mV, which is also the resting potential for pregnant rat myometria, the contribution of the Na⁺/Ca²⁺ exchangers to Ca²⁺ extraction is significantly more dominant than the contribution of the Ca²⁺ pumps. When the holding potential is lower, the contribution of the Na⁺/Ca²⁺ exchangers declines, and, for a holding potential of −80 mV, it is less dominant than the contribution of the Ca²⁺ pumps. The model's ability to predict changes in C_Ca,i was further tested in response to a more complex stimulation of repetitive depolarization that was used in experiments (18). The predicted changes in C_Ca,i during stimulation of the cell by a voltage pulse train, consisting of 10 voltage pulses from a holding potential of −80 mV, having a pulse interval of 100 ms and an interpulse interval of 330 ms, closely match the experimental data and are depicted in Fig. 6. During the first voltage pulses, the level of C_Ca,i is slightly overestimated by the model, but the match between simulated C_Ca,i and measured C_Ca,i improves during the last five voltage pulses and the C_Ca,i decay.

Sensitivity Analysis

A sensitivity analysis was performed to find which parameters have the most significant affect on the characteristics of the cell reaction to a voltage pulse of 200 ms from a holding potential of −50 to 0 mV. Each parameter was changed within the reasonable physiological limits for both pregnant and nonpregnant myometrial cells. The results are presented in Table 5. The most significant changes in maximal and final (at t = 3 s) C_Ca,i resulted from changes in the cell volume. When the volume of the cell was changed to its value in the nonpregnant uterus, 5 × 10⁻¹² liters (34), a volume slightly higher than half the volume previously used for simulating C_Ca,i changes in SMCs from pregnant myometria, the maximal C_Ca,i increased by a factor of 4 and the final C_Ca,i by a factor of 4.5. These values are higher than the physiological range, indicating that in the nonpregnant cell the change of volume must be accompanied by a significant reduction in the Ca²⁺ channels conductivity or increased Ca²⁺ buffering by cytoplasmic proteins, reducing the increase in C_Ca,i. An increase in V_Ca,1/2 may also contribute to a reduction in C_Ca,i. Experimental measurements are required to find which of these explanations is in line with the physiological differences between pregnant and nonpregnant myometria.

The analysis also showed that changes in K_Ca,1/2 significantly affect the final C_Ca,i without significant changes in maximal C_Ca,i. This is caused by a change in the Ca²⁺ flux through the channels when the cell is held at the constant holding potential, due to the change in the form of the VOCC activation function (Eq. 3b), which indicates the fractional amount of open VOCCs, as shown in Fig. 7. A decrease in K_Ca,1/2 causes the channels to close at the holding potential of −50 mV, thereby decreasing C_Ca,i in the simulation, whereas an increase in K_Ca,1/2 causes an increase in the flux through the channels at the holding potential of −50 mV. Changes in the maximal possible fluxes through the Ca²⁺ pumps (V_pmax) and Na⁺/Ca²⁺ exchangers (G_Na/Ca) caused significant changes in the ratio of the fluxes through the exchangers and the pumps but less significant changes in C_Ca,i. Since the exchangers are more dominant than the pumps in Ca²⁺ extraction in this simulation, the change in the maximal flux through the exchangers (G_Na/Ca) had a more significant effect on C_Ca,i than the change in the maximal flux through the pumps (V_pmax).

The intracellular sodium concentration (C_Na,i) has a significant effect on C_Ca,i, since changes in C_Na,i change the reversal potential of the Na⁺/Ca²⁺ exchangers and thus may invert the direction of the Ca²⁺ flux through the exchangers as well as change its magnitude. The effect of changes in C_Na,i on the C_Ca,i response to a voltage pulse to 0 mV is shown in Fig. 8A. When C_Na,i increases, the Nernst potential for sodium decreases and thus also the reversal potential of the exchangers. When the reversal potential drops below the membrane potential, the Ca²⁺ efflux is inverted and Ca²⁺ begins to flow through the exchangers into the cell. For C_Na,i higher than 20 mM, the influx of Ca²⁺ through the exchangers is higher than the Ca²⁺ efflux through the pumps; thus, there is no C_Ca,i decay. The low affinity of the Na⁺/Ca²⁺ exchangers to C_Na,i is not considered in Eq. 5, which describes the Ca²⁺ flux through the exchangers. Therefore, it is possible that Ca²⁺ inflow into the cell, which is simulated by the model at high C_Na,i, will not occur in the real myometrial SMCs at the same conditions. Accordingly, the results of simulations with C_Na,i > 20 mM should be viewed with limited confidence.

The effects of changes of C_Na,i and C_Ca,i on J_Na/Ca when the membrane voltage is −50 mV are shown in Fig. 8B. The flux turns positive at a C_Na,i of ∼20 mM. When C_Ca,i increases, the changes in J_Na/Ca as a function of C_Na,i become more steep; accordingly, the mechanism is most sensitive to C_Na,i when Ca²⁺ must be extracted from the cell to return to resting conditions. To this date, only C_Ca,i changes during contractions have been measured in vitro. To validate the suggested effects of C_Na,i on establishing C_Ca,i changes, experimental simultaneous measurements of both are required.

None of the parameters had a significant affect on the time to peak in C_Ca,i; accordingly, it can be inferred that this time is set almost solely by the duration of the voltage pulse applied to the cell.

MLC Phosphorylation and Stress Production by the Contracting Cell

Using simultaneous measurements of C_Ca,i, MLC phosphorylation and force in human myometrium triggered to contract by electrical field stimulation (30) as well as force development and relaxation during a stretch-induced contraction, the parameters characterizing the third process described by the model, were set. The relative amount of phosphorylated MLC in reaction to a rise in C_Ca,i was calculated and compared with the experimental data, as shown in Fig. 9A. The simulated MLC phosphorylation shows an initial rapid increase at a similar rate to the measured initial phosphorylation increase rate. This is followed by an immediate stabilization to a plateau, whereas the measured response showed more oscillations. However, the final phosphorylation values are very similar. In the same experiment, the normalized force was also calculated and compared with the experimental data, as shown in Fig. 9B. A very close match was reached between the measured and simulated relative force produced by the muscle. The average difference between measured and simulated force in this simulation is ∼5%, whereas the difference in phosphorylation may reach 20%. However, when viewing these results, a variability of 10% in experimental measurements of both force and phosphorylation should be taken into account (30).

The simulated and measured relative force during stretch-induced contraction and relaxation are presented in Fig. 10. The parameters used for these simulations are listed in Table 3. The value of C_Ca,1/2MLCK was found to be 257 nM, within the range of 250–260 nM that was measured in the bovine tracheal SMC, when MLCK was allowed to be inhibited by CaM kinase 2 (23, 26). The Hill coefficient that was found to match the relation between C_Ca,i and MLCK activity (n_M = 8.76) is much higher then the coefficients that were found for the same process in the bovine trachealis SMC, 1.5–2.7 (25, 26). This indicates a steeper relation between Ca²⁺ and MLCK activation in myometrial myocytes. Accordingly, MLCK activity in the myometrial myocyte is very sensitive to minor changes in C_Ca,i near the half-activation concentration of 257 nM and reaches saturation at lower C_Ca,i compared with other SMCs. The rate of dephosphorylated cross-bridge detachment found (0.037) was similar to the rate of phosphorylated cross-bridge detachment (0.035), indicating that latch bridges may not have a significant role in force maintenance in the myometrial SMC.

The maximal measured and simulated force values are similar, although the measured force development slightly precedes the simulated force. During the relaxation phase, the measured and simulated forces show a similar decay rate. The significant difference in time scale between this stretch-induced contraction and the previous electrically induced contraction should be taken into account when viewing these results.

Simulation of Ca²⁺ Control and Stress Production

A complete solution of the entire model was used to simulate stress production in response to a calculated elevation of C_Ca,i due to membrane depolarization. The results for a voltage pulse of 1 s from a holding potential of −80 to 0 mV are depicted in Fig. 11A. The results for a series of pulses of a combined duration of 1 s (i.e., 10 pulses of 0.1 s on followed by 0.33 off) from a holding potential of −80 mV are shown in Fig. 11B. The developed stress was calculated for 20 s of contraction and relaxation. The maximal force produced in reaction to the voltage-pulse train (32.1%) was higher than the maximal force produced in reaction to the single voltage pulse (23.8%), although the single voltage pulse produced a significantly higher rise in C_Ca,i. This suggests that the rate of the rise in C_Ca,i has a more significant affect on force production than its extent. The delay between the time at which maximal C_Ca,i and maximal force were reached was 3.26 s for the pulse train and 4.63 s for the single pulse. It can thus also be inferred that repetitive spiking is more efficient for force production than prolonged depolarization. Accordingly, repetitive firing of action potentials that has been recorded in pregnant myometria, and is sometimes superimposed on a sustained depolarization (13), may produce force more effectively than single spike firing, even if the single spikes are of longer duration.

We have further tested the model ability to predict myometrial SMC behavior by simulating the response of a human pregnant myometrial cell to a plateau-type action potential. The changes in C_Ca,i and stress production in response to the measured membrane potential changes during the action potential firing (14) were computed. The results are shown in Fig. 12, where the simulated stress production is compared with the measured stress. All stresses were normalized to the maximal stress measured during this contraction. Despite differences between the properties of the human pregnant myometrium, from which the data was recorded, and the pregnant rat and nonpregnant human myometria, for which the model parameters were set, the model reproduces the measured behavior fairly accurately.

DISCUSSION

A model describing the relations between the transmembrane voltage, intracellular Ca²⁺ concentration, and force produced by a single myometrial SMC has been developed. We demonstrated that the model successfully reproduces experimental results, even though it accounts only for the most dominant mechanisms of C_Ca,i control in the cell. Thus, it can be further used to predict the contribution of different mechanisms to the cell response to changes in its transmembrane voltage as well as results of future experiments.

The present simulations (Fig. 5) show that C_Na,i is an important factor in setting C_Ca,i, although not much attention has been devoted to measurement of C_Na,i in myometrial SMCs. In addition, we have shown (Fig. 11) that repetitive spiking may enable force production more efficiently, with smaller delays and lower C_Ca,i compared with single spikes. The low conductivity of VOCCs that was found to fit experimental measurements indicates that significant buffering occurs in the myometrial cells. Experimental data are required to validate the ratio of 150 that was found between the amount of Ca²⁺ entering the cell and the amount Ca²⁺ contributing to the increase in C_Ca,i.

The parameters for the Na⁺/Ca²⁺ exchangers and for the Ca²⁺ pumps were separately determined; hence, a correct representation of the relative contribution of each of these mechanisms to the total extraction of Ca²⁺ from the cell was ensured. In simulations of C_Ca,i decay the model allows identification of the dominant mechanism responsible for Ca²⁺ extraction. There is an inconsistency in the literature regarding this matter. While some authors claim that the contribution of the Ca²⁺ pumps to Ca²⁺ extraction is lower than that of the Na⁺/Ca²⁺ exchangers (19), others claim that the Ca²⁺ pumps are more dominant (10). The results of the present study (Fig. 5 and Table 4) show that the relative contribution of the mechanisms depends on the level of C_Na,i, which, in turn, is set by the membrane voltage.

In simulations of C_Ca,i decay at two distinct holding potentials (Figs. 2 and 3), two different levels of C_Na,i were used to match experimental data (Table 3). Thus, the characteristic reversal potential of the Na⁺/Ca²⁺ exchangers was also different, which, in turn, changed the relative contribution of the Na⁺/Ca²⁺ exchangers to Ca²⁺ extraction from the cell. The value of C_Na,i = 16.55 mM, which best described the response of a cell held at −80 mV (19), was significantly higher than C_Na,i = 2.9836 mM, which best described the response of a cell held at −50 mV (18). The decrease of C_Na,i when the transmembrane voltage increases causes an increase in the Nernst potential for Na⁺, and, thus, the reversal potential of the exchangers increases as well (Eqs. 5b and 5c). As a result of the increased difference between the membrane voltage and the reversal potential, the Ca²⁺ flux through the exchangers (as given by Eq. 5a) increases when the holding potential increases. Therefore, simultaneous measurement of C_Na,i and C_Ca,i in voltage-clamp experiments, or an expansion of the model to include Na⁺ control mechanisms, may advance our understanding of how the level of C_Ca,i is controlled.

Further examination of the fluxes presented in Fig. 5B shows that the maximal flux through VOCC is ∼1,400 nM/s. Using the relation I = J·V_cell·Z_Ca·F, the maximal Ca²⁺ current through the channels is 0.00567 pA. The surface area of the myometrial SMC in the late pregnant rat is 7,600 μm² (34). Accordingly, the maximal simulated current density, given by the ratio of the total current and the surface area, is 0.0746 μA/cm². Since the measured density of I_Ca in the myometrial myocyte was found to vary between 3 and 11 μA/cm², where the higher values were measured in pregnant myometria, we find that the effect of buffering reduces the contribution of Ca²⁺ entry into the cell to C_Ca,i by a factor of ∼150, as expected in SMCs (15). Furthermore, the maximal Ca²⁺ conductance of all VOCCs in the cell was found to be 0.046842 nS. If buffering by a factor of ∼150 is taken into account, the total conductivity of all the channels is 7.03 nS. Dividing this by the conductance of a single L-type VOCC, 29 pS in the pregnant human myometrium (7), we find that ∼250 VOCCs can conduct current at the same time.

In simulations of voltage-clamp experiments C_Ca,ex is assumed to be constant, since this concentration is significantly higher than the change in C_Ca,i, because of Ca²⁺ entry into the cell. However, in myometrial tissue at in vivo conditions, the entry of Ca²⁺ into the cell may affect the local extracellular concentration of Ca²⁺. A local decrease in C_Ca,ex will lower the Nernst potential for calcium (E_Ca) and, thus, may limit Ca²⁺ entry through L-type VOCCs. In addition, the decreased E_Ca will cause an increase in the reversal potential of the Na⁺/Ca²⁺ exchangers, thereby increasing Ca²⁺ extraction by this mechanism. Accordingly, for simulations of myometrial SMCs in in vivo conditions, changes in C_Ca,ex should be accounted for.

An underestimation of C_Ca,i in simulations of the cell response to various stimulations might be expected because of exclusion of Ca²⁺ release from intracellular Ca²⁺ stores in the present model. However, its effects are probably compensated for by the overestimation of Ca²⁺ entry through the VOCCs as well as by neglecting the time-dependent changes in the conductivity of VOCCs and in C_Na,i. Shmigol et al. (18) demonstrated that in the response of a myometrial SMC to a series of voltage pulses, because of the contribution of CICR at the initial stages, the measured C_Ca,i increase is relatively constant during the first voltage pulses and drops at the last ones. Conversely, when only the VOCCs are active, the C_Ca,i increase is expected to drop with every consequent voltage pulse. In our simulations, the C_Ca,i increase due to voltage pulses is decreased with every consequent voltage pulse, because of the properties of the VOCCs. Accordingly, the simulated changes in C_Ca,i must be overestimated at the initial stages to coincide with the measured changes at the advanced stages of the simulation.

The Ca²⁺-dependent MLC phosphorylation and force production calculation was based on the model of Hai and Murphy (5), which was originally developed to describe the contraction of vascular SMCs. The MLC phosphorylation and force production by the myometrial SMC were successfully simulated using the model, despite the differences in properties of vascular and myometrial SMCs and the different time scales of simulated contractions (Figs. 9, 10, and 12).

The description of relaxation is limited by the assumption of a constant rate of myosin dephosphorylation by MLCP (K₂ and K₅). In addition, the rate constant for the detachment of latch bridges (K₇), which is computed by fitting to experimental data of Word et al. (30), was found to be similar to the rate constant for detachment of phosphorylated attached cross bridges (K₄). Since latch bridges have a low detachment rate, this indicates a limited role for latch bridges in the maintenance of the generated force by the myometrial SMCs.

The successful prediction of the cell response to basic types of simulations, such as voltage pulses and trains, cannot be extrapolated to prediction of the response to more complex types of stimulation, when the model assumptions are no longer valid. The model will probably underestimate the level of C_Ca,i when Ca²⁺ release from intracellular stores is significant. During prolonged stimulations, the level of C_Ca,i is expected to be slightly overestimated since the time-dependent inhibition of the L-type VOCCs is not taken into account. In addition, the estimation of a constant level of C_Na,i may also lead to an error in the calculation of Ca²⁺ extraction by Na⁺/Ca²⁺ exchangers. Accordingly, to describe the response to complex changes in membrane voltage, such as plateau potentials or bursts of spike potentials superimposed on plateau potentials, as observed in the pregnant rat myometrium (9), the Ca²⁺ control mechanisms must be described in more detail. This includes accounting for additional Ca²⁺ control mechanisms and for time-dependent changes in the properties of existing mechanisms. It will also be beneficial to include mechanisms of controlling additional ions concentrations, such as Na⁺ and K⁺ control mechanisms. Furthermore, to more accurately describe the time-dependent changes in the rates of the reactions involved in the process of contraction and relaxation, the dephosphorylation rate cannot be constant and should be modified. However, it has been shown that force production in response to a typical action potential recorded in the pregnant myometrium can be predicted using the model.

The parameters used for the simulations were optimized to reach the best agreement with experimental measurements by bringing the MSE between simulation and measurement to a local minimal value. It is possible that the global minimal value was not reached and that there are additional sets of parameters that can yield good predictions of the experimental data. From the existence of multiple parameters explaining the same cell response, it can be inferred that the cell can present a given reaction using various combinations of operation levels of its Ca²⁺ control and other mechanisms.

In conclusion, the mathematical model that has been developed in the present study successfully describes the basic processes of excitation and contraction, as well as relaxation of the myometrial SMC. The model can be utilized to study the operation of important cellular mechanisms of Ca²⁺ control, such as L-type VOCCs, Ca²⁺ pumps, and Na⁺/Ca²⁺ exchangers. In addition, it can be used to predict the stress that the muscle will produce, as it depends on the level of C_Ca,i. Since controlled experiments allow to inhibit the operation of specific mechanisms, while model simulations allow to both find many optional operation modes as well as exclude the operation of every possible combination of mechanisms, by combining the two it will be possible to improve the understanding of the operation of various mechanisms together, in the cell, at various conditions, to produce a variety of physiologically important cell behaviors.

FOOTNOTES

• The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.