Steady-state approximations for Hodgkin-Huxley cell models: Reduction of order for uterine smooth muscle cell model

Shawn A. Means; Mathias W. Roesler; Amy S. Garrett; Leo Cheng; Alys R. Clark

doi:10.1371/journal.pcbi.1011359

Abstract

Multi-scale mathematical bioelectrical models of organs such as the uterus, stomach or heart present challenges both for accuracy and computational tractability. These multi-scale models are typically founded on models of biological cells derived from the classic Hodkgin-Huxley (HH) formalism. Ion channel behaviour is tracked with dynamical variables representing activation or inactivation of currents that relax to steady-state dependencies on cellular membrane voltage. Timescales for relaxation may be orders of magnitude faster than companion ion channel variables or phenomena of physiological interest for the entire cell (such as bursting sequences of action potentials) or the entire organ (such as electromechanical coordination). Exploiting these time scales with steady-state approximations for relatively fast-acting systems is a well-known but often overlooked approach as evidenced by recent published models. We thus investigate feasibility of an extensive reduction of order for an HH-type cell model with steady-state approximations to the full dynamical activation and inactivation ion channel variables. Our effort utilises a published comprehensive uterine smooth muscle cell model that encompasses 19 ordinary differential equations and 105 formulations overall. The numerous ion channel submodels in the published model exhibit relaxation times ranging from order 10⁻¹ to 10⁵ milliseconds. Substitution of the faster dynamic variables with steady-state formulations demonstrates both an accurate reproduction of the full model and substantial improvements in time-to-solve, for test cases performed. Our demonstration here of an effective and relatively straightforward reduction method underlines the particular importance of considering time scales for model simplification before embarking on large-scale computations or parameter sweeps. As a preliminary complement to more intensive reduction of order methods such as parameter sensitivity and bifurcation analysis, this approach can rapidly and accurately improve computational tractability for challenging multi-scale organ modelling efforts.

Author summary

Mathematical modeling of physiological organ systems encompassing intracellular to organ-wide behaviour grapple with intrinsically multi-scaled systems in both space and time. Significant computational challenges arise with their numerical solution, often substantially constraining feasibility of in silico investigations. Naturally, streamlining and reducing the mathematical complexity of these models while maintaining accurate reproduction of experimental data is a persistent concern. We present a straightforward method for exploiting multiple time scales reducing the order of models by way of steady-state approximations at the cellular level. Applied to a uterine smooth muscle cell model, we obtain substantial improvements in time to solve by around a factor of two while accurately reproducing full model results. Our successful demonstration here highlights the importance and relative ease of the method—that is not altogether unknown but often overlooked—with potential application across a wide variety of multi-scale organ modeling efforts.

Citation: Means SA, Roesler MW, Garrett AS, Cheng L, Clark AR (2023) Steady-state approximations for Hodgkin-Huxley cell models: Reduction of order for uterine smooth muscle cell model. PLoS Comput Biol 19(8): e1011359. https://doi.org/10.1371/journal.pcbi.1011359

Editor: Andrew D. McCulloch, University of California San Diego, UNITED STATES

Received: December 21, 2022; Accepted: July 14, 2023; Published: August 30, 2023

Copyright: © 2023 Means et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Both the Full Tong, 2011 and Reduced Tong models analysed here are available at the Physiome Repository via the following links: https://models.physiomeproject.org/workspace/261 and https://models.physiomeproject.org/workspace/8bc, respectively. They are ready to run with the OpenCOR environment available for free download at https://opencor.ws/.

Funding: This work was funded, in part, by grants from the Ministry of Business Innovation and Employment’s Catalyst: Strategic fund. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. There was no additional external funding received for this study. See https://www.mbie.govt.nz/science-and-technology/science-and-innovation/funding-information-and-opportunities/investment-funds/catalyst-fund/catalyst-strategic-auckland-bioengineering-institute-12-labours-project/ for more.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Multi-scale mathematical models simulating electrical activity of muscular organs, such as the cardiac [1], gastro-intestinal [2] and uterine [3] systems, rely heavily on accurate descriptions of cellular functions. Striking a balance between realistic representations at the cell level and computational feasibility at the organ scale is a critical consideration in multi-scale modelling. Most biophysical cell models represent complex arrays of ion currents and transporters for comprehensive descriptions of their cellular functions [4]. However, scaling to the organ level often compels simplification of cell level models, with a challenge of retaining key functional contributors in a model without excessive complexity and computational intractability.

A classical example of model simplification is the FitzHugh-Nagumo (FN) reduced order version of the detailed Hodgkin-Huxley (HH) giant squid axon action potential model [5, 6]. The HH model exhibits multiple and wide-ranging time scales and FitzHugh exploited these time scale differences for the FN reduction. Other formal reduction of order methods such as sensitivity and bifurcation analysis are powerful and effective techniques, but often involve considerable effort. Typically encompassing high dimensional parameter spaces, local or global sensitivity analysis methods such as the Sobol or partial rank correlation coefficient [7] combined with dynamic reduction such as that presented by Parthimos, et al. [8] are often non-trivial.

Several multi-scale mathematical studies investigated organ-wide uterine function built on individual uterine smooth muscle cell (uSMC) modeling efforts—as recently reviewed by our group [9]. These cell models range from detailed and complex arrays of numerous ion currents and transporters [10–12] to minimal methods using the classic FN model [13]. Recent trends are toward comprehensive cell level models, but this can significantly complicate assembly and numerical solution when extended to full multi-scale organ models. Particularly, investigations of longer time scale phenomena such as uterine contractions during labour as compared to rapid time scales of cellular behaviour suggest the potential for simplification.

Ion channel submodels in uSMC models exhibit relaxation time scales ranging from less than ten milliseconds to hundreds of seconds, yet organ scale behaviour stretches out to orders of seconds and longer. In this study, we follow an approach inspired by that of FitzHugh [5] for simplifying a complex uSMC model into a more readily solvable system. Full dynamic variables with relatively fast relaxation times are exploited and replaced with steady-state approximations, greatly mitigating computational expense. This method is demonstrated on a comprehensive uSMC model published by Tong et al. [12]. Given the prevalence of HH-type cell models, however, this approach may be applied to a wide range of cell types such as gastro-intestinal [14], cardiac [15], or of course neuronal [5]. Recent multi-scale models for these systems could moreover benefit from consideration of time-scales for these easily overlooked steady-state approximations and corresponding improvements to computational expense [1, 2].

Methods

Reduction of order methodology

We describe here our approach to simplifying a HH-type ion channel cell model based on the extensive effort of Tong, et al., both the initial 2011 and later updated 2014 versions [12, 16]. The complete published version of Tong, et al., 2011 we denote as the ‘Full Tong Model’ (FTM), and our reduced version as simply the ‘Reduced Tong Model’ (RTM). The schematic in Fig 1 shows all components of the FTM model, and the mechanisms that are ultimately removed in our analysis for the RTM. Avoiding exhaustive sensitivity and bifurcation analyses, the approach outlined here can inform other multi-scale organ modeling efforts, particularly for ion channel models exhibiting multiple time scales. The procedure utilised is two-fold: (i) eliminate channels in the cell model based on conductance. This occurs where FTM mechanisms are easily substituted by other mechanisms with modest modifications while maintaining qualitative reproduction of the full model; and (ii) timescale based reduction, in which we examine and replace the ion channel submodel dynamic variables remaining after the reduction of (i) exhibiting relatively fast time scales with their corresponding steady-state approximations. Note, we do not provide all formulations and expressions out of the full suite of 105 equations of the FTM; see Tong, et al. 2011 for details [12]. As the original FTM performed qualitative reproductions of experimental uSMC data; we do not overall apply quantitative metrics for comparing each output from our reduction effort. However, we do present both model results with relative differences upon completion of building the complete RTM.

Download:

Fig 1. Schematic of Full Tong uSMC Model (FTM).

All ion channel, exchanger and pump mechanisms of the 2011 model are presented, with the later added mechanisms of the 2014 version designated as shown. Note, only Ca²⁺ion concentrations dynamically vary in the Tong models; Na⁺, K⁺ and Cl⁻ are held fixed throughout. Ca²⁺concentrations modulate channel behaviour for the Cl⁻ (ANO1) and the K⁺ (BK) as well as the contractile force generating submodel illustrated in the centre. Mechanisms excluded from our final Reduced Tong Model (RTM) are noted with red crossouts.

https://doi.org/10.1371/journal.pcbi.1011359.g001

Step (i) hinges on the observation that conductance of some channels in the FTM are substantially larger than others; see Table 1. Experimentation with a select reduced suite of channel submodels and variations of conductances suggested elimination of several K⁺ channels in the FTM; we detail this below in the section on Elimination of Channels based on Conductance.

Download:

Table 1. Conductances for the FTM.

Items as taken from published values given in Tong, et al., 2011 including channels in the Tong, et al. 2014 version [16]. Ranges provided where suitable or used in original publication(s). Values in the CellML implementation of the model for reproduction of Tong, et al.’s Fig 12 are also shown, as well as equation numbers for the 2011 model (see the physiome repository). Final column indicates inclusion in the RTM; see text for further details.

https://doi.org/10.1371/journal.pcbi.1011359.t001

Step (ii) is inspired by the observed range of time scales. Ion channel time constants associated with the FTM span six orders of magnitude, from 10⁻¹ milliseconds (the Ca²⁺ -modulated BK-channel) to 100 seconds (the K₁ channel), illustrated in Fig 2A. This suggests, that inspired by the FN reduction of HH [5] there is scope to replace dynamical variables—Ordinary Differential Equations (ODEs)—in the FTM with their steady-state formulations. This aspect hinges on the form of channel models we inherit from the seminal work of HH, whereby the ion current is computed thus for a generic ion ‘x’: (1) for a conductance of g_x, where activation and inactivation variables are m_x and h_x, respectively. Driving force is given by the difference of membrane voltage (V_m) from the, typically, Nernst equilibrium for ion ‘x’ (E_x), and valence n along with the gas constant R, temperature T and Faradays constant, F. The dynamic ODE variables m_x and h_x determining the current are classically described as: (2) the activation variable, m_x, and similarly for h_x. These variables relax at some time scale, τ_m, to the steady-state value m_∞, for instance. Often the steady-state and relaxation times are dependent on membrane voltage, V_m, as with the Na⁺ channel and possibly ions or ligands such as Ca²⁺ in the Ca²⁺-modulated BK K⁺ channel.

