Computational physiology of uterine smooth muscle

Abstract

Pregnancy can be accompanied by serious health risks to mother and child, such as pre-eclampsia, premature birth and postpartum haemorrhage. Understanding of the normal physiology of uterine function is essential to an improved management of such risks. Here we focus on the physiology of the smooth muscle fibres which make up the bulk of the uterine wall and which generate the forceful contractions that accompany parturition. We survey computational methods that integrate mathematical modelling with data analysis and thereby aid the discovery of new therapeutic targets that, according to clinical needs, can be manipulated to either stop contractions or cause the uterine wall muscle to become active.

Introduction

The onset of labour is a complex process that requires a host of regulatory pathways working in concert. In particular, the timing of delivery is critical in relation to the state of development of the foetus. A human pregnancy is defined as reaching early term at 37 weeks, whereas full term occurs at 39 weeks and late term by 41 weeks; the pregnancy is considered overdue by 42 weeks. ¹

The foetus develops throughout gestation and gradually acquires the ability to survive outside the uterine environment. If all goes well, the foetal lungs are well developed at term, which will allow survival outside of the womb with no or minimal intervention. On the other hand, delivery before 37 weeks is deemed preterm and associated with serious morbidity and mortality; it is the leading cause of death among newborns and the second highest among children under 5 years old. ² Moreover, preterm labour presents a major risk factor for developmental delay and neurological impairment. ²

To treat preterm labour, the clinician seeks to disrupt the processes that lead to uterine contractions. This can be accomplished through the reduction of excitability of the smooth muscle cells that make up the uterine muscle, or else by disrupting the coordination between the contractions in individual muscle fibres, so that the global synchronicity of muscle activity is diminished.

Reduction of cellular excitability is effective since the contraction of a smooth muscle cell is coupled to, and governed by, the electrical activity at the level of the cell membrane: quenching the latter prevents the former. ³ Disruption of electrical communication between muscle fibres may still allow partial contractions, but lack of coordination at the whole-organ level means that these partial contractions happen at different points in time, greatly diminishing the ability to generate the levels of hydrostatic pressure inside the uterine cavity that are required to expel the foetus. Any agent that has such effects is called a tocolytic, a term whose Greek etymology suggests ‘dissolution of labour’. An agent with the opposite effect, that of inducing labour, is a tocotropic. ⁴

Electrical currents through the cell membrane are conducted by entities known as ion channels. ⁵ Such ion channels represent a target of choice for pharmaceutical disruption, blockade or stimulation. However, these ion channels work together in a complex manner and, moreover, exhibit interactions with various intracellular processes as well as with processes that transcend the cellular level. This complexity necessitates a comprehensive characterisation of the dynamics of this ensemble of processes, if we are to understand in full the functional consequences of any pharmaceutical intervention. This review will survey an integrated approach that brings together acquisition of data by a variety of experimental methods and several mathematical and computational techniques.

Structure of the myometrium

The myometrium is the middle layer of the uterine wall, with the inner layer being the endometrium, which lines the uterine cavity, and the outer layer being the perimetrium. ³ As the name indicates, the myometrium consists of muscle cells (myo- referring to muscle and mētra meaning womb). ⁴ The uterine muscle cells are smooth muscle cells (uterine smooth muscle cells (USMCs); Figure 1) and accordingly have the typical characteristics associated with smooth muscle cells, namely, slow, but strong and prolonged, contraction over a wide range of movement, which is precisely what is required of the myometrium during delivery. ⁴

Figure 1.

The non-pregnant human uterus: gross morphology and anatomy, with the inset showing histological section (haematoxylin and eosin stain) of myometrial tissue.

The USMCs connect to form fibrous structures, which in turn are organised into small bundles, fasciculi. ⁶ As shown in Figure 2, in silico reconstruction of the tissue microstructure indicates that these fasciculi tend to encircle the uterine cavity in the deeper layers (i.e. closest to the endometrial lining) of the body and fundus of the uterus (i.e. the portion distal from the cervix), whereas the superficial layers exhibit longitudinally orientated fibres (green in Figure 2(c)) as well as transversal fibres (blue in Figure 2(c)); the most superficial fibres (to a depth of about 5 mm) are organised in sheets rather than bundles. Nearer to the cervical canal, the deep fibres form a circular ‘cuff’ (red in Figure 2(c)), whereas the superficial fibres appear to run in various directions.^7,8 The ability to map the three-dimensional microarchitecture of the myometrial tissue means that we can use simulations to work out the consequences at the whole-organ level of any phenomena that we observe (and may be able to represent in the mathematical form) at the molecular and cellular levels. Thus, in silico reconstruction of histological data forms a bridge between quite disparate levels of biological organisation.

Figure 2.

The non-pregnant human uterus. (a) Exterior aspect, with the cervical end at the bottom and the cranial extremity (fundus) at the top. (b) Computer-enhanced visualisation of the three-dimensional myometrial structure of the smooth muscle fibre network, based on digitisation of histological sections, followed by automated analysis of fibre direction. Superficial fibres tend to run from top to bottom, whereas deep fibres near the cervix run around the uterine cavity. (c) False colour rendering of an oblique slice through the digitised and reconstructed three-dimensional muscle fibre microarchitecture where colours in the diagram represent the bundle direction as follows: red representing left–right, green representing up–down and blue coming out of the page. Colour intensity is proportional to nuclear density.

Activity of the myometrium

The myometrium remains quiescent through much of pregnancy, but enters a more excitable state as the pregnancy approaches term. ⁴ These changes can be understood in terms of the resting membrane potential (RMP) across the cell membrane, which refers to a difference in voltage (electric potential energy) between the cellular interior (cytosol) and exterior (interstitial fluid). ⁹ This membrane potential undergoes rapid changes during an excitation event, in particular a depolarisation (also known as action potential (AP)), which is concomitant with an influx of calcium ions that enter the cell from the interstitium that surrounds it.

Whereas the RMP is always unstable in the sense that perturbations will induce an AP, the minimal magnitude of the perturbation required to do so varies. It is prohibitively large when the cells are ‘quiescent’, whereas a smaller perturbation suffices later when the pregnancy nears full term. In fact, this smaller perturbation can be readily effected by ionic currents induced by the ‘labour hormone’ oxytocin (a term whose etymology suggests acute birth). ⁴ The required perturbation is never infinitesimal; in other words, there is a finite, non-zero threshold for excitation. The threshold is large during the course of pregnancy and gradually lowers as the time for labour draws near. ⁴

The influx of calcium ions that accompanies an AP is absorbed by various buffering systems and intracellular stores, but, as these take some time to respond, a temporary increase in intracellular calcium occurs. In fact, the stores release calcium into the cytosol in response to this rise (as will be detailed below). Moreover, calcium operates on ion channels in the cell membrane, which increases the influx of calcium. Both these mechanisms exacerbate the calcium rise and effectively constitute a feedforward loop.

The rise in cytosolic calcium is instrumental to USMC function, as it ultimately leads to the contraction of the cell. Electrical excitation of the cell leads to contraction, which is referred to as excitation–contraction (E-C) coupling. ¹⁰ Calcium stimulates phosphorylation of the myosin light chain, which interacts with actin filaments to generate mechanical force. ¹¹

