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

1Introduction¶

The model Bursztyn et al., 2007 describes how the intracellular $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 $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 $Ca^{2+}$ control mechanisms: voltage-operated $Ca^{2+}$ channels, $Ca^{2+}$ pumps and $Na^{+}$/$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 $Ca^{2+}$ ions, including the flux through the L-type voltage dependent $Ca^{2+}$ channels $J_{VOCC}$, the efflux through $Ca^{2+}$ pump $J_{Ca,pump}$, and the flux through the $Na^{+}/Ca^{2+}$ exchangers $J_{Na/Ca}$. The computation of current through the L-type voltage dependent $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 $K_1$ and $K_6$ are regulated by the intracellular $Ca_{i}^{2+}$ concentration. Hence, the link between the excitation and contraction is established through the intracellular $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 $[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 $Ca^{2+}$ channel activation as a function of $V_m$;
4. the membrane potential simulation experiment is used to simulate the excitation and contraction response to a recorded human plateau potential;
5. the $Na_i^{+}$ experiment is used to evaluate the variation in the flux through the $Na^{+}/Ca^{2+}$ exchangers as a function of $[Na_i^{+}]$;
6. the $Ca_{i}^{2+}$ experiment 1 is used to simulate the reaction in human nonpregnant myometrium to an increase in $[Ca_{i}^{2+}]$; and
7. the $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 $Ca_{i}^{2+}$ Concentration¶

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

In the simulation of Figure 3 from the primary paper, the holding voltage is $-50$ mV and the corresponding $[Na_{i}]$ is $2.9836$ mM. At time $t=0$, the initial value of $[Ca_{i}^{2+}]$ is set to $130$ nM, while a $200$-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.

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 $Ca^{2+}$ flux through $Na^{+}/Ca^{2+}$ exchangers and $Ca^{2+}$ pumps during $[Ca_{i}^{2+}]$ decay at a holding voltage of $-80$ mV, while plot B shows $Ca^{2+}$ flux through $Ca^{2+}$ channels and $Ca^{2+}$ extraction mechanisms during $[Ca^{2+}]$ rise and decay in response to a $200$-ms voltage pulse to $0$ mV from holding potential of $-50$ mV.

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

3.2Sensitivity Analysis¶

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

Figure 6A shows the variation in $[Ca_{i}^{2+}]$ decay following a $200$-ms voltage pulse to $0$ mV due to changes in $K_{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 $[Ca_{i}^{2+}]$ rise and decay following a $200$-ms voltage pulse to $0$ mV due to changes in $[Na_i]$, while Figure 7B shows the variation in the flux through $Na^{+}/Ca^{2+}$ exchangers.

3.3MLC Phosphorylation and Stress Production by the Contracting Cell¶

We use a curve fitting method to generate a $[Ca_{i}^{2+}]$ profile close to the measured values. This $[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 $[Ca_{i}^{2+}]$. The model response to the piecewise linear function is shown in Figure 9.

3.4Simulation of $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$ mV to a pulse potential of $0$ mV, while Figure 10A shows the response to a train of 10 pulses from a holding potential of $-80$ mV to a pulse potential of $0$ mV, which has a duration of $100$ ms and the interval between pulses is $330$ 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$ mV.

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

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.

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

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