Download:

Fig 2. Relaxation time scales ion channel submodels.

Panel (A). Suite of τ values as given for the RTM family of ion channels and the corresponding activation and inactivation variables. Note extraordinary range on logarithmic scale over six orders of magnitude from 10⁻¹ milliseconds for the BK-channel xb variable (teal trace) to as high as 100 seconds for the K⁺ channel variable r₁ (dashed black trace). Each variation over V_m are computed from the functional dependencies as given in the FTM, and selected channels in the RTM version appear here as described in text. The rough threshold for consideration of steady-state approximation partitioning the suite of these variables is at 10 ms (open circles), that is also sensitive to the range of V_m encountered during simulation. Panel (B). Corresponding suite of τ values for the FTM but excluded mechanisms from the final RTM. Note extraordinary range of time scales but shifted up around two orders of magnitude compared with the RTM suite.

https://doi.org/10.1371/journal.pcbi.1011359.g002

Typically, numerical simulations over time scales of seconds or longer set time step sizes to around dt = 0.01 seconds or 10 milliseconds—depending on the dynamics of the system of interest—and can rise given enough numerical stability. However, organ-wide uterine functions are dramatically longer [9]. Therefore we considered substituting the dynamic variables (i.e., m_x) relaxing at or below a threshold (i.e., τ_m ≤ dt = 0.01 s) with their corresponding steady-state values (i.e., m_∞(V_m)). We examined these substitutions in turn for multiple ion channel submodels as noted below in section Steady State Approximation.

Summary of the Tong uSMC model

The complete FTM encompasses a total of nineteen ODEs, describing activation and inactivation variables for nine ion channels. Of these channels, the FTM includes two Ca²⁺ (L-type CaL, T-type CaT), one voltage-gated Na⁺, one Cl⁻ (the ANO1) and five K⁺ channels (K₁, K₂, K_a, BK and b, or background current) with two additional channels permeable to either both Na⁺ and K⁺ (the ‘h’ hyperpolarisation) or Ca²⁺ as well (the NSCC with no ODEs). Four additional K⁺ channels appears in the 2014 Tong, et al. model [16]. Membrane voltage, V_m, is then computed by (3) where C_m is membrane capacitance, I_ion is the current through all ion channels, I_exchanger is the current through exchangers (the Na⁺-Ca²⁺ exchanger, NCX, and the Na⁺-K⁺ exchanger, N/K), and I_app is the applied stimulus. Current through each of the channels/exchangers in Table 1 is represented by I_x, where x is that channel’s identifier. Note the FTM does not dynamically track Na⁺ or K⁺ concentrations; the only influence of K⁺ and Na⁺channels is by way of currents in the expression for V_m. Alternatively, [Ca²⁺]_i is tracked dynamically through action of channel influx and pump or exchanger efflux thus: (4) where extruder action such as the plasma-membrane Ca²⁺ ATP-ase (PMCA) and the NCX remove cytosolic Ca²⁺. Free cytosolic [Ca²⁺]_i is a small fraction of actual Ca²⁺ influx due to wide arrays of intracellular Ca²⁺ buffering mechanisms such as the mobile protein calmodulin or uptake via the Sarco/Endoplasmic Reticulum ATP-ase (SERCA) into organelle Ca²⁺ reservoirs. β represents this effect with a simple constant fraction and is set to 0.015, which they reduced from the cited source of Standen & Stanfield, 1982 [17]. Presumably, Tong, et al. lowered the Standen & Stanfield neuronal β for improved fits to uSMC data; we also consider β a free parameter for Ca²⁺ fitting purposes.

Standard Simulation Protocol (SSP)

We compare our RTM to the FTM by using what we define as the standard simulation protocol (SSP). We target the stimulated behaviours shown in Fig 12A of Tong et al. [12], which aimed to replicate experimental observations reported in Fig 2A of Okabe, et al. [18]. This SSP is an applied current of I_app = −0.5 pA/pF over 2 seconds and simulations of model responses out to a t_max of 10 seconds. All tests shown are with this strength and duration I_app as well as t_max except where otherwise stated. The FTM equations and parameters were provided by a CellML implementation of the Tong, et al. 2011 published model in the Physiome Model Repository (PMR [19]). The PMR provides curated versions of published cell models that reproduce results presented in the original publication. Table 1 lists conductances presented in both the original Tong, et al. 2011 publication and the CellML values utilised for reproduction of Fig 12A in Tong, et al. We utilised two numerical platforms: OpenCOR for importing the CellML model descriptions and full model simulations and comparisons between the FTM and RTM, and MATLAB for generation of I-V curves of the individual ion channel submodels. Numerical solutions of the FTM and RTM system of ODEs performed in OpenCOR [20] utilised CVODE with the backwards differentiation formula (BDF) time integrator and tolerances (absolute and relative) set to 1e⁻⁷ and a maximum adaptive dt set to 0.1 ms. MATLAB assembly of I-V curves utilised ode23s and both tolerances set to 1e⁻¹⁰, with maximum T as noted. All single-cell simulations presented were performed on a macbook pro laptop (2.8 GHz Intel i7 with 4 cores and hyperthreading).

Pilot simulations for 2D & 3D tissue geometries were performed on a 144 CPU compute workstation (3.1 GHz Intel Xeon 6254) using a single-thread and the Chaste [21] implementation of the CVODE solver with default options and adaptive time stepping. We utilised a simple 2D ‘tissue’ (1 mm square; 121 nodes with 200 quadrilateral elements) and a hollow 3D cylindrical tube approximating a rat uterine horn (length 10 mm and 0.3 mm thick with 1910 nodes in a 4074 tetrahedral elements). The SSP stimulus of -0.5 pA/pF was applied in a sub-region: 0.2 x 0.01 mm (2D) and 1.0 mm along the tubular axis (3D). Electrical signal diffusion was set at a rate of 0.013 mS/cm (2D) and 250.0 mS/cm (3D).

Model reduction

Elimination of channels based on conductance

Table 1 shows the K⁺ channels in the FTM. One variant (termed K⁺ Channel ‘1’, or K₁) clearly dominates the overall current: with nS/pF. The K⁺ Channel ‘a’ (K_a) is a distant second with nS/pF. Tong noted the K₁ as the most influential of all the K⁺ channels and Knock, et al. observed the K₁ exhibits 67% of total K⁺ current in the myometrium [22], while the K₂ at around 23%, and the Ca²⁺-triggered BK rounding out the remaining 10%. Cells were observed to deploy either the K₁ or the K₂, but overall the K₁ emerged as the more typically expressed in uSMC [22]. Atia et al. also predicted a significant role in uSMC for a functionally similar channel to the K1, known as the Kv2.1 [23], reinforcing the report of Knock and the model of Tong. We thus explore the potential for the K₁ (or Kv2.1) as a representative K⁺ channel for the entire K⁺ current, or I_K of Eq (3), in the FTM for our reduction effort. Note, we do not alter the dynamics of the K₁ channel model given to us by Tong et al.—only the conductance values are varied.

Eliminating K⁺ channel ‘a’.

We initially test disabling of the K⁺ current, K_a due to prominence of over remaining K⁺ channels without K₁ (see Table 1) and compare the simulated result of the FTM for the SSP (See Fig 3A). Deactivation of K_a presents a similar result to the FTM but with noticeable deviations. This includes an accelerated frequency of action potentials (APs), diminished V_m amplitude and higher peak [Ca²⁺]_i. Modifying the conductances for I_K1 and I_Na we obtain an improved fit, with slightly overshot peaks and valleys of amplitudes but excellent temporal alignment of APs. For I_K1, we increased above the CellML value to g_K1 = 0.598 nS/pF as well as raising g_Na to 0.07 nS/pF—still well within the range given in Table 1. Additionally, raising the constant fraction value for free unbuffered calcium, β, by a mere 2.5% tightly aligns the [Ca²⁺]_i trace during the test stimulus peaks and out to steady-state as shown.

Download:

Fig 3. K⁺ current inactivation.

(Panel A) Disabling K_a current. FTM result with identical test protocol to the SSP reproduced by the CellML implementation (blue trace). Both [Ca²⁺]_i and V_m (upper and lower Panes, respectively; note varying time scales) display similarity to FTM with K_a disabled and no modifications otherwise (black trace). Altered conductances for K₁ and I_Na of 0.598 and 0.07 nS/pF, respectively, along with increased fraction of free (unbuffered) Ca²⁺ by 2.5% (β = 0.0154) presents improved qualitative fit to FTM (red trace, dotted). (Panel B) Disabling K₂, K_a and h utilising the SSP as in Panel A. Disabling only K₂ (black trace), exhibits significant rise in [Ca²⁺]_i (upper pane), and truncated APs settling into elevated plateau V_m (lower pane) compared to FTM result (blue trace). Disabling K₂, K_a and h shows similar behaviour (red trace). Alterations to K₁ and I_Na of 0.655 and 0.05125 nS/pF, respectively, with β = 0.0146 improve qualitative fit to FTM (pilot-RTM, orange dotted). Experimental trace adapted from Tong, 2011 Fig 12A originally via Okabe, et al. [18] (inset, Panel B) showing pilot-RTM reproduces observed V_m peaks above 0 mV over 1 second. ‘Pilot-RTM’ here refers to interim reduced model; see text for more.