Humoral triggering of E-C coupling

Oxytocin is a peptide hormone consisting of nine amino acids, ¹² which is a major driver of the excitation (and thereby contraction) of the myometrium. It is secreted by the posterior pituitary gland, as well as by the uterus itself; there are thus two modes of blood-borne chemical communication, referred to as endocrine and paracrine, respectively. ⁴ During labour, secretion of oxytocin occurs as part of the Ferguson reflex (Figure 3), a positive feedback loop by which contraction of the uterus exerts pressure on the cervix, stimulating sensory neurones that drive secretion by the pituitary gland. ⁴ This in turn increases the force and frequency of contractions, leading to increased signalling. In this manner, the contractions become more frequent and forceful as labour progresses. ⁴

Figure 3.

The Ferguson reflex. Hypothalamic neurones (1) release the peptide hormone oxytocin into the bloodstream via the vascular bed of the posterior part of the pituitary (2). This oxytocin is carried by the bloodstream around the body (3). Smooth muscle cells in the uterine wall express oxytocin receptors (4), which mediate the positive effect of the hormone on the excitability of these cells. Uterine contractions lead to a build-up of pressure inside the womb, which is registered by stretch receptors in the cervix (5). The sensory neurone forms a synapse with a neurone of the ascending pathway (6) which conveys the stimulus to the hypothalamus (1), thus closing the positive feedback loop.

Oxytocin binds to the oxytocin receptor (OTR), which is a G-protein-coupled receptor. ¹³ This interaction leads to activation of phospholipase C (PLC) via a G protein.^14,15 The activated PLC cleaves the phospholipid phosphatidylinositol 4,5-bisphosphate (PIP2) into diacylglycerol (DAG) and inositol triphosphate (IP₃); whereas DAG remains bound to the plasma membrane, IP₃ is water-soluble and diffuses away from the plasma membrane. ¹⁶ Both of these products are major second messengers: DAG modulates conductances such as the BK potassium channel, ¹⁷ and IP₃ stimulates calcium release from the cell’s internal stores. These calcium ions operate on an anionic conductance in the plasma membrane, which moderately depolarises the plasma membrane. This initial depolarisation leads to the opening of the L-type voltage-gated calcium channel, causing a large calcium influx and further depolarisation resulting in an AP. ¹⁸ These events in turn engage the BK channel,^19–22 which opens and thus heralds repolarisation of the cell. The pathway is summarised in Figure 4.

Figure 4.

The phospholipase-C signalling cascade, through which the electrical excitation of a uterine muscle cell is modulated. The first messenger oxytocin (1) binds to the oxytocin receptor (2), which activates PLC via the activated alpha subunit of the G-protein (3); PLC then cleaves PIP₂ into DAG and IP₃ (4); DAG activates PKC (5), which leads to closure of the BK potassium channel (6), whereas IP₃ binds its receptor in the membrane of the sarcoplasmic reticulum (7), where the majority of intracellular calcium is stored. This interaction leads to a flow of calcium from these stores into the cytoplasm (8), which activates an anionic conductance (9). The changes in potassium and anion conductance alter the membrane potential, thus activating voltage-gated calcium channels, through which extracellular calcium can now enter the cell (10). Both the depolarisation and the rise in cytosolic calcium raise the open probability of the BK channel (11). Active mechanisms transport calcium ions into the sarcoplasmic reticulum (12).

We thus observe positive feedback loops at various levels of biological organisation: at the molecular/cellular level, as a calcium rise reinforces itself at the cellular level through calcium-mediated calcium release and calcium acting on ion channels, and also at the whole-organism level, via the Ferguson reflex. In either case, the fundamental dynamic character of a positive feedback loop is to elicit a rapid all-or-none type response in a system. This is appropriate as the uterus, geared for quiescence as the child develops, suddenly needs to commit to full activation at all levels when the time comes for the child to be born.

Ion channels and the myometrial conductome

The ‘conductome’ of a cell, or its ‘conductance repertoire’, is the term given to the collection of protein species that permit the movement of charged particles across the membrane. ²³ The cell membrane itself is phospholipid bilayer, which is highly impermeable to charged particles. Embedded within the plasma membrane are ion channels, pores that are selectively permeable to specific ions and typically composed of several proteins. ⁵

The numbers of ion channels present at the membrane and their gating parameters (such as ligand presence) can be regulated by the cell, allowing it to fine-tune ionic fluxes across the membrane. Many different ion channels exist, and the composition of the conductance repertoire varies significantly, depending upon cell type and function. ⁵

Ion channels are broadly categorised in terms of the main ion carried through the channel’s pore as well as its gating variables, that is, the physicochemical factors that affect its conductivity. ⁹ Ion channels generally allow the passive movement of ions down their electrochemical gradient when open. The cation channels are selective, allowing calcium (Ca²⁺), potassium (K⁺) or sodium (Na⁺), while the anion channels allow non-specific anion flux; the chloride ion (Cl⁻) is the most abundant anion. Besides ion channels, which allow ion fluxes between the intracellular and extracellular compartments, there also occur hemichannels that form the so-called gap junctions; these are cell-to-cell contacts that permit electrical continuity of the cytoplasm and cytoskeleton between adjacent cells, and thereby enable intercellular signalling. ²⁴ This type of communication is essential to the propagation of electrical impulses through excitable tissues such as the myometrium.^25,26

Changes occurring within the myometrium through gestation

The myometrial conductome undergoes numerous changes throughout pregnancy which correlate with, first, quiescence and, subsequently, global synchronisation. ⁴ Complete lack of excitation will prevent contraction altogether, while impaired electrical signalling will ultimately lead to disordered contraction, prolonging delivery, which carries with it increased foetal and maternal risk of complications. ⁴ Conversely, excitability of the myometrium prior to term may increase the likelihood of preterm delivery. ⁴ Even at term, if the myometrium begins to contract before the cervix if fully ripened, labour may be prolonged, increasing the risks of maternal and foetal complications. ⁴

Throughout gestation, the uterus maintains quiescence by high expression of the BK family of potassium channels. The potassium current carried by this channel acts to hyperpolarise the membrane, leading to a reduced likelihood of stimuli being able to attain the threshold for opening the L-type calcium channel and triggering an AP. ²⁷

In the event that a ‘nominally quiescent’ cell does become excited, the propagation of depolarisation events to nearby cells is prevented by the low abundance of connexins, which make up the hemichannels, throughout much of gestation. ²⁸ Connexin 43 expression is increased in preparation for labour, while the expression of connexin 26, another hemichannel component, peaks before delivery and then falls. By altering the composition of gap junctions, the spread of the electrical excitation through the myometrium is modulated. ²⁶ Indeed, numerous changes to the myometrium occur up until delivery, with significant differences being observed in the transcriptomes of gestationally matched myometrial biopsies. ²⁹

Ion channels and currents across the cell membrane

The quantitative datasets produced by patch-clamp electrophysiological experimentation lend themselves well to mathematical modelling. An individual channel may be considered to be in an open or closed state, which determines whether it is able to conduct a current. The magnitude of the current that will flow through a single channel in the open state is determined by its unitary conductance ${(\mathrm{g}}_{\mathrm{i}})$ ; wider pores in the membrane will allow a more substantial ion flux than the narrower ones, and the mechanisms by which channels select for specific ions contribute to electrical resistance. As these factors are different for each particular type of channel, each accordingly has a different unitary conductance.^5,9

