Mathematical model of excitation-contraction in a uterine smooth muscle cell

Abstract

The Bursztyn et al. (2007) paper proposes a mathematical model of excitation-contraction in a myometrial smooth muscle cell (SMC). The model incorporates processes of intracellular Ca2+Ca^{2+} concentration control, myosin light chain (MLC) phosphorylation and stress production. We create a modularized CellML implementation of the model, which is able to simulate these processes against the original data.

Keywords:cellular calcium control mechanismsmyometrial contractionscomputational modelCellML

1Introduction

The model Bursztyn et al., 2007 describes how the intracellular Ca2+Ca^{2+} concentration increases due to membrane depolarization, which leads to an increase in the myosin phosphorylation rate and therefore the generation of muscle contraction. Here, we divide the mathematical model into distinct sub-modules encoded in CellML enabling reuse of the various sub-modules in future studies and models. We also reuse a few previously developed sub-modules at https://models.physiomeproject.org/workspace/6bc, such as Nernst potential computation and stimulation protocols, to reduce the implementation workload and potentially improve the model quality. We successfully reproduce control of intracellular Cai2+Ca_{i}^{2+} concentration, MLC phosphorylation, stress production and sensitivity analysis. From the primary paper we extracted data using the Engauge digitizer software Mitchell et al., 2020 to compare the current simulation results against those in the primary publication.

2Model description

2.1Primary Publication

The model Bursztyn et al., 2007 incorporates three Ca2+Ca^{2+} control mechanisms: voltage-operated Ca2+Ca^{2+} channels, Ca2+Ca^{2+} pumps and Na+Na^{+}/Ca2+Ca^{2+} exchangers, which employ the mathematical formulation proposed in Parthimos et al. (1999). The cross-bridge model of Hai & Murphy (1988) is used to describe the processes of myosin light chain (MLC) phosphorylation and stress production, which is essentially a deterministic multi-state Markov model (MM) Sakmann & Neher, 1996. The primary paper Bursztyn et al., 2007 has shown that the model is able to reproduce the results against experimental data obtained in pregnant rats and measurements in human nonpregnant myometrial cells. There is no publicly available code for this model, making it difficult to reuse this work. We therefore present in this article a CellML implementation for researchers to reuse in further developments based on this model.

2.2Modularized CellML model

The modularized CellML implementation is available in the Physiome Model Repository (PMR) at https://models.physiomeproject.org/workspace/6bb and the documentation can be found in the corresponding exposure at https://models.physiomeproject.org/e/742. In this manuscript we focus on reproducibility and reusability. The main components of this model and the performed simulation experiments are summarized as follows.