https://doi.org/10.1371/journal.pcbi.1011359.g003

Eliminating K⁺ channel ‘2’.

Although the conductance of K_a is substantially higher than that for the K₂ channel in the FTM (0.2 versus 0.04 nS/pF), the impact of disabling the K₂ is noticeably more significant (Fig 3B). Without the I_K2, under the SSP, APs are damped into a higher steady-state V_m and [Ca²⁺]_i rises considerably higher to around 450 nM. Inactivating I_Ka and I_h as well enhances these effects without substantial overall change; the K₂ channel is thus more influential despite the order of magnitude less conductance than the K_a. This is explained by the time scales for current activity of the K₂ and the K_a channels as presented in the Tong, 2011 current data (see their Figures 6E and 7C [12]). The K₂ remains active over the course of seconds whereas the K_a shuts down within 50-100 ms; hence the K₂’s apparently greater impact we see here. Nevertheless, modest modifications to channel conductances and the β produce the qualitatively good alignments for this ‘pilot-RTM’ between AP peaks and the corresponding stair-step rise in [Ca²⁺]_i concentration. We only present the interim RTM here including only conductance variations without any steady-state modifications as done below; this we denote the ‘pilot-RTM’. Moreover, the pilot-RTM results shown exclude the N/K and NCX exchanger influence on V_m (see crossouts in Fig 1), although J_NCX still actively extrudes [Ca²⁺]_i. Modest adjustments to the I_K1 and I_Na thus compensate for and reproduce the overall result of the FTM without the full suite of K⁺ mechanisms.

The experimental trace of Okabe et al. is shown for comparison (Panel B, inset Fig 3), as utilised in Fig 12A of Tong, et al. Our pilot-RTM result qualitatively reproduces well the reported Okabe trace [18]. A similar spike train of APs and peak AP levels exceed 0 mV over the 1 second duration of their stimulus test with similar conditions, whereas the FTM V_m levels all fall under 0 mV after the initial AP. However, temporally the AP peaks are slightly delayed for the pilot-RTM from the Okabe trace; this alignment is improved below with steady-state approximations.

Other K⁺ mechanisms.

Remaining K⁺ currents in the 2011 FTM such as the background, I_b, and the Ca²⁺-modulated BK sit at opposite ends of the conductance range. I_b is smaller than all others by an order of magnitude and the BK is at 30% of total (compare Table 1). Despite the marginal magnitude of g_b, it maintains resting V_m at around -58 mV acting as a collective K⁺ channel representative [12]. Disabling I_b slowly depolarises V_m until the L-type channel is triggered and bursting occurs. Alternatively, the NSCC also maintains resting V_m but in the opposite direction, by preventing V_m from falling to the equilibrium potential for I_b at E_K = −84 mV. Both functional forms of I_b and NSCC are computationally straightforward, and we thus retain both mechanisms for resting V_m maintenance in the RTM.

The BK, however, is a rather different concern. Deactivation of I_BK produces similar results as inactivation of I_K2, with truncated APs leading to a depolarised steady-state V_m during stimulus, and a significant rise in free [Ca²⁺]_i over the FTM under the SSP, that can be mitigated with adjustments to g_K1. However, the BK channel dynamic variables as displayed in Fig 2A relax faster than any other mechanism in the FTM. Both the τ_xa and τ_xb scales are well under 10 ms over the entire range of -100 to 50 mV—beyond the physiological V_m for uSMC. This rapid relaxation appears ripe for exploitation and retention of the BK in our RTM. Particularly, of great interest to multi-scale modeling efforts are Ca²⁺ dynamics in the uSMC and Ca²⁺ modulation of the BK, the ANO1, contractile force generation, and myriad other elements subject to Ca²⁺ fluctuation. Assembly of an RTM that captures essential behaviour while reducing computational complexity by elimination of ODE variables in the model achieves this aim quite effectively, which we describe in more detail next with Steady-state Approximation.

Steady-state approximation

Uterine organ-wide functions in clinical settings occur over seconds to hours and months [24] that vastly eclipse behaviour of individual channel subunits. Organ scale models typically attempt to reflect clinical measurements, and so the timescales of importance in these models are seconds to hours and beyond.

The K⁺ BK channel.

The Ca²⁺-modulated K⁺ channel, or ‘big’ current BK, adapts rapidly to rising V_m with relaxation τ for the BK dynamic variables, x_α and x_β, reaching at most about 4 ms (V_m ≤ 50 mV, see Fig 2A). These τ are further independent of [Ca²⁺]_i and the steady-state (SS) depictions of the BK current are dependent only on V_m: (5) for fractional activation of the two variables x_α and x_β by the proportions p_α and p_β, respectively, and a Ca²⁺ dependency, z_i, not detailed here. Dynamic ODEs tracking these two activation variables, x_α and x_β, are given by the standard formulation (6) for some time span t ∈ [0, T]. We omit numerous details such as with z_i computed in ; see Tong, et al. for complete formulations. Clearly, the ratio of T and τ_i determine accuracy of approximating x_i with . Utilising instead of the full ODE when simulating beyond about T = 10 ms should be adequate for seconds-scaled simulations given dt at around 0.01 seconds, and the extraordinarily rapid relaxation times of τ_α and τ_β under 10 ms.

Numerical scans solving the ODE of Eq 6 over V_m ∈ [−100, 60] out to T = 10 ms display excellent agreement between full ODE solutions for both x_α and x_β and their corresponding SS approximations (see Fig 4A), computing current, I_BK, with a normalised conductance, i.e., as given in Eq 5. Deviation for I_BK (ΔI_BK) over the extended V_m range out to 200 mV (well above experimental observation) is at most about 0.2 (relative, normalised), and machine-ϵ errors for physiological V_m below about 20 mV. FTM results for the SSP utilising these reproducing the result of Tong, Fig 12 are virtually identical to the full BK ODE model, with relative ΔV_m under 1% throughout for all combinations of SS approximations. The τ for both BK variables that are comfortably under 10 ms within physiological ranges of V_m ∈ [−80, 20] mV for uSMC as noted in Table 2 enable this result.

Download:

Fig 4. Steady-state approximation comparisons with full ODE solutions for I − V curves and V_m results from FTM simulations; BK (Panel A) and T-type (Panel B) channels.

Full ODE solutions displayed with black traces (labeled ‘FTM’) and hybrid SS-ODE or full SS solutions with dotted and/or coloured traces as noted. For currents of channels, numerically-solved ODEs and SS approximation computations for each channel performed over a range of V_m for t ∈ [0, T] with T = 10 or 100 ms as noted, and corresponding currents, I_x, shown, with conductances, all set to 1 (Panels A, B, upper panes). Resulting I_x for each channel scaled by maximum current computed with ODE solutions: outward for BK and inward for T-type. SSP simulated V_m presented for both full ODE and SS variables of the FTM (Panels a, b, lower inset panes) with ΔV_m for each result as shown. (Panel a) BK channel variables x_α, x_β computed either with ODE solver (black trace) or SS approximation for current (scaled I_bk, upper pane) and relative difference in current (ΔI_bk, upper pane inset) for comparison. Note legend traces apply to both panes: blue trace for both x_α and x_β approximated with SS, and red and cyan dotted for each SS-approximated variable in solo. Lower pane presents ΔV_m for SSP results (V_m shown in inset) with combinations of full ODE (FTM, black trace) and SS-approximations as annotated. (Panel b) CaT channel variables b and g. I_cat (scaled) in upper pane over scanned V_m range with T = 10 ms and T = 100 ms (inset) for full ODE (black trace, FTM) and SS-approximation for both b and g (blue dotted) or each variable SS-approximated in solo as noted (red, cyan). Legend traces apply to both panes. Lower pane: ΔV_m for SSP simulation with CaT channel either full ODE (FTM, black trace) or combinations of SS-approximations as noted; corresponding V_m results for SSP of each ODE-SS configuration displayed in inset. See SSP description for numerical details.

https://doi.org/10.1371/journal.pcbi.1011359.g004

Download:

Table 2. RTM Channel subunit relaxation time scales per variable.

Two ranges of V_m presented: beyond the physiological range of uSMC ([-100,200]) and within ([-80, 20]). Maximum values given correspond to plots shown in Fig 2.

https://doi.org/10.1371/journal.pcbi.1011359.t002

The T-Type Ca²⁺ channel.

Relative to the L-type , T-type conductance is substantially lower (see Table 1) manifesting currents at around two orders of magnitude less for CaT than the L-type during FTM SSP simulations. Doubling still holds I_CaT at one-tenth of I_CaL. Interestingly, setting the already low to zero effectively eliminates two AP pulses from the SSP FTM simulation and reduces [Ca²⁺]_i considerably by ∼100 nM. The CaT maintains a trickling current reversing roles of significance at resting V_m, with I_CaT ∼ −0.05 pA/pF and I_CaL ∼ −0.007 pA/pF. Thus, the CaT shapes behaviour at resting and stimulated V_m, and we hence explore its steady-state simplification.

Analogously for the CaT channel with activation variable b and inactivation g, with their corresponding SS and τ dependencies on V_m: (7) the τ_b remains below 10 ms over a wide range of V_m (Fig 2A; Table 2). Hence, the SS-approximation for activation b computing I_CaT (with normalised conductance) and the full dynamic model is an excellent fit (Fig 4B). However, below resting V_m at around -55 mV, SS approximated currents deviate significantly for the inactivation variable g either in solo or combined with the activation variable b. Given τ_g is well over 100 ms and rises even higher at hyperpolarised V_m, this is unsurprising (Fig 2). Rather unexpected, however, is the agreement for full simulations of the FTM with all these combinations of SS and ODE variables (Fig 4B, lower pane) where ΔV_m is at most about 0.1. Extending solution scans over V_m of the inactivation variable out to T = 100 ms sharply improves the fit for I_CaT at or above resting V_m with all SS approximation combinations (upper pane inset, Fig 4B). Even the combined b and g SS approximation—with no ODEs dynamically simulating the T-type behaviour—we obtain excellent reproduction of I_CaT from resting V_m and up with ΔI_CaT falling dramatically from O(0.1) for T = 10 ms to O(0.01) with T = 100 ms. Permitting an adaptive time-stepper to increase dt to 0.1 seconds during full FTM solutions, we thus obtain the overall good representation of V_m for the SSP as shown.