The challenge for a mathematical model is to determine the likelihood that a channel of species $\mathrm{i}$ will be open at a particular time ${(\mathrm{P}}_{\mathrm{i}}\left(\mathrm{t}\right))$ and also to determine how quickly the channel will transition between open and closed states. The factors that cause a channel to move between an open and a closed state are termed gating variables.^5,9 In constructing mathematical models, these variables may be measured empirically (such as membrane potential, ${\mathrm{E}}_{\mathrm{m}}$ ) or estimated by means of a calculation. Although gating variables are naturally associated with an interval scale, limitations in measurement techniques sometimes constrain one to work a weaker scale.

When numerous copies of a channel type are considered together, for example, for an entire cell, the open probability ${\mathrm{P}}_{\mathrm{i}}\left(\mathrm{t}\right)$ can be taken to represent the fraction of the total channel population that is in the open state; this is a simple corollary of the law of large numbers. Thus, if the membrane density of channels of the species can be determined $\left({\mathrm{\kappa }}_{\mathrm{i}}\right)$ , the total current that will be carried ${(\mathrm{I}}_{\mathrm{i}})$ can be approximated as follows

$${\mathrm{I}}_{\mathrm{i}}={\mathrm{\kappa }}_{\mathrm{i}}{\mathrm{g}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{i}}\left({\mathrm{E}}_{\mathrm{m}}-{\mathrm{E}}_{\mathrm{rev},\mathrm{i}}\right)$$

The reversal potential ${\mathrm{E}}_{\mathrm{rev},\mathrm{i}}$ may be calculated from the Nernst equation

$${\mathrm{E}}_{\mathrm{rev},\mathrm{i}}=\frac{\mathrm{RT}}{{\mathrm{z}}_{\mathrm{i}}\mathrm{F}}\ln \frac{\mathrm{intracellular}\mathrm{concentration}\mathrm{of}\mathrm{ion}\mathrm{carried}\mathrm{by}\mathrm{i}}{\mathrm{extracellular}\mathrm{concentration}\mathrm{of}\mathrm{ion}\mathrm{carried}\mathrm{by}\mathrm{i}}$$

where $\mathrm{R}$ is the gas constant, $\mathrm{T}$ the absolute temperature, ${\mathrm{z}}_{\mathrm{i}}$ the valency of the ion species conducted by channel type $\mathrm{i}$ and $\mathrm{F}$ is Faraday’s constant. ⁹ This formula represents the thermodynamic equilibrium between the electrical forces experienced by an ion and its random Brownian motion leading to a net movement of the ion down its concentration gradient. When these two processes oppose each other across a selectively permeable membrane, but are in equilibrium, there is no net movement of the ion. The value of ${\mathrm{E}}_{\mathrm{m}}$ at which this happens is ${\mathrm{E}}_{\mathrm{rev},\mathrm{i}}$ . ⁹ The RMP can be intuitively thought of as a weighted compromise between the ${\mathrm{E}}_{\mathrm{rev},\mathrm{i}}$ values for the various ionic species $\mathrm{i}$ .

The cytosolic and interstitial fluids contain electrolytes which endow these media with rapid Debye equilibration, ⁹ which compels any charge imbalance to accumulate atomistically close to cell membrane; accordingly, the latter acts as a capacitance ${\mathrm{C}}_{\mathrm{m}}$ . If we treat Debye equilibration as infinitely rapid, the potential difference ${\mathrm{E}}_{\mathrm{m}}$ between cytosol and interstitial fluid follows immediately from charge conservation

$${\mathrm{E}}_{\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{C}}_{\mathrm{m}}^{-1}\sum _{\mathrm{i}}\underset{-\infty }{\overset{\mathrm{t}}{\int }}{\mathrm{I}}_{\mathrm{i}}\left(\tau \right)\mathrm{d}\tau$$

which makes clear that we will have obtained a mathematically closed system as soon as ${\mathrm{I}}_{\mathrm{i}}$ has been described for all channel species $\mathrm{i}$ extant within the conductome.

The pioneering work by Hodgkin and Huxley

Hodgkin et al.^30–34 modelled the AP in a seminal series of papers. Not only did they develop the technique of patch-clamp electrophysiology, but they also used it to derive mathematical models that describe ionic currents across the membrane in an approach that remains relevant nearly 70 years later. This feat is made all the more impressive by the fact that the existence of the ion channels was not confirmed until years later.

In the Hodgkin–Huxley (HH) model, an ion channel is envisaged as a pore in the membrane containing a series of barriers, each of which can bar the flow of ions by being in the closed state, and all of which must be open to allow the pore to be conducting. This conceptual model is shown schematically in Figure 5. The assumption of statistical independence between these movable barriers (‘gates’) translates into the proposition that the probability of the channel being in a conducting state can be written as the product of the corresponding probabilities for the individual gates. In particular, when all the gates are of the same type, one would simply take a power of the individual gate probability to give the overall pore probability (the exponent being the number of gates that are arranged in series within the pore of the channel).

Figure 5.

Representation of the Loligo sodium conductance model as described by Hodgkin and Huxley. The model consisted of four independent gates: 3 m-gates (blue) and 1 h-gate (pink), all of which must be open in order to allow current to flow. At the resting potential, the channel is unlikely to be open, as one or more m-gates will be closed. Following an initial depolarisation of the membrane, the m-gates are opened, allowing sodium ion flux. After sustained depolarisation, the slower h-gate closes, inactivating the current.

The genius of Hodgkin and Huxley consisted in translating a series of voltage-clamp measurements into a dynamical model of these gates. In particular, the probability for each gate is assumed to follow linear relaxation behaviour towards a final value, with the key proviso that both this final value and the relaxation rate are taken to be functions of ${\mathrm{E}}_{\mathrm{m}}$ . It is this functional dependence that makes ${\mathrm{E}}_{\mathrm{m}}$ the gating variable (other examples of gating variables are, for example, concentrations of intracellular or extracellular ions or molecules, which might be metabolites such as ATP or secondary messengers such as IP₃; cf. Figure 4). The model’s non-linear properties, notably its ability to reproduce the AP, derive from the presence of voltage dependence in the relaxation rates and instantaneous equilibrium values of the gating variables.

Empirical state transition graphs

The state transition kinetics for a single gate in the HH model is particularly simple

$$\mathrm{C}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC}}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO}}\end{array}\mathrm{O}$$

where $\mathrm{C}$ is the closed state, $\mathrm{O}$ the open state, m_oc ${\mathrm{\mu }}_{\mathrm{CO}}$ the transition rate (in the sense of a hazard function on the exponentially distributed lifetime of the closed state) from closed to open and m_oc ${\mathrm{\mu }}_{\mathrm{OC}}$ is the corresponding rate from open to closed. The stroke of brilliance was not only to hypothesise that these rates both depend on the reigning membrane potential ${\mathrm{E}}_{\mathrm{m}}$ , but also to devise a technique to reconstruct empirical formulas for these dependencies. The above diagram can be viewed as a state transition graph with two nodes (open and closed) and two directed edges (it is thus, technically speaking, a ‘directed graph’ or digraph but we will just use the friendlier term graph). The state transition graph of the entire channel is then defined as the Cartesian product over the graphs for each of the gates, where the assumption of statistical independence allows us to carry over the transition rates from the factor graphs in the obvious way; Figure 6 illustrates this procedure for the product of two simple gates. Three or more gates can be multiplied together in the same manner, although the resulting hypercube graph is perhaps less readily visualised as the numbers become larger.

