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

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 C a 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


Introduction
The model (Bursztyn et al., 2007) describes how the intracellular C a 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 C a 2+ i 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.

Model description 2.1 Primary Publication
The model (Bursztyn et al., 2007) incorporates three C a 2+ control mechanisms: voltage-operated C a 2+ channels, C a 2+ pumps and N a + /C a 2+ exchangers, which employ the mathematical formulation proposed in (Parthimos et al., 1999). The cross-bridge model of Hai and 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 and 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.

Modularized 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 C a 2+ ions, including the flux through the L-type voltage dependent C a 2+ channels J V O C C , the efflux through C a 2+ pump J C a,pump , and the flux through the N a + /C a 2+ exchangers J N a/C a . The computation of current through the L-type voltage dependent C a 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 and Murphy (1988), where the rate constants of MLC phosphorylation K 1 and K 6 are regulated by the intracellular C a 2+ i concentration. Hence, the link between the excitation and contraction is established through the intracellular C a 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 [C a 2+ i ] 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 C a 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 N a + i experiment is used to evaluate the variation in the flux through the N a + /C a 2+ exchangers as a function of [N a + i ]; 6. the C a 2+ i experiment 1 is used to simulate the reaction in human nonpregnant myometrium to an increase in [C a 2+ i ]; and 7. the C a 2+ i 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 and 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,

2/8
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).

Model 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.

Control of Intracellular C a 2+
i Concentration To simulate the intracellular C a 2+ concentration decay, we need to increase the [C a 2+ i ] 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 [C a 2+ i ]. In our simulation of Figure  2A, we set the [C a 2+ i ] value at time t = 0 to the experimentally measured one, i.e., 980.63 nM, and simulate the subsequent C a 2+ i decay due to extrusion of C a 2+ i by the N a + /C a 2+ exchangers only. Figure 2B describes a second simulation under similar conditions but with both C a 2+ i extrusion mechanisms (exchanges and pumps) active. The membrane potential throughout the simulation holds at −80 mV and the corresponding [N a i ] is 16.55 mM.
In the simulation of Figure 3 from the primary paper, the holding voltage is −50 mV and the corresponding [N a i ] is 2.9836 mM. At time t = 0, the initial value of [C a 2+ i ] 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 C a 2+ flux through N a + /C a 2+ exchangers and C a 2+ pumps during [C a 2+ i ] decay at a holding voltage of −80 mV, while plot B shows C a 2+ flux through C a 2+ channels and C a 2+ extraction mechanisms during [C a 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.

Sensitivity Analysis
The effect of parameter K C a,1/2 and intracelluar N a + concentration [N a + 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 ρ V C a in Figure 6B. While the linear [N a + i ] from 1 mM to 46 mM is created to compute the flux J N a/C a in Figure 7B, the membrane potential holds at −50 mV. Figure 6A shows the variation in [C a 2+ i ] decay following a 200-ms voltage pulse to 0 mV due to changes in K C a,1/2 , while the variation in the activation function in Figure 6B indicates the fractional amount of open VOCCs.     Fig. 6 in Bursztyn et al. (2007)). Figure 7A shows the variation in [C a 2+ i ] rise and decay following a 200-ms voltage pulse to 0 mV due to changes in [N a i ], while Figure 7B shows the variation in the flux through N a + /C a 2+ exchangers.

MLC Phosphorylation and Stress Production by the Contracting Cell
We use a curve fitting method to generate a [C a 2+ i ] profile close to the measured values. This [C a 2+ i ] is used to drive the cross-bridge model of Hai and 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 [C a 2+ i ]. The model response to the piecewise linear function is shown in Figure 9.

Simulation of C a 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 5/8 Figure 6. Sensitivity analysis for K C a,1/2 . A: variation in [C a 2+ i ] decay B: variation in the activation function (c.f., Fig. 7 in Bursztyn et al. (2007)).

Figure 11. Simulation of changes in [C a 2+
i ] and stress in response to a plateau potential V m (c.f., Fig. 12 in Bursztyn et al. (2007)).
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

Discussion
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 Parameters primary paper current CellML V p,max 5.1449e-4 mM · s −1 · mV −1 5.1449e-4 mM · s −1 G N a/C a 5.7297e-5 mM · s −1 · mV 5.7297e-5 mM · s −1 · mV 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:

Acknowledgement
The Freifeld lab is supported by the Zuckerman STEM Leadership Program and the Taub family foundation.
There are details/documentation on how the source code was compiled There are details on how to run the code in the provided documentation The initial conditions are provided for each of the simulations Details for creating reported graphical results from the simulation results Source code: a declarative language is used (e.g. SBML, CellML, NeuroML) The algorithms used are defined or cited in previous articles The algorithm parameters are defined Post-processing of the results are described in sufficient detail

Executable model provided:
The model is executable without source (e.g. desktop application, compiled code, online service) There are sufficient details to repeat the required simulation experiments

The model is described mathematically in the article(s):
Equations representing the biological system There are tables or lists of parameter values There are tables or lists of initial conditions

Machine-readable tables of initial conditions
The simulation experiments using the model are described mathematically in the article: Integration algorithms used are defined Stochastic algorithms used are defined Random number generator algorithms used are defined Parameter fitting algorithms are defined The paper indicates how the algorithms yield the desired output