The L-Type Ca²⁺ channel.

Essential to both uSMC function and hence any uSMC model, the L-type Ca²⁺ channel submodel has three variables, one activation (d) and two inactivation (f₁, f₂). Time constants for activation d and one inactivation f₁ are at or below 10 ms; f₂ meanwhile remains well above to nearly 100 ms over a wide range of V_m (Fig 2A). Currents, I_CaL, are produced with the expressions (8) with τ_f1 fixed at 12 ms, and f_Ca a Ca²⁺-modulating function desensitising the L-type by shifting V_m-dependent activation to higher V_m with increasing [Ca²⁺]_i (see Tong, et al. for full formulations). Both inactivation variables relax to the same steady-state function f_ss at their respective given proportions. f₂ with its slow τ_f2 around resting V_m should give a poor SS representation and we thus only tested f₂ approximated by f_ss in combination with the other two variables—that indeed performs inadequately at T = 10 ms (see Fig 5A). Only the activation variable, d, at the shorter time-span T = 10 ms displays reasonable SS approximation that further improves with T extended to 100 ms. However, the f₂ SS comparison retains unacceptable deviations from the ODE-solved I_CaL even at this longer timescale, and thus informs the FTM ODE-SS variants tested next.

Download:

Fig 5. Steady-state approximation comparisons with ODE solutions for I − V curves and V_m results from FTM simulations; L-type (Panel A) and Na⁺ (Panel B) channels.

Full ODE solutions displayed with black traces (labeled ‘FTM’) and hybrid SS-ODE or full SS solutions with dotted/coloured traces as noted. Presented plots similar as for BK and T-type, with upper panes displaying currents (scaled to maximum inbound) solved out to T = 10 or 100 ms (inset), and lower panes V_m or ΔV_m as labeled. (Panel A) CaL channel, d activation and f₁, f₂ inactivation variables. I_CaL shown (upper pane) for full ODE (black trace, FTM) and SS approximated variables as noted (blue-, red- and cyan-dotted). SS approximations for only f₂ not performed. V_m presented (lower pane) and not ΔV_m due to obvious deviations; note same trace legend as upper pane (black FTM and coloured SS combinations). SSP simulations with varied conductances (inset) for two SS approximations and improved V_m fit to FTM: d and f₁ both SS (cyan trace) with g_CaL reduced; only d approximated with SS (red trace) with g_CaL reduced as well as g_Na and g_K1 both increased (see text for details). (Panel B). Na⁺ channel, m and h activation/inactivation variables. Upper pane presents currents for full ODE, hybrid SS and full SS solutions out to T = 10 and 100 ms (inset), and lower pane result for SSP showing ΔV_m (relative) and V_m (inset) with same ODE and SS combinations as noted in upper legend: full ODE (black trace), both m and h activation/inactivation variables SS (blue dotted) and m or h SS (red and cyan dotted, respectively). See SSP description for numerical details.

https://doi.org/10.1371/journal.pcbi.1011359.g005

Modest variations in conductances for the better performing combinations (d alone or combined with f₁) improved the reproduction of the full FTM result (Fig 5A lower pane, inset), but required introduction of non-zero Na⁺ conductance. The CellML implementation reproducing Tong, et al.’s Fig 12 runs with an inactivated Na⁺ channel (see Table 1); substantial variations to g_CaL and g_K1 stubbornly fail at improving the SS approximations as shown. Alternatively, we do obtain considerable improvements in the qualitative fits with modest variations to conductances in two variant combinations of SS variables: (i) activation d alone with g_Na = 0.0703, g_K1 = 0.535 and g_CaL = 0.39 nS/pF; and (ii) d and inactivation f₁ with g_Na = 0.0703, g_K1 = 0.60 nS/pF and g_CaL at the baseline CellML value of 0.6 nS/pF. All these varied conductances still fall well within physiological ranges as discussed by Tong, et al.

The voltage-gated Na⁺ channel.

The conductance of the voltage-gated Na⁺g_Na is at most around 1/6 of g_CaL yet shapes the model response to stimulus current, similar to the T-type channel. Raising g_Na from zero (set in the CellML implementation) to the range observed by Tong, et al. increases the number of APs then damps V_m response to a depolarised level with successively larger increases of g_Na up to 0.1 nS/pF. We thus include the Na⁺ channel in our RTM, and given relaxation times of the activation variable m of at most 7.2 ms (during depolarisation), a SS approximation for its action looks promising. Current computed for Na⁺ channel includes m and the inactivation variable h: (9) where τ_h sharply peaks to nearly 1,000 ms around resting V_m, and diminishes prospects for a successful SS approximation of h. Current solutions comparing the full ODE and SS evaluations with either both variables or each in solo reflect these relaxation time scales, with m SS closely mimicking the full ODE I_Na, but not h for T out to 10 ms (Fig 5B, upper pane). Extending T to 100 ms improves all approximations, yet still leaves a deviation around resting V_m for either combination with h – as expected given the significant rise for τ_h at that V_m value.

Regardless, deploying these SS evaluations in substitute of the full ODE I_Na calculation for SSP comparisons with the FTM shows excellent qualitative agreement with all combinations of m and h (Fig 5B, lower pane inset). Relative deviations from V_m, or ΔV_m, is more revealing, however, of errors rising to order 0.1 for all variations of SS-ODE m and h. Closer examination of the SSP result shows peaks varying somewhat for each SS-ODE combination, with m-SS alone slightly temporally advanced and h-SS combinations slightly delayed from the full ODE FTM simulated APs. These slight time-shifts are enough to result in the rather high ΔV_m observed. Given the qualitative agreement using both SS m and h, we thus include h as a SS approximation in the RTM.

Thus far, of the four channel submodels considered here for SS approximation, the BK, T-type, L-type and V_m-gated Na⁺, all include nine activation and inactivation variables. Of these variables, only the L-type inactivation f₂ appears unsuitable for SS approximations with its relaxation time reaching up to about 100 ms and the inactivation h for the Na⁺ channel peaking to nearly 1 second (Fig 2). Preliminary results support exclusion of f₂ but inclusion of h for SS substitutions, however, likely due to the peak of τ_h confined to V_m around resting (Fig 5). Alternatively, notice τ_f2 remains well above 10 ms throughout the AP event and hence performs poorly with a SS representation. We thus apply a SS approximation to eight of the dynamic variables and eliminate in turn eight ODEs from the FTM in our RTM. This results in a reduction from 19 ODEs in the FTM into 5 for the RTM thanks to six variables removed using I_K1 as representative of three other K⁺ channels and SS approximations eliminating eight. Twelve less dynamic variables in the system should hold significant impact on the computational performance of the RTM which we consider next.

Results

RTM performance in reproducing experiments

Individually, each modification to the FTM presented so far is promising, requiring either no or modest alterations to the model parameters listed in Table 1 as provided with the CellML implementation. It is the collective behaviour that is of most import, however, and here we present each modification explored above into a final assembled RTM aimed at reproducing the FTM behaviour but simplified and with reduced computational expense.

The SSP.

Initially, we show effectiveness of the RTM in reproducing the SSP so far utilised as a test comparison with the full FTM, and present V_m, [Ca²⁺]_i and force generation for the RTM (Fig 6). The FTM calculates force as dynamically dependent on [Ca²⁺]_i with a familiar ODE tracking deviation from a steady-state: (10) with half-maximal activation at around 160 nM [Ca²⁺]_i, and Hill coefficient about 3.6 (via the CellML implementation). Unfortunately, the time scale, τ_ω, is quite long at 4 seconds and higher as [Ca²⁺]_i rises during stimulus, precluding application of a SS approximation for ω here.

Download:

Fig 6. SSP RTM reproduction and comparison with FTM.

(Panel A) Complete RTM (red dotted) with all modifications detailed in text presented in comparison with FTM (blue trace), showing [Ca²⁺]_i (upper pane), V_m (middle pane) and computed Force (lower pane) as dependent on [Ca²⁺]_i (see text). (Panel B) Relative difference for results presented in Panel A. Note semilog scales for all. Key RTM parameters include: g_CaL = 0.6, g_Na = 0.0895, g_K1 = 0.7 pA/pF and β = 0.0169.

https://doi.org/10.1371/journal.pcbi.1011359.g006

Nevertheless, we obtain qualitatively excellent reproduction of the FTM with our RTM for this SSP instance (Fig 6A), with reasonable relative deviations throughout (Panel B). The RTM again excludes three K⁺ channels (Sec. Elimination of Channels Based on Conductance) as well as the N/K exchanger, and eight ODEs using SS approximations (Sec. Steady-state Approximation) including the inactivation h for the Na⁺ channel. Modest alterations to the conductances listed in Table 1 as noted in Fig 6 produced this result. Although V_m peaks and valleys of the RTM moderately exceed or overshoot the FTM during stimulus and recovery, resulting impact on [Ca²⁺]_i where the RTM falls below the FTM is mitigated with a straightforward alteration to β. Raising the free-buffer fraction of [Ca²⁺]_i up from 0.015 to 0.0169 enabled the tight correspondence for the [Ca²⁺]_i presented, and provides the superb result for force generation as computed from [Ca²⁺]_i. This raised β moreover still remains lower than that used by the Tong, et al. cited model of Standen & Stanfield with their β set to 0.05 for neuronal axons [17].

Plateaus, bursting and sustained APs.