Figure 6.

Cartesian multiplication of two gate graphs. The box symbol on the left represents the operation of graph multiplication. The vertices in the product graph (on the right) contain the open/closed information regarding the gates in the two factor graphs on the left. The rate labels carry over in the obvious manner, as a result of statistical independence.

The vertices in the state transition graph correspond to minima in the so-called ‘energy landscape’ of the ion channel. The regular hypercube arrangement of the HH model is seldom encountered in nature. Individual gates often have more than two states and are not necessarily statistically independent. ⁵ Two representative examples of empirically determined state transition graphs are shown in Figure 7. In addition to closed and open states, there are inhibited states as well, closed states with long persistence times or exit rates that depend on special gating.

Figure 7.

Examples of empirical state transition graphs. Top: Kv4.3 channel; bottom: hERG channel.

Several remarks are in order since the issue appears to invite numerous common misunderstandings. First, the HH model channels are often discussed as if they were antonymic to ‘Markov chain models’, but this is a distinction without a difference. Since state transition graphs are isomorphic to continuous-time Markov chains (i.e. there is a unique way of going from an instance of the former to an instance of the latter, and vice versa), both have equal call to be regarded as such.

Second, models based on empirical state transition graphs are often erroneously believed to be more complex than HH models; the truth is that models such as those typically encountered in the literature (e.g. Figure 7) have rather simpler graphs than the HH model depicted in Figure 5, as may be attested by anyone who has attempted to sketch on a flat page, say, the four-dimensional hypercube defined by the graph product for the classic HH Loligo giant axion sodium channel

$${\mathrm{C}}_{\mathrm{m}}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC},\mathrm{m}}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO},\mathrm{m}}\end{array}{\mathrm{O}}_{\mathrm{m}}\square {\mathrm{C}}_{\mathrm{h},1}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC},\mathrm{h},1}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO},\mathrm{h},1}\end{array}{\mathrm{O}}_{\mathrm{h},1}\square {\mathrm{C}}_{\mathrm{h},2}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC},\mathrm{h},2}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO},\mathrm{h},2}\end{array}{\mathrm{O}}_{\mathrm{h},2}\square {\mathrm{C}}_{\mathrm{h},3}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC},\mathrm{h},3}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO},\mathrm{h},3}\end{array}{\mathrm{O}}_{\mathrm{h},3}$$

where the symbol □ represents the Cartesian graph multiplication operator; see Figure 6. It should be admitted, however, that the energy landscape of a proteinaceous entity as large as an ion channel will have hundreds if not thousands of (mostly extremely shallow) local minima; an ultrahigh-resolution empirical state transition graph could accordingly grow to enormous sizes, as the scope and complexity of the graph is a function of the accuracy with which the energy landscape of a given ion channel has been probed.

Third, HH models and models based on empirical state transition graphs are often contrasted in terms of ‘macroscopic’ currents versus ‘microscopic’ currents. The ‘macro’ currents are currents that are in reality carried by two or more distinct ion channel species, but represented in the model as though it were carried by a single type of pore (or at any rate, by fewer pore types than actually present), whereas ‘micro’ currents are purported to each correspond to a single ion channel species in the conductome. For example, the large complement of potassium channels in an USMC (Figure 8) may be represented as only a handful so-called ‘pooled conductances’.

Figure 8.

Schematic representation of the uterine smooth muscle cell, with its extended conductome, that is, all electrogenic entities predicted on the basis of tissue-level transcriptomics. Potassium channels are indicated in red, calcium channels in blue and chloride channels in yellow.

In brief, empirical state transition models are definitely intended to correspond to individual channel species (up to the caveat related to the accuracy with which the energy landscape is known), whereas the HH-type state transition graph is vanishingly unlikely to be the empirical graph of any existing ion channel, given its exceptionally regular hypercube topology (even if this is not logically impossible, and it would be an interesting challenge to design a proteinaceous ion channel that obeys HH rules). At the same time, the HH model is a time-honoured routine tool to approximate ‘macroscopic’ currents. Finally, the empirical state transition models are often loosely said to be ‘superior’ – in view of the preceding three remarks one may appreciate both why this impression easily arises and why it is misleading at best.

Stochastic versus deterministic currents at the cellular level

While there is no fundamental logical distinction between HH models and those predicated on empirical state transition graphs, there do exist important differences between the ways in which either can be treated from a computational point of view. To appreciate this point, let us consider the simplest non-trivial state transition graph

$$\mathrm{A}\stackrel{{\mathrm{\mu }}_{\mathrm{AB}}}{\to }\mathrm{B}$$

where mAB is a constant. Individual systems subject to this kinetics behave stochastically. In a numerical context, continuous-time Monte Carlo simulation can be employed. ³⁵ If the system (also referred to in this context as particle) starts out in state B, there is nothing to be said and the particle just stays there forever. If the particle is in state A to begin with, at time zero, it transitions to state B at a point in time $\mathrm{T}$ units after time zero, where $P[\mathrm{T}>\mathrm{t}]=\exp \{-{\mathrm{\mu }}_{\mathrm{AB}}\mathrm{t}\}$ (i.e. $\mathrm{T}=-{\mathrm{\mu }}_{\mathrm{AB}}^{-1}\mathrm{lnU}$ , where $P[\mathrm{U}>\mathrm{u}]=\mathrm{u}$ , a useful thing to know if drawings from the standard uniform distribution are all one has at one’s disposal), and it will remain in state B henceforth. If one has ${\mathrm{n}}_{\mathrm{A}}\left(0\right)$ particles in state A at time zero, the number of particles at time $\mathrm{t}$ is subject to the following approximation

$${\mathrm{n}}_{\mathrm{A}}\left(\mathrm{t}\right)\approx {\mathrm{n}}_{\mathrm{A}}\left(0\right)\exp \left\{-{\mathrm{\mu }}_{\mathrm{AB}}\mathrm{t}\right\}$$

which becomes more accurate as ${\mathrm{n}}_{\mathrm{A}}(0)\to \infty$ by the law of large numbers. Mathematicians make the idea of accuracy more precise by describing probabilistic bounds on the error of the approximation as a function of ${\mathrm{n}}_{\mathrm{A}}(0)$ .

The number of ion channels present in a given cell’s conductome will never be infinite, so the question becomes whether this number is sufficiently large to warrant the use of this ‘large numbers’ approximation (also known as the approximation ‘in expectation’ or the ‘mean field’ approximation). In practice, the criterion ought to be whether the approximation incurred is negligible compared to other sources of error that affect the quantitative prediction.

The basics of the mean field approximation

Let us consider the basic gate

$$\mathrm{C}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC}}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO}}\end{array}\mathrm{O}$$

We have the following expression for the open probability ${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)$