The J_Ca component defines the fluxes of Ca2+Ca^{2+} ions, including the flux through the L-type voltage dependent Ca2+Ca^{2+} channels JVOCCJ_{VOCC}, the efflux through Ca2+Ca^{2+} pump JCa,pumpJ_{Ca,pump}, and the flux through the Na+/Ca2+Na^{+}/Ca^{2+} exchangers JNa/CaJ_{Na/Ca}. The computation of current through the L-type voltage dependent Ca2+Ca^{2+} channels reuses the imported ionic current component ( https://models.physiomeproject.org/workspace/6bc).

The EC_uSMC component combines the excitation description J_Ca component, the reversal potential defined in the Nernst potential component, and the imported CB4HM component ( https://models.physiomeproject.org/workspace/6bc). The CB4HM component defines the four-state cross-bridge model of Hai & Murphy (1988), where the rate constants of MLC phosphorylation K1K_1 and K6K_6 are regulated by the intracellular Cai2+Ca_{i}^{2+} concentration. Hence, the link between the excitation and contraction is established through the intracellular Ca2+Ca^{2+} concentration. The component can simulate MLC phosphorylation and stress production process provided various voltage stimuli.

To test individual processes, we built the Unit_uSMC component, which decouples the connection between the excitation and contraction, thus enabling the use of an arbitrary [Cai2+][Ca_{i}^{2+}] profile to validate the response of phosphorylation and stress development.

Following best practices, our CellML implementation separates the mathematics from the parameterisation of the model. The mathematical model is imported into a specific parameterised instance in order to perform numerical simulations. The parameterisation would include defining the stimulus protocol to be applied. We have seven sets of simulation experiments and corresponding simulation results to reproduce the corresponding figures in the primary publication:

  1. the single stimulation experiment is used to validate the excitation and contraction response to a single voltage pulse with various durations and levels;
  2. the multiple stimulation experiment is used compare the reaction to the voltage-pulse train with a single pulse;
  3. the voltage ramp experiment is used to compute the Ca2+Ca^{2+} channel activation as a function of VmV_m;
  4. the membrane potential simulation experiment is used to simulate the excitation and contraction response to a recorded human plateau potential;
  5. the Nai+Na_i^{+} experiment is used to evaluate the variation in the flux through the Na+/Ca2+Na^{+}/Ca^{2+} exchangers as a function of [Nai+][Na_i^{+}];
  6. the Cai2+Ca_{i}^{2+} experiment 1 is used to simulate the reaction in human nonpregnant myometrium to an increase in [Cai2+][Ca_{i}^{2+}]; and
  7. the Cai2+Ca_{i}^{2+} experiment 2 is used to simulate the reaction during stretch-induced phasic contraction of human myometrium.

Simulation settings and detailed solver information are encoded in SED-ML documents for execution of the simulation experiments Waltemath et al., 2011. The Python scripts used with OpenCOR Snapshot 2021-05-19 Garny & Hunter, 2015 to perform simulations and produce the figures in the primary paper are also included in the folder ./Simulation/src. The name of the simulation and plot scripts indicates the figure number in the primary paper. For example, Fig2_sim.py is used to generate the simulation data and Fig2_plot.py reproduces the graph shown in Figure 2 from Bursztyn et al. (2007).

3Model results

In the presented figures, the dots denote the simulated data extracted from the primary publication, while the solid lines are the simulation results produced by the current CellML implementation.

3.1Control of Intracellular Cai2+Ca_{i}^{2+} Concentration

To simulate the intracellular Ca2+Ca^{2+} concentration decay, we need to increase the [Cai2+][Ca_{i}^{2+}] first. In the experimental paper Shmigol et al., 1998, they used a train of pulses from a holding potential of 80-80 mV to a step potential of 00 mV to generate an increase in [Cai2+][Ca_{i}^{2+}]. In our simulation of Figure 2A, we set the [Cai2+][Ca_{i}^{2+}] value at time t=0t=0 to the experimentally measured one, i.e., 980.63980.63 nM, and simulate the subsequent Cai2+Ca_{i}^{2+} decay due to extrusion of Cai2+Ca_{i}^{2+} by the Na+/Ca2+Na^{+}/Ca^{2+} exchangers only. Figure 2B describes a second simulation under similar conditions but with both Cai2+Ca_{i}^{2+} extrusion mechanisms (exchanges and pumps) active. The membrane potential throughout the simulation holds at 80-80 mV and the corresponding [Nai][Na_{i}] is 16.5516.55 mM.

In the simulation of Figure 3 from the primary paper, the holding voltage is 50-50 mV and the corresponding [Nai][Na_{i}] is 2.98362.9836 mM. At time t=0t=0, the initial value of [Cai2+][Ca_{i}^{2+}] is set to 130130 nM, while a 200200-ms pulse voltage is applied. The potential levels vary from plot A to C. Figure 4 uses the similar settings with different pulse potentials as well.

The experiment setting is summarized in Table 1 and the simulation results are shown in Figures 1, 2 and 3.

Table 1:Experiment setting

FigTest voltage (mV)Pulse duration (ms)[Nai][Na_i] (mM)Initial [Cai][Ca_i] (nM)
2ANANA16.55980.63
2BNANA16.55980.63
3A02002.9836130
3B102002.9836130
3C-102002.9836130
4A-202002.9836130
4B202002.9836130
Simulations of [Ca_{i}^{2+}] decay (in nM). A: under inhibition of Ca^{2+} pumps. B: in control conditions (c.f., Fig. 2 in ).

Figure 1:Simulations of [Cai2+][Ca_{i}^{2+}] decay (in nM). A: under inhibition of Ca2+Ca^{2+} pumps. B: in control conditions (c.f., Fig. 2 in Bursztyn et al. (2007)).

Simulation of [Ca_{i}^{2+}] rise and decay following a 200-ms voltage pulse of 0 mV (A), 10 mV (B), and -10 mV (C) (c.f., Fig. 3 in ).

Figure 2:Simulation of [Cai2+][Ca_{i}^{2+}] rise and decay following a 200200-ms voltage pulse of 00 mV (A), 1010 mV (B), and 10-10 mV (C) (c.f., Fig. 3 in Bursztyn et al. (2007)).

The experiment setting to reproduce primary Figure 5A and 5B are the same as the one used in primary Figure 2B and Figure 3A, respectively. The simulation results are shown in Figure 4. Plot A shows Ca2+Ca^{2+} flux through Na+/Ca2+Na^{+}/Ca^{2+} exchangers and Ca2+Ca^{2+} pumps during [Cai2+][Ca_{i}^{2+}] decay at a holding voltage of 80-80 mV, while plot B shows Ca2+Ca^{2+} flux through Ca2+Ca^{2+} channels and Ca2+Ca^{2+} extraction mechanisms during [Ca2+][Ca^{2+}] rise and decay in response to a 200200-ms voltage pulse to 00 mV from holding potential of 50-50 mV.

Simulation of [Ca_{i}^{2+}] rise and decay following a 200-ms voltage pulse of -20 mV (A), 20 mV (B) (c.f., Fig. 4 in ).

Figure 3:Simulation of [Cai2+][Ca_{i}^{2+}] rise and decay following a 200200-ms voltage pulse of 20-20 mV (A), 2020 mV (B) (c.f., Fig. 4 in Bursztyn et al. (2007)).

Simulations of Ca^{2+} fluxes through various Ca^{2+} control mechanisms, including Ca^{2+} entry and extraction from the cell (c.f., Fig. 5 in ).

Figure 4:Simulations of Ca2+Ca^{2+} fluxes through various Ca2+Ca^{2+} control mechanisms, including Ca2+Ca^{2+} entry and extraction from the cell (c.f., Fig. 5 in Bursztyn et al. (2007)).

Figure 5 shows the model reaction to a train of pulse potentials.

Simulation of [Ca_{i}^{2+}] rise and decay in response to a train of 10 voltage pulses of 100 ms duration from a holding potential of -80 mV to a pulse potential of 0 mV, with an interval of 330 ms between the pulses (c.f., Fig. 6 in ).

Figure 5:Simulation of [Cai2+][Ca_{i}^{2+}] rise and decay in response to a train of 10 voltage pulses of 100100 ms duration from a holding potential of 80-80 mV to a pulse potential of 00 mV, with an interval of 330330 ms between the pulses (c.f., Fig. 6 in Bursztyn et al. (2007)).

3.2Sensitivity Analysis

The effect of parameter KCa,1/2K_{Ca,1/2} and intracelluar Na+Na^{+} concentration [Nai+][Na_i^+] is shown in Figure 6 and Figure 7. Figure 6A and Figure 7A simulations are performed on a holding potential of 50-50 mV and the experiment setting is listed in Table 2. We provide a voltage ramp ranging from 100-100 mv to 6060 mV to compute the activation function ρVCa\rho_{VCa} in Figure 6B. While the linear [Nai+][Na_i^{+}] from 11 mM to 4646 mM is created to compute the flux JNa/CaJ_{Na/Ca} in Figure 7B, the membrane potential holds at 50-50 mV.

Figure 6A shows the variation in [Cai2+][Ca_{i}^{2+}] decay following a 200200-ms voltage pulse to 00 mV due to changes in KCa,1/2K_{Ca,1/2}, while the variation in the activation function in Figure 6B indicates the fractional amount of open VOCCs.

Figure 7A shows the variation in [Cai2+][Ca_{i}^{2+}] rise and decay following a 200200-ms voltage pulse to 00 mV due to changes in [Nai][Na_i], while Figure 7B shows the variation in the flux through Na+/Ca2+Na^{+}/Ca^{2+} exchangers.

Table 2:Experiment setting

FigTest voltage (mV)Pulse duration (ms)[Nai][Na_i] (mM)Initial [Cai][Ca_i] (nM)
7A02002.9836130
8A0200shown in the legend130
Sensitivity analysis for K_{Ca,1/2}. A: variation in [Ca_{i}^{2+}] decay B: variation in the activation function (c.f., Fig. 7 in ).

Figure 6:Sensitivity analysis for KCa,1/2K_{Ca,1/2}. A: variation in [Cai2+][Ca_{i}^{2+}] decay B: variation in the activation function (c.f., Fig. 7 in Bursztyn et al. (2007)).

Sensitivity analysis for [Na_{i}^{+}]. A: variation in [Ca_{i}^{2+}] rise and decay. B: variation in the flux through Na^{+}/Ca^{2+} exchangers (c.f., Fig. 8 in ).

Figure 7:Sensitivity analysis for [Nai+][Na_{i}^{+}]. A: variation in [Cai2+][Ca_{i}^{2+}] rise and decay. B: variation in the flux through Na+/Ca2+Na^{+}/Ca^{2+} exchangers (c.f., Fig. 8 in Bursztyn et al. (2007)).

3.3MLC Phosphorylation and Stress Production by the Contracting Cell

We use a curve fitting method to generate a [Cai2+][Ca_{i}^{2+}] profile close to the measured values. This [Cai2+][Ca_{i}^{2+}] is used to drive the cross-bridge model of Hai & Murphy (1988). The simulation result is shown in Figure 8. To simulate the stretch-induced contraction and relaxation, we use a piecewise linear function to construct the input [Cai2+][Ca_{i}^{2+}]. The model response to the piecewise linear function is shown in Figure 9.

Simulation of MLC phosphorylation (A) and force production (B) stress in response to an increase in [Ca_{i}^{2+}] (C) (c.f., Fig. 9 in ).

Figure 8:Simulation of MLC phosphorylation (A) and force production (B) stress in response to an increase in [Cai2+][Ca_{i}^{2+}] (C) (c.f., Fig. 9 in Bursztyn et al. (2007)).

Stress development and relaxation during stretch-induced phasic contraction of human myometrium (c.f., Fig. 10 in ).

Figure 9:Stress development and relaxation during stretch-induced phasic contraction of human myometrium (c.f., Fig. 10 in Bursztyn et al. (2007)).

3.4Simulation of Ca2+Ca^{2+} Control and Stress Production

The whole model simulation is performed by providing pacing voltage pulses mimicking the membrane depolarization. Figure 10A shows the response to a single 1-s voltage pulse from a holding potential of 80-80 mV to a pulse potential of 00 mV, while Figure 10A shows the response to a train of 10 pulses from a holding potential of 80-80 mV to a pulse potential of 00 mV, which has a duration of 100100 ms and the interval between pulses is 330330 ms.

We use a piecewise linear function to simulate a recorded human plateau potential of pregnant myometrium. Figure 11 shows the reaction of the cell model to this action potential following a holding potential of 50-50 mV.

The experiment setting for Figure 10 and Figure 11 is summarized in Table 3

Table 3:Experiment setting

FigTest Pulse[Nai][Na_i] (mM)Initial [Cai][Ca_i] (nM)
11Asingle 1-s pulse of 0 mV16.55130
11Bten 100-ms pulses of 0 mV with 300-ms interval16.55130
12simulated action potential2.983639
Simulations of changes in [Ca_{i}^{2+}] and stress in response to pacing pulses (c.f., Fig. 11 in ).

Figure 10:Simulations of changes in [Cai2+][Ca_{i}^{2+}] and stress in response to pacing pulses (c.f., Fig. 11 in Bursztyn et al. (2007)).

Simulation of changes in [Ca_{i}^{2+}] and stress in response to a plateau potential V_m (c.f., Fig. 12 in ).

Figure 11:Simulation of changes in [Cai2+][Ca_{i}^{2+}] and stress in response to a plateau potential VmV_m (c.f., Fig. 12 in Bursztyn et al. (2007)).

4Discussion

We implemented the model Bursztyn et al., 2007 using CellML in a modularized style which can be reused in future. We have successfully replicated all the figures and summarized detailed experiment settings for simulation in Section 3. In doing this, we noticed trivial typographical errors in parameter units and references in Table 3 of Bursztyn et al. (2007). Hence, we correct these in Table 4 to remove potential confusion.

Table 4:Correction of primary Table 3

Parametersprimary papercurrent CellML
Vp,maxV_{p,max}5.1449e-4 mMs1mV1mM\cdot s^{-1} \cdot mV^{-1}5.1449e-4 mMs1mM\cdot s^{-1}
GNa/CaG_{Na/Ca}5.7297e-5 mMs1mVmM\cdot s^{-1} \cdot mV5.7297e-5 mMs1mV1mM\cdot s^{-1} \cdot mV^{-1}
[Nai][Na_{i}]16.55 mMmM (18)2.9836 mMmM (18)
[Nai][Na_{i}]2.9836 mMmM (19)16.55 mMmM (19)

Another implementation note is that the primary paper used “efflux” to indicate the direction of Jca,pumpJ_{ca,pump} described in Equation (4) with a positive sign. In CellML implementation, the direction information is included in the equation. That is,

JCa,pump(t)=Vpmax[Cai(t)]n[Kph]n+[Cai(t)]nJ_{Ca,pump}(t)=-V_{pmax}\dfrac{[Cai(t)]^n}{[K_{ph}]^n+[Cai(t)]^n}

Consequently, CellML model uses a positive sign before the term Jca,pumpJ_{ca,pump} in the intracellular calcium concentration equation:

dCai(t)dt=JVOCC(t)+JCa,pump(t)+JNa/Ca(t)\frac{\mathrm{d}Cai(t) }{\mathrm{d} t}=J_{VOCC}(t)+J_{Ca,pump}(t)+J_{Na/Ca}(t)

Acknowledgments

The Freifeld lab is supported by the Zuckerman STEM Leadership Program and the Taub family foundation.

References
  1. Bursztyn, L., Eytan, O., Jaffa, A. J., & Elad, D. (2007). Mathematical model of excitation-contraction in a uterine smooth muscle cell. American Journal of Physiology-Cell Physiology, 292(5), C1816–C1829.
  2. Mitchell, M., Muftakhidinov, B., Winchen, T., Wilms, A., Schaik, B. van, badshah400, Mo-Gul, Badger, T. G., Jędrzejewski-Szmek, Z., kensington, & kylesower. (2020). markummitchell/engauge-digitizer: Nonrelease. Zenodo. 10.5281/zenodo.3941227
  3. Parthimos, D., Edwards, D. H., & Griffith, T. (1999). Minimal model of arterial chaos generated by coupled intracellular and membrane Ca2+ oscillators. American Journal of Physiology-Heart and Circulatory Physiology, 277(3), H1119–H1144.
  4. Hai, C.-M., & Murphy, R. A. (1988). Cross-bridge phosphorylation and regulation of latch state in smooth muscle. American Journal of Physiology-Cell Physiology, 254(1), C99–C106.
  5. Sakmann, B., & Neher, E. (1996). Single-channel recording. General Pharmacology: The Vascular System, 27(6), 1078. https://doi.org/10.1016/0306-3623(96)90073-7
Funding and support by:
International Union of Physiological SciencesAuckland Bioengineering InstituteDigital ScienceVPH Institute