The promising fit to the SSP with the RTM may be due to the SSP utilised as a test case for each modification considered in turn above. We thus exposed the RTM to additional scenarios, specifically, Tong et al.’s Fig 11 of 2011 [12] and Fig 8 of Tong, et al.’s 2014 [16]. Varying I_app from a strong down to a mild current produces the array of AP shapes and responses shown in Fig 7. For high-enough I_app, a plateau-type AP forms with V_m settling into an elevated steady state comparable to observation of Wilde & Marshall [25] (Fig 7A). Lowering I_app then generates a burst of APs that plateau at depolarised V_m analogous to Meisheri, et al. [26], that also emerges with a depolarising ‘pulse’ of initially elevated extracellular K⁺ and no applied stimulus (Fig 7B). Further drops of stimulus strength demonstrates sustained APs for the duration of I_app either over 10 or 50 secs (Fig 7C and 7D), analogous to Wilde & Marshalls report [25]. Each conductance and parameter of the RTM is moreover identical to those utilised for the SSP comparison result except for the depolarising K⁺ pulse. Only one other exception occurs, however, with the longer-duration I_app over 50 s requiring a higher g_K1 of 0.75 pA/pF.

Download:

Fig 7. RTM: Variation of I_app.

Comparison results from Tong, 2011 Fig 11A, Fig 11B and Fig 11C and Tong, 2014 Fig 8D. Stimulus strength varied as noted aimed at reproducing experimental results (insets, Panels A-C) as given in Tong, 2011 Fig 11, or extended bursting (Panel D). Each I_app applied up to t = 10 s starting at t = 1 s, except Panel D out to t = 50 s. Additional simulation with elevated extracellular initial condition K⁺ at 10 mM for comparison and no applied stimulus (blue dotted, Panel B). Conductances for key currents: g_CaL = 0.6, g_Na = 0.0895, g_K1 = 0.7 pA/pF (identical to SSP values, Panels A-C), with g_K1 = 0.75 for Panel D. Experimental traces (insets) adapted from Tong, et al. 2011 Fig 11.

https://doi.org/10.1371/journal.pcbi.1011359.g007

Over a longer duration I_app, the 2011 FTM failed to generate the AP bursting profile shown here produced by the RTM in Fig 7D. This complication with the 2011 model version inspired later additions of the four K⁺ channels noted as such in Fig 1 by Tong, et al. with their 2014 model [16]. Increasing the g_K1 as described by Tong, et al. from the CellML value of 0.52 to 0.64 pA/pF improved the FTM response, elongating an otherwise short burst of APs before damping into the plateau phase (see Fig 2 of [16]). However, the FTM required rescaling the already quite slow K1 time constants by substantial factors, i.e., 2 × τ_q, 20 × τ_r1 and 20 × τ_r2—pushing them as high as 5e6 ms (see Table 2)—in order to generate the longer AP burst we obtain here by simply increasing g_K1 to 0.75 pA/pF in our RTM.

Computational performance

For solutions of the FTM in SSP conditions and each variant of reduction detailed in methods—substituting one K⁺ channel for three others and SS approximations for ion channel submodels—and the full RTM encompassing all these reductions, we present performance details in Table 3. Over a sample of n = 10 runs, the average time-to-solve per variant improves modestly for the representative K₁ variant, and increases for each channel submodel with SS approximations all to around the same improvement factor of 1.7. Combined together into the final assembled RTM, we observe that the RTM is, on average, 2.1 times faster to solve for a single cell than the FTM. The individual factors of improvement are not cumulative. This is likely due to internal mechanisms beyond user control with the CVODE solver.

Download:

Table 3. Comparison time-to-solve the SSP.

FTM performance and each variant of reduction including conductance replacement (K₁ representative of K2, Ka and h), individual ion channel submodel SS approximations as noted, and the full RTM with all reduced items presented. Solutions of the single cell SSP performed in OpenCOR using CVODE with maximum dt = 0.1, tolerances (relative/absolute) at 1e⁻⁷, with n = 10 trials each on a macbook pro laptop. Solutions of 2D and 3D tissue SSP trials performed using CHASTE with CVODE (dt_max = 0.1 with default settings) on a single-thread of Intel Xeon 6254 CPU with n = 10 trials. Average time to solve (μ) and standard deviation for the set of trials (σ) provided as noted. Percent improvement calculated as change relative to results with FTM, i.e., abs(t_FTM − t_RTM)/t_FTM.

https://doi.org/10.1371/journal.pcbi.1011359.t003

Pilot simulations on both 2D and 3D geometries comparing the FTM and RTM in a monodomain tissue model were conducted to assess scalability—as described in Methods. The RTM again proves faster with factor improvements of 2.6 and 3.7 speedups in the 2D and 3D geometries, respectively (see Table 3), using the default CVODE options in CHASTE [21]. This follows from the reduced suite of ODEs in the RTM. For the 2D simulation, a total of 121 nodes × 5 ODEs gives 605 unknowns per time step versus 121 × 19 or 2299 per time step in the FTM. Similarly for the 3D simulation, the RTM solved 9550 unknowns versus the FTM’s 36,290 per time step.

Discussion

We presented our steady-state approximation method for reducing the order of a uSMC model by considering the multiple time scales inherent in a published model by Tong, et al., 2011 [12]. Our final assembled RTM demonstrates both accurate reproduction of experimental results while substantially improving computational performance. Arriving at this RTM does require careful adjustments of parameters such as conductance strengths or buffering fractions and numerical experimentation. Nevertheless, it is a straightforward process, particularly compared with formal sensitivity and bifurcation analysis [8]. Evaluating overall time-scales for each submodel within the context of the system of interest is essential. Importantly, in the uSMC model considered here, the activity of the physiological system is over multiple seconds and longer while individual channel dynamics may settle to steady-state at orders of magnitude faster [9].

Observing multiple time scales in a model system and exploiting steady-states for relatively fast sub-mechanisms is not unknown and indeed formed the basis for many systems including the widely-used FN neuronal AP model [5]. The seminal ventricular cardiac model of Luo and Rudy, for instance, incorporated a steady-state activation function for their time-independent K⁺ channel [27, 28]. With a τ of 0.7 ms, it operates considerably faster than other channels in their model with inactivations relaxing up to 100 ms. Moreover, the Luo & Rudy sequence of papers simulate cardiac behaviour over hundreds of milliseconds which is quite suitable for a cell model expressing ion channel time scales of 10 to 100 ms. Subsequent multi-scale models of the heart built on the Luo & Rudy cell model investigate similar temporal ranges [29], as do other multi-scale efforts founded on different cell models. Recent organ models such as from Margara, et al. [1] deploy the Tomek, et al. ventricular myocyte model [15] itself based on the Grandi, et al. model [30] that includes an ODE with a τ for a Na⁺ channel relaxing within 0.17 ms. Resolving fully three orders of magnitude faster than the simulated scales of the organ model, such an ion channel mechanism is a ripe candidate for steady-state approximation.

Cardiac electrophysiology models using bidoman/monodomain methods [31] laid the groundwork for their application to other organs, such as the gastrointestinal system [2]. Coupling partial differential equations for electrical diffusion with systems of intracellular ion dynamics, these models are computationally demanding; any simplifications are eagerly sought and cell-level models are a natural target. For instance, the Corrias & Buist model of the gastric pacemaker, the interstitial cells of Cajal (ICC), incorporates relaxation times spanning from just under 10 ms (for Na⁺and K⁺ channels) up to nearly 100 ms (L-type Ca²⁺) [14]. Du. et al. simplified the Corrias & Buist model in their later organ-scale model of gastric electrophysiology, investigating slow-waves of activity over scales of seconds up to minutes [32]. Yet, the ICC Na⁺ channel activation variable with a τ of 10 ms remained in the full multi-scale model [2]. Approximating such rapid millisecond-scaled dynamic mechanisms with their steady-state evaluations in a seconds-to-minutes simulation is easily overlooked. As demonstrated here, such approximations—for suitable time scales of interest—hold considerable potential for improving computational performance.

Sensitivity & bifurcation analysis sifting out model components while maintaining essential dynamics as applied to the uSMC model of Rihana, et al. [11] are certainly important and valid approaches. Reduced from the original ten dynamical variables, Laforet, et al. utilised the reduction of Parthimos et al.’s arterial model [8] and the Rihana, et al. reduction [33] to produce a three-variable uSMC cell model [34] slated for later use in that group’s full 3D multi-scale model of the uterus [3]. The relative simplicity of their minimal uSMC cell model enabled greater flexibility to explore otherwise computationally infeasible simulations—but only after substantial exertion for the cell model analysis. Of note, we considered performing a steady-state reduction of order with the Rihana, et al. model, but neither that nor subsequent versions published parameter values for the relaxation τ’s [11].

With the Tong, et al. 2011 model, after removing several K⁺ channels, using the K₁ as a representative K⁺ channel, and the steady-state approximations shown, we eliminated eight dynamical variables out of the total thirteen in our RTM. The relative conductance strengths listed in Table 1 informed selection of candidates for elimination starting with the second-strongest conductance channel, K_a, and progressing through the rest. Interestingly, despite nearly an order of magnitude stronger conductance than , K₂ emerged as the more relevant channel (compare Fig 3B with Fig 3A). K₂ relaxes over the course of seconds while the K_a does within about 50 ms (compare Fig 6E and 7C, [12]). Hence, both conductances and time scales dictate suitability for this representative channel approach.

Observing the relative conductance strengths in any suitable system—with an eye to relative dynamic behaviour—permits reproduction of this representative channel approach in the following manner: (i) Enumerate conductances according to strength; (ii) Evaluate relative contribution to dynamics; and (iii) Adjust dominant or representative conductance. Certainly, this method may be combined with a sensitivity and bifurcation analysis to reduce our RTM even further, with the corresponding effort involved.