$${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)={\overline{\mathrm{P}}}_{\mathrm{O}}+\left({\mathrm{P}}_{\mathrm{O}}\left(\mathrm{0}\right)-{\overline{\mathrm{P}}}_{\mathrm{O}}\right)\exp \left\{-\left({\mathrm{\mu }}_{\mathrm{OC}}+{\mathrm{\mu }}_{\mathrm{CO}}\right)\mathrm{t}\right\}$$

where ${\overline{\mathrm{P}}}_{\mathrm{O}}={(1+{\mathrm{\mu }}_{\mathrm{OC}}/{\mathrm{\mu }}_{\mathrm{CO}})}^{-1}$ . The mean field approximation amounts to positing that, out of the $\mathrm{N}$ instances of this system that are present, a number ${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)\mathrm{N}$ is open. The variance about this expectation is ${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)(1-{\mathrm{P}}_{\mathrm{O}}(\mathrm{t}\left)\right)\mathrm{N}$ . If we divide the standard deviation by the expectation, we obtain the coefficient of variation, which gives an indication of the relative error. This quantity works out as $\sqrt{(1/{\mathrm{P}}_{\mathrm{O}}(\mathrm{t})-1)/\mathrm{N}}$ , which can be made arbitrarily small by taking $\mathrm{N}$ large enough – a fairly elementary observation which nonetheless forms the foundation for the mean field approach.

For the ‘product system’ in Figure 6, we have

$$\begin{array}{c}{\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)=\left({\overline{\mathrm{P}}}_{\mathrm{O},1}+\left({\mathrm{P}}_{\mathrm{O},1}\left(\mathrm{0}\right)-{\overline{\mathrm{P}}}_{\mathrm{O},1}\right)\exp \left\{-\left({\mathrm{\mu }}_{\mathrm{OC},1}+{\mathrm{\mu }}_{\mathrm{CO},1}\right)\mathrm{t}\right\}\right)\\ \times \left({\overline{\mathrm{P}}}_{\mathrm{O},2}+\left({\mathrm{P}}_{\mathrm{O},2}\left(\mathrm{0}\right)-{\overline{\mathrm{P}}}_{\mathrm{O},2}\right)\exp \left\{-\left({\mathrm{\mu }}_{\mathrm{OC},2}+{\mathrm{\mu }}_{\mathrm{CO},2}\right)\mathrm{t}\right\}\right)\end{array}$$

where ${\overline{\mathrm{P}}}_{\mathrm{O},i}={(\mathrm{1}+{\mathrm{\mu }}_{\mathrm{OC},i}/{\mathrm{\mu }}_{\mathrm{CO},i})}^{-1}$ for $\mathrm{i}=1,2$ . In the special case where the rates are equal ( ${\mathrm{\mu }}_{\mathrm{OC},1}={\mathrm{\mu }}_{\mathrm{OC},2}={\mathrm{\mu }}_{\mathrm{OC}}$ and ${\mathrm{\mu }}_{\mathrm{CO},1}={\mathrm{\mu }}_{\mathrm{CO},2}={\mathrm{\mu }}_{\mathrm{CO}}$ ), this reduces to an expression that is even nicer

$${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)={\left({\overline{\mathrm{P}}}_{\mathrm{O}}+\left({\mathrm{P}}_{\mathrm{O}}\left(0\right)-{\overline{\mathrm{P}}}_{\mathrm{O}}\right)\exp \left\{-\left({\mathrm{\mu }}_{\mathrm{OC}}+{\mathrm{\mu }}_{\mathrm{CO}}\right)\mathrm{t}\right\}\right)}^2$$

Such neat factorisation is the hallmark of HH models, but it does not generically obtain for any given empirical state transition graph. If the graph has $\mathrm{n}$ vertices, the probabilities for each of these states obey a non-autonomous linear differential equation

$$\frac{\mathrm{d}}{\mathrm{dt}}\mathbf{p}=\mathbf{M}\left(\mathbf{u}\left(\mathrm{t}\right)\right)\cdot \mathbf{p}$$

where $\mathbf{p}\in {[0,1]}^{\mathrm{n}}$ is a vector for which the ith element denotes the probability of the system being in the ith state (graph vertex), $\mathbf{u}\in R^{\mathrm{m}}$ is a gating input (of some dimension $\mathrm{m}$ ; when there is a single gating variable, such as membrane potential, we have $\mathrm{m}=1$ ; if the channel is, in addition, calcium gated, we would have $\mathrm{m}=2$ and so on, according to whether the channel is sensitive to several other gating inputs) and $\mathbf{M}$ is an $\mathrm{n}\times \mathrm{n}$ matrix which depends on $\mathbf{u}$ , generally in a non-linear fashion, whereas $\mathbf{u}$ depends on time. Since the channel has to be in some state and can only be in one state at a time, it follows that the elements of $\mathbf{p}$ sum to $1$ . As we saw above, the dynamics of ${\mathrm{E}}_{\mathrm{m}}$ depends itself on the currents that flow, thus ‘closing the loop of causality’.

The open probability is given by the expression

$${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)=\mathbf{c}\cdot \mathbf{p}\left(\mathrm{t}\right)$$

where $\mathbf{c}$ is a $1\times \mathrm{n}$ matrix storing the unitary conductances associated with each of the $\mathrm{n}$ states. Usually, $\mathbf{c}$ can be factored into a scalar which is the unitary conductance of the channel, and a row vector consisting of 0s and 1s according to whether the corresponding state is open.

If we fix $\mathbf{u}(\mathrm{t})$ at some constant value ${\mathbf{u}}^{*}$ (i.e. we ‘clamp’ the gating input in the parlance of the electrophysiologist), the above differential equation becomes autonomous. We have the standard solution for a constant linear system ³⁶

$${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)={\overline{\mathrm{P}}}_{\mathrm{O}}^{*}+\sum _{\mathrm{k}=1}^{\mathrm{n}-1}{\mathrm{\kappa }}_{\mathrm{k}}^{*}\exp \left\{-{\lambda }_{\mathrm{k}}^{*}\mathrm{t}\right\}$$

where $0<{\lambda }_1^{*}<{\lambda }_2^{*}<\cdots <{\lambda }_{\mathrm{n}-1}^{*}$ by a Lyapunov-type argument on closed compartmental systems (which also accounts for the celebrated H-theorem) ³⁶ – we assume here for the sake of simplicity that the eigenvalues of $\mathbf{M}({\mathbf{u}}^{*})$ are all distinct, as is generically warranted. The asterisk indicates dependence on ${\mathbf{u}}^{*}$ , that is, ${\overline{\mathrm{P}}}_{\mathrm{O}}^{*}={\overline{\mathrm{P}}}_{\mathrm{O}}({\mathbf{u}}^{*})$ . In intuitive terms, this result says that ${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)$ will relax towards ${\overline{\mathrm{P}}}_{\mathrm{O}}^{*}$ , apart from no more than $\mathrm{n}-1$ ‘wobbles’. In this sense, the behaviour of the gating-clamped system can be regarded as fairly simple. Put differently, much of the dynamic interest of a conductance really derives from its gating dependency.