Meanwhile, the relative simplicity of comparing time scales for individual sub-model variables with the overall time scale of interest also suggests the possibility of automation. Algorithmically, determining suitability of steady-state approximations within a modelling framework appears a pragmatic possibility, perhaps as a precursor to a full reduction analysis, and may be similar to the following procedure:

Set time scale of interest, T_i, (here in seconds);
Set membrane voltage range of interest, V_span, (here ∼ [−50, 50] mV);
Identify potential variables for SS approximation:
- If τ_x ≤ 0.01 × T_i for V_m ∈ V_span
- Test corresponding variable, X, for SS performance in channel model:
  - Select time step size as appropriate for comparison; i.e., dt = 0.01 ms;
  - Perform full ODE simulation channel model;
  - Compare ODE result with SS variable approximation.

The above procedure simply replicates what we conducted here in our SS explorations, including dynamic variables that performed well with SS approximations despite relaxation rates well exceeding our ‘rule of thumb’ guide for 10 ms. As seen in the plots for the Na⁺ variable h and the T-type variable g in Fig 2, this is due to the peak τ levels sitting at the lower end of V_span. Hence, by specifying these temporal and voltage spans of interest, an automated system may identify potential dynamic variables for SS approximations and evaluate their suitability with full ODE and SS comparisons comparable to those conducted here.

Limitations

Representative channel usage such as the K₁ in our RTM here, does entail some disadvantages, however. Although for our purposes we reproduced the test cases shown for the SSP and variations of I_app, there may be scenarios where the RTM will fail without the full bank of K⁺ channels in the FTM. By the same token, there may also be scenarios where the FTM would fail given the extraordinary array of K⁺ channels discovered in the exhaustive transcription maps of Atia, et al., for instance [23]. These possibilities are simply unknown, particularly given the uncertainties and rather interesting questions as to why so many ion channels with their subtle or dramatic differences are expressed in uSMC, or any cell type for that matter.

Modeling implications: Realism and relevance

Our reduction reinforces how models replicating measured data do not require every mechanism known to reside in a cellular system. Note the expanded 2014 FTM utilising nine K⁺ channels still omitted numerous examples of the 22 channels reported by Atia et al. [23]. Regardless, the FTM—without this full suite of channels—reproduced the measured data that we also accomplished here with even fewer mechanisms in the RTM. Certainly, excluding channels observed in a cell reduces the molecular fidelity of any mathematical model regardless of the cellular system under study. However, fidelity must necessarily be balanced against tractibility and relevance.

Mathematical models of the biological system are constrained by the available measured data; our ‘minimal model’ presented here reproduces the available data. Of course, as with any model, if in the future the model failed to reproduce new data, molecular mechanisms could be reintroduced systematically to determine which drive the particular novel behaviour. This is a particular benefit of a reduced modeling strategy, in that it provides important insights into the required mechanisms for different cell behaviours, teasing out their relative influences.

An example of a mechanism that may be important for some behaviours and not others is the ryanodine receptor Ca²⁺ channel (RyR). It is clearly expressed in the uSMC, but also appears to be functionally inert for the behaviours observed [35, 36]. Reports of the RyR’s appearance in the uSMC likely inspired its incorporation in a Ca²⁺ dynamical uSMC model that also excluded the other intracellular Ca²⁺ channel, the inositol trisphosphate receptor (IP3R) [37]. Unlike the RyR, the IP3R not only appears in the uSMC, it does appear influential on uterine contractile activity [9]. Striking a careful balance between accuracy and tractability, ‘minimal models’ aid illuminating the relevance of resident cellular mechanisms but of course are always constrained by the experimental scenarios providing the data.

Time scales and performance

Suitability of the steady-state approximation approach is of course dependent on the time scales of interest, or T_i. For overall behaviour of an organ-scale model with relatively similar time scales on the ion channel level as, say, with the Luo & Rudy cardiac model, a steady-state substitution may fail to perform. Alternatively, with longer organ-scale phenomena as with the uterus over multiple seconds and more, the extraordinarily rapid ion channel relaxation times as in the FTM prove excellent candidates. Careful consideration of the relative time scales and particularly the context for the overall system of interest naturally dictates whether elimination of dynamic variables will succeed—as described above.

If indeed the relative time scales prove suitable as demonstrated here with the RTM, steady-state replacements of full dynamical variables combined with modest calibrations to parameters in the model readily provide substantial improvements in computational speed while maintaining accuracy. For a single cell model, a factor of two reduction in time to solve may be negligible—150 ms versus 300 (Table 3). However, distributed over significant numbers of grid points in a full organ-scale model this difference is amplified considerably, potentially producing significant reductions in absolute computation times. Hence, we trialed pilot simulations of a monodomain tissue model on simple 2D and 3D geometries (square and hollow cylindrical tube, respectively) with the improved times-to-solve as noted exceeding the single-cell simulation improvements (Table 3). Of course, different implementations, solvers and settings will produce different results. Regardless, reducing the number of unknowns per time step by roughly a factor of 4, as here with the RTM, will reduce computational load and enable greater flexibility in parameter sensitivity analyses and in silico investigations.

Conclusion

A reduction of order was performed for the Tong, et al. uSMC model utilising steady-state approximations for rapidly resolving dynamical variables. The relatively slow functions of uSMCs over the scale of seconds and longer provide the key context for eliminating whole ODEs with millisecond-level time scales. Examining these relative time scales and testing suitability for reproducing observed uSMC yields the RTM presented here. Demonstrating the effectiveness and relative ease of a steady-state reduction of order approach with this uSMC model test case underlines its utility for the numerous cell models deployed in myriad multi-scale representations, and holds great potential for expanding the computational possibilities for simulating the many intrinsically multi-scaled organs and systems.

Acknowledgments

The authors warmly thank Prof. James Sneyd, Mathematics, University of Auckland for informative and helpful discussions.

References

1. Margara F, Wang ZJ, Levrero-Florencio F, Santiago A, Vázquez M, Bueno-Orovio A, et al. In-silico human electro-mechanical ventricular modelling and simulation for drug-induced pro-arrhythmia and inotropic risk assessment. Progress in Biophysics and Molecular Biology. 2021;159:58–74. pmid:32710902
- View Article
- PubMed/NCBI
- Google Scholar
2. Du P, O’Grady G, Cheng LK. A theoretical analysis of anatomical and functional intestinal slow wave re-entry. Journal of Theoretical Biology. 2017;425:72–79. pmid:28450068
- View Article
- PubMed/NCBI
- Google Scholar
3. Yochum M, Laforêt J, Marque C. Multi-scale and multi-physics model of the uterine smooth muscle with mechanotransduction. Computers in Biology and Medicine. 2018;93:17–30. pmid:29253628
- View Article
- PubMed/NCBI
- Google Scholar
4. Fall CP, Marland ES, Wanger JM, Tyson JJ, editors. Computational Cell Biology. 1st ed. Interdisciplinary Applied Mathematics. Springer New York; 2002.
5. Fitzhugh R. Thresholds and plateaus in the Hodgkin-Huxley nerve equations. J Gen Physiol. 1960;43(5):867–896. pmid:13823315
- View Article
- PubMed/NCBI
- Google Scholar
6. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117(4):500–544. pmid:12991237
- View Article
- PubMed/NCBI
- Google Scholar
7. Ali MA, Means SA, Ho H, Heffernan J. Global sensitivity analysis of a single-cell HBV model for viral dynamics in the liver. Infect Dis Model. 2021;6:1220–1235. pmid:34786526
- View Article
- PubMed/NCBI
- Google Scholar
8. Parthimos D, Edwards DH, Griffith TM. Minimal model of arterial chaos generated by coupled intracellular and membrane Ca2+ oscillators. Am J Physiol. 1999;277(3):H1119–44. pmid:10484436
- View Article
- PubMed/NCBI
- Google Scholar
9. Garrett AS, Means SA, Roesler MW, Miller KJW, Cheng LK, Clark AR. Modeling and experimental approaches for elucidating multi-scale uterine smooth muscle electro- and mechano-physiology: A review. Front Physiol. 2022;13:1017649. pmid:36277190
- View Article
- PubMed/NCBI
- Google Scholar
10. Bursztyn L, Eytan O, Jaffa AJ, Elad D. Mathematical model of excitation-contraction in a uterine smooth muscle cell. Am J Physiol Cell Physiol. 2007;292(5):C1816–29. pmid:17267547
- View Article
- PubMed/NCBI
- Google Scholar
11. Rihana S, Terrien J, Germain G, Marque C. Mathematical modeling of electrical activity of uterine muscle cells. Med Biol Eng Comput. 2009;47(6):665–675. pmid:19301052
- View Article
- PubMed/NCBI
- Google Scholar
12. Tong WC, Choi CY, Kharche S, Holden AV, Zhang H, Taggart MJ. A computational model of the ionic currents, Ca2+ dynamics and action potentials underlying contraction of isolated uterine smooth muscle. PLoS One. 2011;6(4):e18685. pmid:21559514
- View Article
- PubMed/NCBI
- Google Scholar
13. Sheldon RE, Baghdadi M, McCloskey C, Blanks AM, Shmygol A, van den Berg HA. Spatial heterogeneity enhances and modulates excitability in a mathematical model of the myometrium. J R Soc Interface. 2013;10(86):20130458. pmid:23843249
- View Article
- PubMed/NCBI
- Google Scholar
14. Corrias A, Buist ML. Quantitative cellular description of gastric slow wave activity. Am J Physiol Gastrointest Liver Physiol. 2008;294(4):G989–95. pmid:18276830
- View Article
- PubMed/NCBI
- Google Scholar
15. Tomek J, Bueno-Orovio A, Passini E, Zhou X, Minchole A, Britton O, et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. Elife. 2019;8. pmid:31868580
- View Article
- PubMed/NCBI
- Google Scholar
16. Tong WC, Ghouri I, Taggart MJ. Computational modeling of inhibition of voltage-gated Ca channels: identification of different effects on uterine and cardiac action potentials. Frontiers in Physiology. 2014;5. pmid:25360118
- View Article
- PubMed/NCBI
- Google Scholar
17. Standen NB, Stanfield PR. A binding-site model for calcium channel inactivation that depends on calcium entry. Proc R Soc Lond B Biol Sci. 1982;217(1206):101–110. pmid:6131420
- View Article
- PubMed/NCBI
- Google Scholar
18. Okabe K, Inoue Y, Kawarabayashi T, Kajiya H, Okamoto F, Soeda H. Physiological significance of hyperpolarization-activated inward currents (Ih) in smooth muscle cells from the circular layers of pregnant rat myometrium. Pflugers Arch. 1999;439(1-2):76–85. pmid:10651003
- View Article
- PubMed/NCBI
- Google Scholar
19. Yu T, Lloyd CM, Nickerson DP, Cooling MT, Miller AK, Garny A, et al. The Physiome Model Repository 2. Bioinformatics. 2011;27(5):743–744. pmid:21216774
- View Article
- PubMed/NCBI
- Google Scholar
20. Garny A, Hunter PJ. OpenCOR: a modular and interoperable approach to computational biology. Frontiers in Physiology. 2015;6. pmid:25705192
- View Article
- PubMed/NCBI
- Google Scholar
21. Cooper FR, Baker RE, Bernabeu MO, Bordas R, Bowler L, Bueno-Orovio A, et al. Chaste: Cancer, Heart and Soft Tissue Environment. Journal of Open Source Software. 2020;3(47):1848. pmid:37192932
- View Article
- PubMed/NCBI
- Google Scholar
22. Knock GA, Smirnov SV, Aaronson PI. Voltage-gated K+ currents in freshly isolated myocytes of the pregnant human myometrium. J Physiol. 1999;518 (Pt 3)(Pt 3):769–781. pmid:10420013
- View Article
- PubMed/NCBI
- Google Scholar
23. Atia J, McCloskey C, Shmygol AS, Rand DA, van den Berg HA, Blanks AM. Reconstruction of Cell Surface Densities of Ion Pumps, Exchangers, and Channels from mRNA Expression, Conductance Kinetics, Whole-Cell Calcium, and Current-Clamp Voltage Recordings, with an Application to Human Uterine Smooth Muscle Cells. PLoS Comput Biol. 2016;12(4):e1004828. pmid:27105427
- View Article
- PubMed/NCBI
- Google Scholar
24. Bakker PCAM, Van Rijswijk S, van Geijn HP. Uterine activity monitoring during labor. J Perinat Med. 2007;35(6):468–477. pmid:18052832
- View Article
- PubMed/NCBI
- Google Scholar
25. Wilde DW, Marshall JM. Effects of tetraethylammonium and 4-aminopyridine on the plateau potential of circular myometrium from the pregnant rat. Biol Reprod. 1988;38(4):836–845.
- View Article
- Google Scholar
26. Meisheri KD, McNeill JH, Marshall JM. Effect of isoproterenol on the isolated pregnant rat myometrium. Eur J Pharmacol. 1979;60(1):1–6. pmid:520410
- View Article
- PubMed/NCBI
- Google Scholar
27. Luo CH, Rudy Y. A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circ Res. 1991;68(6):1501–1526. pmid:1709839
- View Article
- PubMed/NCBI
- Google Scholar
28. Luo CH, Rudy Y. A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes. Circ Res. 1994;74(6):1071–1096. pmid:7514509
- View Article
- PubMed/NCBI
- Google Scholar
29. Colli Franzone P, Pavarino LF, Taccardi B. Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models. Math Biosci. 2005;197(1):35–66. pmid:16009380
- View Article
- PubMed/NCBI
- Google Scholar
30. Grandi E, Pasqualini FS, Bers DM. A novel computational model of the human ventricular action potential and Ca transient. J Mol Cell Cardiol. 2010;48(1):112–121. pmid:19835882
- View Article
- PubMed/NCBI
- Google Scholar
31. Clayton RH, Bernus O, Cherry EM, Dierckx H, Fenton FH, Mirabella L, et al. Models of cardiac tissue electrophysiology: progress, challenges and open questions. Prog Biophys Mol Biol. 2011;104(1-3):22–48. pmid:20553746
- View Article
- PubMed/NCBI
- Google Scholar
32. Du P, Gao J, O’Grady G, Cheng LK. A simplified biophysical cell model for gastric slow wave entrainment simulation. Annu Int Conf IEEE Eng Med Biol Soc. 2013;2013:6547–6550. pmid:24111242
- View Article
- PubMed/NCBI
- Google Scholar
33. Rihana S, Santos J, Mondie S, Marque C. Dynamical analysis of uterine cell electrical activity model. Conf Proc IEEE Eng Med Biol Soc. 2006;2006:4179–4182. pmid:17946608
- View Article
- PubMed/NCBI
- Google Scholar
34. Laforet J, Rabotti C, Terrien J, Mischi M, Marque C. Toward a multiscale model of the uterine electrical activity. IEEE Trans Biomed Eng. 2011;58(12):3487–3490. pmid:21968708
- View Article
- PubMed/NCBI
- Google Scholar
35. Malik M, Roh M, England SK. Uterine contractions in rodent models and humans. Acta Physiol (Oxf). 2021;231(4):e13607. pmid:33337577
- View Article
- PubMed/NCBI
- Google Scholar
36. Matsuki K, Takemoto M, Suzuki Y, Yamamura H, Ohya S, Takeshima H, et al. Ryanodine receptor type 3 does not contribute to contractions in the mouse myometrium regardless of pregnancy. Pflugers Arch. 2017;469(2):313–326. pmid:27866274
- View Article
- PubMed/NCBI
- Google Scholar
37. Testrow CP, Holden AV, Shmygol A, Zhang H. A computational model of excitation and contraction in uterine myocytes from the pregnant rat. Sci Rep. 2018;8(1):9159. pmid:29904075
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Margara F, Wang ZJ, Levrero-Florencio F, Santiago A, Vázquez M, Bueno-Orovio A, et al. In-silico human electro-mechanical ventricular modelling and simulation for drug-induced pro-arrhythmia and inotropic risk assessment. Progress in Biophysics and Molecular Biology. 2021;159:58–74. pmid:32710902
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Du P, O’Grady G, Cheng LK. A theoretical analysis of anatomical and functional intestinal slow wave re-entry. Journal of Theoretical Biology. 2017;425:72–79. pmid:28450068
View Article
PubMed/NCBI
Google Scholar

[6] View Article

[7] PubMed/NCBI

[8] Google Scholar

[ref3] 3. Yochum M, Laforêt J, Marque C. Multi-scale and multi-physics model of the uterine smooth muscle with mechanotransduction. Computers in Biology and Medicine. 2018;93:17–30. pmid:29253628
View Article
PubMed/NCBI
Google Scholar

[10] View Article

[11] PubMed/NCBI

[12] Google Scholar

[ref4] 4. Fall CP, Marland ES, Wanger JM, Tyson JJ, editors. Computational Cell Biology. 1st ed. Interdisciplinary Applied Mathematics. Springer New York; 2002.

[ref5] 5. Fitzhugh R. Thresholds and plateaus in the Hodgkin-Huxley nerve equations. J Gen Physiol. 1960;43(5):867–896. pmid:13823315
View Article
PubMed/NCBI
Google Scholar

[15] View Article

[16] PubMed/NCBI

[17] Google Scholar

[ref6] 6. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117(4):500–544. pmid:12991237
View Article
PubMed/NCBI
Google Scholar

[19] View Article

[20] PubMed/NCBI

[21] Google Scholar

[ref7] 7. Ali MA, Means SA, Ho H, Heffernan J. Global sensitivity analysis of a single-cell HBV model for viral dynamics in the liver. Infect Dis Model. 2021;6:1220–1235. pmid:34786526
View Article
PubMed/NCBI
Google Scholar

[23] View Article

[24] PubMed/NCBI

[25] Google Scholar

[ref8] 8. Parthimos D, Edwards DH, Griffith TM. Minimal model of arterial chaos generated by coupled intracellular and membrane Ca2+ oscillators. Am J Physiol. 1999;277(3):H1119–44. pmid:10484436
View Article
PubMed/NCBI
Google Scholar

[27] View Article

[28] PubMed/NCBI

[29] Google Scholar

[ref9] 9. Garrett AS, Means SA, Roesler MW, Miller KJW, Cheng LK, Clark AR. Modeling and experimental approaches for elucidating multi-scale uterine smooth muscle electro- and mechano-physiology: A review. Front Physiol. 2022;13:1017649. pmid:36277190
View Article
PubMed/NCBI
Google Scholar

[31] View Article

[32] PubMed/NCBI

[33] Google Scholar

[ref10] 10. Bursztyn L, Eytan O, Jaffa AJ, Elad D. Mathematical model of excitation-contraction in a uterine smooth muscle cell. Am J Physiol Cell Physiol. 2007;292(5):C1816–29. pmid:17267547
View Article
PubMed/NCBI
Google Scholar

[35] View Article

[36] PubMed/NCBI

[37] Google Scholar

[ref11] 11. Rihana S, Terrien J, Germain G, Marque C. Mathematical modeling of electrical activity of uterine muscle cells. Med Biol Eng Comput. 2009;47(6):665–675. pmid:19301052
View Article
PubMed/NCBI
Google Scholar

[39] View Article

[40] PubMed/NCBI

[41] Google Scholar

[ref12] 12. Tong WC, Choi CY, Kharche S, Holden AV, Zhang H, Taggart MJ. A computational model of the ionic currents, Ca2+ dynamics and action potentials underlying contraction of isolated uterine smooth muscle. PLoS One. 2011;6(4):e18685. pmid:21559514
View Article
PubMed/NCBI
Google Scholar

[43] View Article

[44] PubMed/NCBI

[45] Google Scholar