If the modes of the system (time scales $\sim 1/{\lambda }_{\mathrm{k}}^{*}$ ) are sufficiently fast, we may use the approximation ${\mathrm{P}}_{\mathrm{O}}\left(\mathrm{t}\right)\approx {\overline{\mathrm{P}}}_{\mathrm{O}}(\mathbf{u}(\mathrm{t}))$ – the mean field conductance is ‘slaved’ to the gating input. Alternatively, using rate quantities of the form

$${\mathrm{M}}_{\mathrm{ij}}^{-1}\frac{\partial {\mathrm{M}}_{\mathrm{ij}}}{\partial {\mathrm{u}}_{\ell }}\frac{{\mathrm{du}}_{\ell }}{\mathrm{dt}}$$

as our guide, we can divide the time axis into a succession of intervals, during each of which the gating input is treated as not changing substantially, and employ a piecewise clamped approximation. That is to say, we pretend that the gating input is constant over intervals and concatenate the corresponding clamped solutions. Once more, some economy might be achieved by discarding one or more ‘fast’ modes.

Yet another option is to consider whether the analytical solution is close to being factorable, that is, whether the system can be approximated as a system of independent gates – the number of which would have to be rather less than $\mathrm{n}-1$ to make this approach worth considering.

Continuous-time Monte Carlo simulation

Whenever the state transition graph is the HH hypercube, it appears that, in practice, the mean field approximation is used. Although this has no a priori justification, habitual practice dovetails with the use of HH as a ‘macro’ model for pooled conductances. On the other hand, in the case of empirical state transition graphs, the subject at hand is a definite species of ion channel, which may be present at such low copy numbers per cell that the mean field approximation is not valid. In this case, continuous-time Monte Carlo simulation is required. We are dealing with the ‘exponentiated’ graph