[ref13] 13. Sheldon RE, Baghdadi M, McCloskey C, Blanks AM, Shmygol A, van den Berg HA. Spatial heterogeneity enhances and modulates excitability in a mathematical model of the myometrium. J R Soc Interface. 2013;10(86):20130458. pmid:23843249
View Article
PubMed/NCBI
Google Scholar

[47] View Article

[48] PubMed/NCBI

[49] Google Scholar

[ref14] 14. Corrias A, Buist ML. Quantitative cellular description of gastric slow wave activity. Am J Physiol Gastrointest Liver Physiol. 2008;294(4):G989–95. pmid:18276830
View Article
PubMed/NCBI
Google Scholar

[51] View Article

[52] PubMed/NCBI

[53] Google Scholar

[ref15] 15. Tomek J, Bueno-Orovio A, Passini E, Zhou X, Minchole A, Britton O, et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. Elife. 2019;8. pmid:31868580
View Article
PubMed/NCBI
Google Scholar

[55] View Article

[56] PubMed/NCBI

[57] Google Scholar

[ref16] 16. Tong WC, Ghouri I, Taggart MJ. Computational modeling of inhibition of voltage-gated Ca channels: identification of different effects on uterine and cardiac action potentials. Frontiers in Physiology. 2014;5. pmid:25360118
View Article
PubMed/NCBI
Google Scholar

[59] View Article

[60] PubMed/NCBI

[61] Google Scholar

[ref17] 17. Standen NB, Stanfield PR. A binding-site model for calcium channel inactivation that depends on calcium entry. Proc R Soc Lond B Biol Sci. 1982;217(1206):101–110. pmid:6131420
View Article
PubMed/NCBI
Google Scholar

[63] View Article

[64] PubMed/NCBI

[65] Google Scholar

[ref18] 18. Okabe K, Inoue Y, Kawarabayashi T, Kajiya H, Okamoto F, Soeda H. Physiological significance of hyperpolarization-activated inward currents (Ih) in smooth muscle cells from the circular layers of pregnant rat myometrium. Pflugers Arch. 1999;439(1-2):76–85. pmid:10651003
View Article
PubMed/NCBI
Google Scholar

[67] View Article

[68] PubMed/NCBI

[69] Google Scholar

[ref19] 19. Yu T, Lloyd CM, Nickerson DP, Cooling MT, Miller AK, Garny A, et al. The Physiome Model Repository 2. Bioinformatics. 2011;27(5):743–744. pmid:21216774
View Article
PubMed/NCBI
Google Scholar

[71] View Article

[72] PubMed/NCBI

[73] Google Scholar

[ref20] 20. Garny A, Hunter PJ. OpenCOR: a modular and interoperable approach to computational biology. Frontiers in Physiology. 2015;6. pmid:25705192
View Article
PubMed/NCBI
Google Scholar

[75] View Article

[76] PubMed/NCBI

[77] Google Scholar

[ref21] 21. Cooper FR, Baker RE, Bernabeu MO, Bordas R, Bowler L, Bueno-Orovio A, et al. Chaste: Cancer, Heart and Soft Tissue Environment. Journal of Open Source Software. 2020;3(47):1848. pmid:37192932
View Article
PubMed/NCBI
Google Scholar

[79] View Article

[80] PubMed/NCBI

[81] Google Scholar

[ref22] 22. Knock GA, Smirnov SV, Aaronson PI. Voltage-gated K+ currents in freshly isolated myocytes of the pregnant human myometrium. J Physiol. 1999;518 (Pt 3)(Pt 3):769–781. pmid:10420013
View Article
PubMed/NCBI
Google Scholar

[83] View Article

[84] PubMed/NCBI

[85] Google Scholar

[ref23] 23. Atia J, McCloskey C, Shmygol AS, Rand DA, van den Berg HA, Blanks AM. Reconstruction of Cell Surface Densities of Ion Pumps, Exchangers, and Channels from mRNA Expression, Conductance Kinetics, Whole-Cell Calcium, and Current-Clamp Voltage Recordings, with an Application to Human Uterine Smooth Muscle Cells. PLoS Comput Biol. 2016;12(4):e1004828. pmid:27105427
View Article
PubMed/NCBI
Google Scholar

[87] View Article

[88] PubMed/NCBI

[89] Google Scholar

[ref24] 24. Bakker PCAM, Van Rijswijk S, van Geijn HP. Uterine activity monitoring during labor. J Perinat Med. 2007;35(6):468–477. pmid:18052832
View Article
PubMed/NCBI
Google Scholar

[91] View Article

[92] PubMed/NCBI

[93] Google Scholar

[ref25] 25. Wilde DW, Marshall JM. Effects of tetraethylammonium and 4-aminopyridine on the plateau potential of circular myometrium from the pregnant rat. Biol Reprod. 1988;38(4):836–845.
View Article
Google Scholar

[95] View Article

[96] Google Scholar

[ref26] 26. Meisheri KD, McNeill JH, Marshall JM. Effect of isoproterenol on the isolated pregnant rat myometrium. Eur J Pharmacol. 1979;60(1):1–6. pmid:520410
View Article
PubMed/NCBI
Google Scholar

[98] View Article

[99] PubMed/NCBI

[100] Google Scholar

[ref27] 27. Luo CH, Rudy Y. A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circ Res. 1991;68(6):1501–1526. pmid:1709839
View Article
PubMed/NCBI
Google Scholar

[102] View Article

[103] PubMed/NCBI

[104] Google Scholar

[ref28] 28. Luo CH, Rudy Y. A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes. Circ Res. 1994;74(6):1071–1096. pmid:7514509
View Article
PubMed/NCBI
Google Scholar

[106] View Article

[107] PubMed/NCBI

[108] Google Scholar

[ref29] 29. Colli Franzone P, Pavarino LF, Taccardi B. Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models. Math Biosci. 2005;197(1):35–66. pmid:16009380
View Article
PubMed/NCBI
Google Scholar

[110] View Article

[111] PubMed/NCBI

[112] Google Scholar

[ref30] 30. Grandi E, Pasqualini FS, Bers DM. A novel computational model of the human ventricular action potential and Ca transient. J Mol Cell Cardiol. 2010;48(1):112–121. pmid:19835882
View Article
PubMed/NCBI
Google Scholar

[114] View Article

[115] PubMed/NCBI

[116] Google Scholar

[ref31] 31. Clayton RH, Bernus O, Cherry EM, Dierckx H, Fenton FH, Mirabella L, et al. Models of cardiac tissue electrophysiology: progress, challenges and open questions. Prog Biophys Mol Biol. 2011;104(1-3):22–48. pmid:20553746
View Article
PubMed/NCBI
Google Scholar

[118] View Article

[119] PubMed/NCBI

[120] Google Scholar

[ref32] 32. Du P, Gao J, O’Grady G, Cheng LK. A simplified biophysical cell model for gastric slow wave entrainment simulation. Annu Int Conf IEEE Eng Med Biol Soc. 2013;2013:6547–6550. pmid:24111242
View Article
PubMed/NCBI
Google Scholar

[122] View Article

[123] PubMed/NCBI

[124] Google Scholar

[ref33] 33. Rihana S, Santos J, Mondie S, Marque C. Dynamical analysis of uterine cell electrical activity model. Conf Proc IEEE Eng Med Biol Soc. 2006;2006:4179–4182. pmid:17946608
View Article
PubMed/NCBI
Google Scholar

[126] View Article

[127] PubMed/NCBI

[128] Google Scholar

[ref34] 34. Laforet J, Rabotti C, Terrien J, Mischi M, Marque C. Toward a multiscale model of the uterine electrical activity. IEEE Trans Biomed Eng. 2011;58(12):3487–3490. pmid:21968708
View Article
PubMed/NCBI
Google Scholar

[130] View Article

[131] PubMed/NCBI

[132] Google Scholar

[ref35] 35. Malik M, Roh M, England SK. Uterine contractions in rodent models and humans. Acta Physiol (Oxf). 2021;231(4):e13607. pmid:33337577
View Article
PubMed/NCBI
Google Scholar

[134] View Article

[135] PubMed/NCBI

[136] Google Scholar

[ref36] 36. Matsuki K, Takemoto M, Suzuki Y, Yamamura H, Ohya S, Takeshima H, et al. Ryanodine receptor type 3 does not contribute to contractions in the mouse myometrium regardless of pregnancy. Pflugers Arch. 2017;469(2):313–326. pmid:27866274
View Article
PubMed/NCBI
Google Scholar

[138] View Article

[139] PubMed/NCBI

[140] Google Scholar

[ref37] 37. Testrow CP, Holden AV, Shmygol A, Zhang H. A computational model of excitation and contraction in uterine myocytes from the pregnant rat. Sci Rep. 2018;8(1):9159. pmid:29904075
View Article
PubMed/NCBI
Google Scholar

[142] View Article

[143] PubMed/NCBI

[144] Google Scholar

Figures

Abstract

Author summary

Introduction

Methods

Reduction of order methodology

Summary of the Tong uSMC model

Standard Simulation Protocol (SSP)

Model reduction

Elimination of channels based on conductance

Eliminating K+ channel ‘a’.

Eliminating K+ channel ‘2’.

Other K+ mechanisms.

Steady-state approximation

The K+ BK channel.

The T-Type Ca2+ channel.

The L-Type Ca2+ channel.

The voltage-gated Na+ channel.

Results

RTM performance in reproducing experiments

The SSP.

Plateaus, bursting and sustained APs.

Computational performance

Discussion

Limitations

Modeling implications: Realism and relevance

Time scales and performance

Conclusion

Acknowledgments

References

Eliminating K⁺ channel ‘a’.

Eliminating K⁺ channel ‘2’.

Other K⁺ mechanisms.

The K⁺ BK channel.

The T-Type Ca²⁺ channel.

The L-Type Ca²⁺ channel.

The voltage-gated Na⁺ channel.