$${\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}^{\square {\mathrm{N}}_{\mathrm{ion}\mathrm{channel}}}=\underset{{\mathrm{N}}_{\mathrm{ion}\mathrm{channel}}\mathrm{copies}}{\underbrace{{\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\square {\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\square \cdots \square {\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}}}$$

where ${\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}$ is the graph of a single ion channel of interest (cf. Figure 7 for examples) and ${\mathrm{N}}_{\mathrm{ion}\mathrm{channel}}$ is the copy number present on the cell. As ${\mathrm{N}}_{\mathrm{ion}\mathrm{channel}}$ grows without bound, the system ${\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}^{{\square \mathrm{N}}_{\mathrm{ion}\mathrm{channel}}}$ behaves more and more deterministically, as prescribed by the mean field approximation. The exponential power ${\mathrm{N}}_{\mathrm{ion}\mathrm{channel}}$ should certainly not be confused with the exponential in an individual HH-type channel with several gates. In the latter, the single channel is conducting only if all the gates in its pore are open; in a cellular population of $\mathrm{N}$ channels, the total current is proportional to the number of channels that are open, as we saw above.

Statistical independence greatly simplifies the Monte Carlo simulation. For instance, if we take the exponentiated system

$${\left(\mathrm{C}\begin{array}{c}{\mathrm{\mu }}_{\mathrm{OC}}\\ \leftrightarrows \\ {\mathrm{\mu }}_{\mathrm{CO}}\end{array}\mathrm{O}\right)}^{\square \mathrm{N}}$$

and consider a point in time $\mathrm{t}$ where we have ${\mathrm{n}}_{\mathrm{A}}(\mathrm{t})$ particles in state A (and hence ${\mathrm{n}}_{\mathrm{B}}(\mathrm{t})=\mathrm{N}-{\mathrm{n}}_{\mathrm{A}}(\mathrm{t})$ in state B), elementary calculations exploiting statistical independence allow us to show that

$$\ln \mathbb{P}\left[\mathrm{no}\mathrm{particle}\mathrm{changes}\mathrm{state}\mathrm{between}\mathrm{t}\mathrm{and}{\mathrm{t}}^{\prime }\right]=\left(\mathrm{N}{\mathrm{\mu }}_{\mathrm{BA}}+{\mathrm{n}}_{\mathrm{A}}\left(\mathrm{t}\right)\left({\mathrm{\mu }}_{\mathrm{AB}}-{\mathrm{\mu }}_{\mathrm{BA}}\right)\right)\left(\mathrm{t}-{\mathrm{t}}^{\prime }\right)$$

for ${\mathrm{t}}^{\prime }\ge \mathrm{t}$ , and that the first particle to change state will be a state-A particle (transitioning to state B) with probability ${(1+({\mathrm{\mu }}_{\mathrm{BA}}/{\mathrm{\mu }}_{\mathrm{AB}})(\mathrm{N}/{\mathrm{n}}_{\mathrm{A}}(\mathrm{t})-1))}^{-1}$ . These results make the setting up of a stochastic simulation virtually trivial. Similarly, if we have a general exponentiated graph ${\mathrm{G}}^{\square \mathrm{N}}$ , we need only to consider a single instance of $\mathrm{G}$ , with ‘buckets’ of particles at each vertex containing the particles, and employ this type of calculation to choose, probabilistically, the next ‘buckets exchanging a particle’ event.

If we picture this system as the single gate in a (rather elementary variety of) ion channel, we may consider the current conducted, which is proportional to the fraction of the pores that are in the open state out of a per-cell population of $\mathrm{N}$ . Such a current is graphed as a function of time in Figure 9, for various values of $\mathrm{N}$ . For $\mathrm{N}=1$ , we observe the current switching abruptly between its minimum and its maximum value, as the gate transitions between the $\mathrm{C}$ and $\mathrm{O}$ states. For larger $\mathrm{N}$ , more intermediate levels of current become available in steps of $1/\mathrm{N}$ of the maximum; moreover, the noisiness gradually diminishes as the mean field approximation improves (cf. the $\sim 1/\sqrt{\mathrm{N}}$ argument alluded to above). In the limit $\mathrm{N}\to \infty$ , these levels become arbitrarily close to one another and the system resides, ‘most of the time’, closest to the level corresponding to the equilibrium open probability ${(1+{\mathrm{\mu }}_{\mathrm{OC}}/{\mathrm{\mu }}_{\mathrm{CO}})}^{-1}$ (statistical mathematics makes this intuitive but vague phrase ‘most of the time’ more precise).

Figure 9.

Current (ordinate) plotted as a function of time (abscissa). Numbers indicate $\mathrm{N}$ in the system ${(\underset{{\mathrm{\mu }}_{\mathrm{CO}}}{\overset{{\mathrm{\mu }}_{\mathrm{OC}}}{\mathrm{C}\leftrightarrows \mathrm{O}}})}^{\square \mathrm{N}}$ . Halfway through the simulation, the opening rate mco ${\mathrm{\mu }}_{\mathrm{CO}}$ is suddenly tripled to mimic a gating effect, as may be observed during voltage-clamp experiments. Abscissa is scaled to 110% of the maximum current carried.

Stochastic disruption

Far from merely being of mathematical interest, the difference between the Monte Carlo and mean field regimes could well be of crucial physiological importance. Consider a number of electrically isolated USMCs, and assume for the sake of the argument that they all have identical conductomes and receive identical stimulation, amounting to a supraliminal perturbation in the mean field. However, all have a potassium channel with high unitary conductance, present in a large enough copy number (identical for every cell) to inhibit the AP when all are open, but small enough that the mean field approximation fails and stochasticity reigns. In these circumstances, what we observe is that some cells are excited, whereas others remain quiescent, even though in the mean field all would be excited.

Let us contrast a large number of potassium channels with small unitary conductances to a small number of potassium channels with large unitary conductances (but otherwise identical gating characteristics). It can clearly be arranged that the two cases would behave identically in the mean field. However, let us suppose that the mean field approximation is only valid for the many channels/small unitary conductance case; the few channels/large unitary conductance case is then revealed to have the ability to desynchronise the myometrium, as soon as stochasticity is taken into account. This effect – stochastic disruption – may help explain why we encounter such a byzantine variety of potassium channels in a single cell type (Figure 8) and why the potassium channel BK plays a major role during the course of pregnancy but not during delivery.

Gating heterogeneity

Statistical independence loomed large in the foregoing arguments, and we would therefore do well to ask if, besides being mathematically convenient, this assumption is well warranted from a biological point of view. The harsh truth is that this is not the case. However, we find a saving grace in the fact that correlations between ion channels are generally mediated by their gating characteristics (an exception occurs when channels are capable of direct physical interaction, bypassing the intermediary of dissolved molecules or the reigning membrane potential). In other words, the channels are conditionally independent given their gating, or, as Bayesian statisticians would put it, ion channels are d-separated on the evidence provided by the gating inputs.

In terms of our notation, the situation is as follows: each channel has its (generically distinct) gating input ${\mathbf{u}}_{\ell }(\mathrm{t})$ , and, instead of the exponentiated graph, the general case involves the more intricately structured Cartesian product

$${\square }_{\ell =1}^{\mathrm{N}}{\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\left({\mathbf{u}}_{\ell }\left(\mathrm{t}\right)\right)={\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\left({\mathbf{u}}_1\left(\mathrm{t}\right)\right)\square {\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\left({\mathbf{u}}_2\left(\mathrm{t}\right)\right)\square \cdots \square {\mathrm{G}}_{\mathrm{ion}\mathrm{channel}}\left({\mathbf{u}}_{\mathrm{N}}\left(\mathrm{t}\right)\right)$$

where the index $\ell$ does not stand for channel type, but ranges over the $\mathrm{N}$ individual molecular copies of a given species of ion channel. Statistical independence only obtains when ${\mathbf{u}}_{\ell }(\mathrm{t}){\equiv }_{\ell }\mathbf{u}(\mathrm{t})$ , that is, all copies of the ion channel ‘see’ the same gating environment. This occurs when gating ligands are well mixed and present at the same concentration everywhere across the cell. Another important example of such gating uniformity is the membrane potential, which equilibrates across the cellular surface at an extremely rapid rate. ⁹

USMCs are connected by gap junctions (Figure 8). As the pregnancy reaches full term, connexins are upregulated, ³⁷ leading to a dramatic decrease of electrical resistance across these junctions. This not only permits the propagation of excitation over larger regions of the organ, but also brings the cytoplasms of neighbouring cells closer to behaving as a single conductor. In mathematical terms, this implies that the mean field approximation becomes better. In other words, the stochastic disruption effect described above is suppressed by the electrical joining up of neighbouring cells. Under these conditions, a fasciculus becomes more deterministic and unified in terms of its electrophysiological behaviour.

As shown in Figure 4, cytosolic calcium is a key intermediary in cellular excitation, its local concentration acting as an important gating variable. Inasmuch as this concentration is spatially heterogeneous even within a single cell, we have a situation where two copies of the same ion channel species, say $\ell$ and ${\ell }^{\prime }$ , can experience distinct gating inputs, that is, ${\mathbf{u}}_{\ell }(\mathrm{t})\ne {\mathbf{u}}_{{\ell }^{\prime }}\left(\mathrm{t}\right)$ .

To describe this heterogeneity, the modeller can opt for different approaches. One is to obtain data regarding these inputs at a high spatiotemporal resolution and feed these directly into the model at the level of $\{{\mathbf{u}}_1(\mathrm{t}),\dots ,{\mathbf{u}}_{\mathrm{N}}(\mathrm{t})\}$ . Another is to formulate mathematical models for these inputs. For instance, one might represent the local calcium concentration as a stochastic process. The latter can take the form of a continuous-time Markov chain with state transition graph $\mathrm{H}$ . The kinetics topology of a calcium-gated channel can then be represented by the Markov chain over the product $\mathrm{G}\Omega \mathrm{H}$ ; the rate labels on the directed edges (‘arrows’) do not carry over as shown in Figure 6, since G’s rates are dependent on the state (i.e. vertex) of $\mathrm{H}$ , and this cross-dependence has to be taken into account (to make matters even more complicated, the rates of $\mathrm{H}$ depend on the global state of the cell).

In the case where the calcium fluctuations are spatially uncorrelated across the cell, the ion channel population can be represented by the exponentiated graph ${(\mathrm{G}\Omega \mathrm{H})}^{\Omega \mathrm{N}}$ ; the latter, finally, may or may not be amenable to a mean field approximation (thus, even if local stochasticity is taken properly into account, the global, cellular behaviour may turn out to be ‘deterministic’).

This construction can also be used when two or more ion channels are in direct physical interaction. Thus, if we have two ion channels $\mathrm{G}$ and ${\mathrm{G}}^{\prime }$ , we study the Markov chain over the product $\mathrm{G}\square {\mathrm{G}}^{\prime }$ . We emphasise that this construction is necessary and proper, whether or not we ultimately work with the mean field approximation or with Monte Carlo simulation, for instance.

As a physiologically relevant application of gating heterogeneity, we consider the cytosolic calcium gradients that are created when calcium diffuses via the gap junctions from an excited cell to its quiescent neighbour. ³⁸ The cytosolic calcium concentration in the neighbouring cell has risen in the process of activation, as described in the foregoing sections, and the concentration difference across the gap junction causes a net diffusion of calcium into the inactive cell. The spindle-shaped USMC can be treated as axially symmetric. The relevant spatial parameter is therefore just the longitudinal distance from the gap junctional zone. If no other calcium-related phenomena are (as yet) happening, such as release from the sarcoplasmic reticulum (SR) or influx through calcium channels, the spatiotemporal description of the cytosolic calcium concentration in the inactive cell can be treated as a standard diffusion problem. Combining the two, one can infer ${\mathbf{u}}_{\ell }(\mathrm{t})$ , the gating input on the $\ell \mathrm{th}$ copy of the channel.

If the cell were to localise calcium-gated potassium channels (such as BK^19,20) near gap junctions, the calcium ions that enter would furnish a drive towards ${\mathrm{E}}_{\mathrm{rev},{\mathrm{Ka}}^{+}}$ , which counteracts the depolarising current conducted via the gap junction and thus effectively shutting down cell-to-cell propagation of activity (even when connexins are being expressed). By contrast, if the cell were to populate the neighbourhood of gap junctions with channels that conduct ions whose ${\mathrm{E}}_{\mathrm{rev}}$ is more positive than the RMP (e.g. Cl⁻, Na⁺, Ca²⁺), this would provide a drive that reinforces the depolarising currents transmitted via the connexins. Calcium-gated ion channels that co-localise with gap junctions could therefore serve as a kind of transistor, modulating the transfer of the exciting stimulus from one cell to the next, both in a positive and negative sense.

The transistor effect indicates that the spatial disposition of calcium-gated ion channels over the cell membrane, relative to the location of the gap junctions, should be carefully observed and taken into account in the mathematical models of individual cells. This effect may be of particular physiological importance in controlling the global spread of activation, that is, the ability of an excitation event to propagate over the entire organ.

Detailed biophysics of myometrial excitability

The USMC’s conductome comprises a vast repertoire of potential conductance entities, where the multimeric nature of ion channels contributes combinatorial richness (Figure 8). Modellers of excitable cells could attempt to tackle this diversity by pooling conductances together into ‘macro’ currents that are more readily characterised using electrophysiological methods. However, this approach has a fundamental flaw, inasmuch as a given putative drug compound generally will have distinct pharmacodynamics for each individual ion channel. This severely limits the utility of any model based on ‘macroscopic’ currents (i.e. ‘pooled conductances’), particularly in the context of evaluating or predicting the effects of putative drugs such as candidate tocolytics or tocotropics.

To overcome these limitations, Atia et al. ²³ constructed a comprehensive ‘micro’ currents model incorporating all potential conductance species that were detected by Chan et al. ²⁹ using RNAseq. This dataset was obtained from human tissue samples, and so is immediately applicable to the human USMC. The potential conductome (conductance repertoire) that could be generated through the combination of each of the constituent subunits expressed was determined. Empirical state transition kinetics for each of these conductances was gleaned from a comprehensive literature search. This exercise left only the expression levels of the ion channels as the free parameters, that is, a set

$$\left\{{\mathrm{\kappa }}_1,{\mathrm{\kappa }}_2,\dots ,{\mathrm{\kappa }}_{\mathrm{n}}\right\}$$

where $\mathrm{n}$ is about 30, to be estimated from further experimental data. The customary approach is to take observations, for example, a series of measurements of ${\mathrm{E}}_{\mathrm{m}}\left(\mathrm{t}\right)$ over an interval of time, that is, $\mathrm{t}\in [{\mathrm{t}}_1,{\mathrm{t}}_2]$ where ${\mathrm{t}}_1$ and ${\mathrm{t}}_2$ are the start and termination times of the period of observations, and to compare the observed values of ${\mathrm{E}}_{\mathrm{m}}\left(\mathrm{t}\right)$ to those obtained by solving the integral equation for ${\mathrm{t}}_1<\mathrm{t}\le {\mathrm{t}}_2$

$${\mathrm{E}}_{\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{E}}_{\mathrm{m}}\left({\mathrm{t}}_1\right)+{\mathrm{C}}_{\mathrm{m}}^{-1}\sum _{\mathrm{i}}\underset{{\mathrm{t}}_1}{\overset{\mathrm{t}}{\int }}{\mathrm{I}}_{\mathrm{i}}\left(\tau \right)\mathrm{d}\tau$$

which typically has to be accomplished by numerical means, as closed-form solutions are generally not available. A time-honoured method to compare data to model predictions makes use of the sum of the squares of the differences; ³⁶ this sum is a function of the parameter vector $\{{\mathrm{\kappa }}_1,{\mathrm{\kappa }}_2,\dots ,{\mathrm{\kappa }}_{\mathrm{n}}\}\in R^{+,\mathrm{n}}$ and the value of the latter is sought for which the error sum attains a global minimum.

However, there is no such unique value; instead, in the generic case, there is a subspace of $R^{+,\mathrm{n}}$ , each element of which is such that the extremum of the minimand obtains. This indeterminacy is a general difficulty that arises for both ‘macro’ and ‘micro’ models, even if in the former case it is sometimes dealt with by pretending it is not there.

A pragmatic if somewhat brute-force approach to this identifiability problem is to take more measurements, that is, observe the system under various conditions as it exhibits a wider range of physiological behaviours, and add the sum of squared errors together into one ‘grand sum’ as the minimand. This will never increase the size of the indeterminate minimiser subspace, but it may well reduce it.

There are now two cases to consider: (1) by thus concatenating a sufficient number of experiments, this dimension will be reduced to zero, that is, a unique value for $\{{\mathrm{\kappa }}_1,{\mathrm{\kappa }}_2,\&,{\mathrm{\kappa }}_{\mathrm{n}}\}$ can be found eventually, even if this may require a huge amount of data, or (2) no matter how many observational data are appended, there is an irreducible core of dimension greater than 0. In the latter case, there is an irreducible functional redundancy in the system.

The indeterminate subspace can be characterised as the kernel of a certain operator, which allows its structure to be explored with the customary tools of linear algebra. ²³ For instance, one can determine subsets of ion channels that are substitutable for other such subsets, in the sense of compensatory changes in surface densities of the latter set when the surface densities of the elements of the former are changed, all the while preserving the behaviour observed (i.e. ${\mathrm{E}}_{\mathrm{m}}\left(\mathrm{t}\right)$ for $\mathrm{t}\in [{\mathrm{t}}_1,{\mathrm{t}}_2]$ ). Life scientists often wax philosophical on the subject of functional redundancy: it should be accorded a central role in theoretical biology and its significance is intuitively clear, even if a precise general definition has always been elusive. Here we have a precise experimental and mathematical method to characterise such functional redundancy.

To resolve the identifiability problem, we may choose a value of $\{{\mathrm{\kappa }}_1,{\mathrm{\kappa }}_2,\&,{\mathrm{\kappa }}_{\mathrm{n}}\}$ within the kernel that satisfies an additional constraint, such as minimising the cost for the cell in expressing the ion channels. Using this approach, a putative density vector can be identified and used to predict the behaviour of the electrophysiology of the whole cell assuming a given pharmacodynamics for a specific ion channel. ³⁹ The redundancy analysis plays a key role here: for instance, suppose that we have a candidate drug that knocks out the potassium channel Kir7.1, and we find that this makes the cell more excitable: a putative tropotopic agent. ³⁹ However, if the redundancy analysis shows that the cell can compensate for the loss of Kir7.1 by a suitable up-expression of other channels, this should give us pause; it is not a given that a USMC will do this just because the linear algebra indicates that it can. For instance, on the one hand we observe that a downshift in the BK channel is readily compensated and, on the other hand, knock-out of the hERG channel calls for substantial adjustments to the conductome, with many channels requiring expression levels beyond physiological limits. These predictions are consistent with experimental evidence; inhibition of BK has no effect on uterine contractility, ⁴⁰ but the inhibition of hERG depolarises the membrane and prevents the generation of APs. ⁴¹

Outlook: a computational systems biology platform

Detailed biophysical models of the uterine smooth muscle provide a platform to simulate, accurately at molecular, histological and whole-organ levels, smooth muscle behaviour (e.g. to assess putatively pharmacologically active compounds): this comprises specific ion channels at the molecular level and effects on the propagation of the contraction wave at the whole-organ level, which can be studied in detail on the basis of comprehensive high-resolution digitisation of muscle fibre microarchitecture. To realise such simulations presents numerous challenges, in terms of data acquisition, mathematical foundations and computational integration of data of various modalities, including transcriptomics, histology and electrophysiology.