Uterine smooth muscle cells remain quiescent throughout most of gestation, only generating spontaneous action potentials immediately prior to, and during, labor. This study presents a method that combines transcriptomics with biophysical recordings to characterise the conductance repertoire of these cells, the ‘conductance repertoire’ being the total complement of ion channels and transporters expressed by an electrically active cell. Transcriptomic analysis provides a set of potential electrogenic entities, of which the conductance repertoire is a subset. Each entity within the conductance repertoire was modeled independently and its gating parameter values were fixed using the available biophysical data. The only remaining free parameters were the surface densities for each entity. We characterise the space of combinations of surface densities (density vectors) consistent with experimentally observed membrane potential and calcium waveforms. This yields insights on the functional redundancy of the system as well as its behavioral versatility. Our approach couples high-throughput transcriptomic data with physiological behaviors in health and disease, and provides a formal method to link genotype to phenotype in excitable systems. We accurately predict current densities and chart functional redundancy. For example, we find that to evoke the observed voltage waveform, the BK channel is functionally redundant whereas hERG is essential. Furthermore, our analysis suggests that activation of calcium-activated chloride conductances by intracellular calcium release is the key factor underlying spontaneous depolarisations.
The human genome contains a large number of distinct ion channels, each of which may be assembled into multimeric complexes that modulate the electrical activity of the cell in a subtly different way [1]. Determination of the total ion channel repertoire of a given cell (its ‘conductance repertoire’) is cumbersome via conventional biophysical techniques. The latter rely on the specificity and availability of suitable pharmalogical agents or protocols, which may or may not be able to differentiate between combinations of conductances that can give rise to similar behaviors at the electrophysiological level. By contrast, transcriptomic analysis accurately surveys the complete complement of mRNA coding for all potential conductances. However, such data sets are semi-quantitative at best, absent a straightforward relationship between mRNA levels and surface expression/functionality of the electrogenic proteins. The classic approach in electrophysiology has been to fit a suitable mathematical model to the data, where the parameters to be estimated represent not only the densities of the electrogenic entities but also their biophysical properties [2–4]. This approach is hampered by the large number of distinct electrogenic species. One option is to combine the contributions from several electrogenic species (which are often oligomeric complexes) into single entities, resulting in a ‘macroscopic current’ model [4–10]. Such models may not be sufficiently detailed for pharmacological purposes and accurate assessment of currents within native cells can be technically challenging. The gold standard, therefore, is to determine the conductance repertoire at the level of individual molecular species. The present paper integrates transcriptomics with the information encoded by the action potential (AP) waveform. Electrophysiological data do not necessarily fix a unique solution; there are infinitely many alternative conductance repertoires, that are all equally compatible with the available electrophysiological data. We should consider all conductance repertoires that are consistent with the data, and focus on the properties that they share as a class. This motivates a formal characterisation of the functional redundancy of the system. Our approach is based on two key ideas. First, we use two kinds of data to constrain the space of possibilities, namely mRNA sequencing to determine which species the cell is capable of expressing, combined with parameter estimation to determine the conductance repertoire from electrophysiological data. Second, we obtain the parameters associated with the gating kinetics from the analysis of published data, which means that the only remaining free parameters are the surface densities of the species on the list compiled from transcriptomics. We focus on a cell type of considerable clinical interest, namely the myometrial smooth muscle cell (MSMC), which is the principal unit of electrical activity in the uterus, normally quiescent except in labor, when spontaneous action potentials are generated by the MSMC [11]. A detailed, comprehensive characterisation of these cells would aid the development of therapeutics to manage both preterm birth and perinatal complications associated with uterine contractility, such as postpartum haemorrhage [12]. From a pharmacological point of view, the most promising targets are electrogenic entities whose contribution to the overall electrophysiological behavior is essential in the sense that the role of such an entity cannot be fulfilled by some combination of other ion channels. We show that our analysis is capable of identifying such entities and verify these predictions by means of additional experiments.
Physiological predictions were made, using the simulation model in its free-running mode, on the basis of the conductance repertoire as per the minimal ℓ1-norm criterion. The predictions were validated against experimental data obtained from voltage-clamp experiments for currents isolated by antagonist application. Sets of currents were measured for control and antagonist-inhibited current, and specific current was obtained by subtracting the current recorded in the presence of antagonist from the total current; measurements were reported as mean current at 450 ms ± SEM. The simulations were entirely independent of these additional data. To inhibit and thereby isolate the Kv2. 1 currents, we applied stromatoxin (ScTx), an inhibitor of Kv2. 1 homomeric channels, Kv2. 1/Kv6. 1, Kv2. 1/Kv9. 3, as well as the A-type channel Kv4. 2. The late outward current was measured to separate Kv2. 1 from A-type currents. The measured ScTx-sensitive late outward current (at t = 450 ms), resulting from voltage-clamp experiment at 40 mV after application of 100 nM ScTx, was 3. 71±0. 41 pA/pF (n = 19). The total peak Kv2. 1 current elicited by the simulation for the same voltage step (40 mV) was 3. 77 pA/pF (the channels Kv2. 1, Kv6. 1, and Kv9. 3 contribute to this current). Dofetilide is a selective inhibitor of hERG channels. The isolated dofetilide-sensitive current at 450 ms observed in voltage-clamp experiments at 40 mV after the application of 1 μM dofetilide was 1. 24±0. 61 pA/pF (n = 7), as compared to 1. 544 pA/pF for the simulated current elicited by hERG at the same voltage. Fig 3 shows the experimental current densities under the voltage-clamp conditions. The simulated values of peak currents at selected voltage steps agree, to within the margins of experimental error, with the currents measured for Kv2. 1 channels and hERG channel. The functional redundancy is encoded by the kernel of a linear transformation. An informative graphical representation of this functional redundancy can be obtained by transforming the basis of the kernel to a suitable echelon form, which is then visualised as a heat map. These redundancy maps, shown in Figs 4 and 5, are to be read one row at a time. A unit change in the density of the channel appearing in the leading (left-most) position of each row can be functionally compensated, by upward or downward shifts represented by the trailing elements (i. e. , those to the right). The redundancy map thus represents a compensation signature, charting the cell’s ability of exhibiting the same electrophysiological behavior despite altered channel densities. This compensation is computed on the basis of the voltage waveform that was used to derive the basis vectors, but it is not unreasonable to expect that the compensatory adjustments prescribed by the redundancy map are applicable for a wider range of behaviors. An additional caveat is that these adjustments are exact only for infinitesimal changes in the leading channel density; for finite perturbations in the leading variable, the corresponding shifts dictated by the trailing elements are only a first-order approximation. This means that there is a finite operating range within which the approximation can be expected to be reasonably accurate. Moreover, the operating range is further constrained by considerations of physiologic feasibility: channel densities cannot be negative and, similarly, there will be a upper bound to the current density that can be physically sustained by the cell membrane (there will also be a physiological upper bound to the quantities of channel proteins that can be synthesised and inserted into the membrane). To estimate such upper bounds we used the procedure detailed in Materials and Methods, “Physiological upper bounds on channel density values. ” Uterine contraction during labor is characterised by its automaticity [43]. The myogenic mechanisms underlying uterine activity are not fully understood, although it is generally accepted that a major contribution is made by the increase in the concentration of intracellular calcium of MSMCs via L-type calcium channels [43]. Here we explore three mechanisms that could result in the initiation of periodic spontaneous contractions. The first mechanism acts through the activation of the purinergic membrane receptor (P2RX4), which is an adenosine triphosphate (ATP) -gated ion channel. The second mechanism depends on the periodic depletion of phosphatidylinositol-4,5 bisphosphate (PIP2), which modulates the activation of a number of plasma membrane potassium ion channels. Closure of these channels would promote depolarisation and subsequent contractions. The third mechanism is based on periodic intracellular Ca2+ release from sarcoplasmic reticular stores, leading to activation of the Ca2+-activated Cl− conductance ANO1. Using free-running mode simulations in which the cell was activated through ATP variation (as outlined above), we evaluated the effect of increasing or decreasing various channel densities. The results are shown in Fig 13. When the BK channel density was set to zero, no significant effect was observed on the waveform, in accordance with our earlier characterisation of BK as functionally redundant. On the other hand, variation of Kv2. 1 density, from around 30 channels/pF (double the estimated channel density of Kv2. 1) to 0 channels/pF (complete block of Kv2. 1), had a substantial effect on the waveform (Fig 13). When Kv2. 1 was completely blocked, there was a more marked overshoot of the AP that rapidly depolarised to a plateau level with increased AP frequency. The slight change in the resting membrane potential (RMP) might imply a minor role for Kv2. 1 in maintaining the RMP. In order to determine what role Kv2. 1 currents might play in the whole tissue and to test the predictions of our model, we examined ∼3 × 10 mm muscle strips from the longitudinal layer of murine myometrium under current-clamp conditions using a sharp microelectrode to record fluctuations in membrane potential over time. Typical examples of the activity observed in the myometrial muscle strips are shown in Fig 14A and 14B, where we observe a spontaneous AP generated from a RMPs of −45 mV to −50 mV under control conditions for both Day 15 (D15) and Day-18 (D18), respectively. The application of 100 nM ScTx elicits an increase in AP spike amplitude in the D15 sample and also in the D18 sample (albeit to a lesser extent). The summary data (n = 4 and n = 6 for D15 and D18, respectively, p ≤ 0. 05) for these experiments are presented in the graphs in Fig 14C, 14D, 14E and 14F, where mean AP spike amplitude was significantly increased from 42. 87±5. 07 mV to 53. 25±3. 78 mV in D15 and from 45. 5±4. 96 mV to 56. 21±3. 37 mV (p ≤ 0. 05) in D18 tissue. Application of 100 nM ScTx also significantly increased AP frequency in D15 samples from a mean of 1. 69±0. 18 AP min−1 to 2. 31±0. 34 AP min−1 and 1. 11±0. 22 AP min−1 to 1. 38±0. 19 AP min−1. A modest yet consistent reduction in the duration of AP was apparent in APs recorded from D15 tissue. The number of experiments carried out in the current study was insufficient to detect a significant change in RMP with ScTx application in either D15 or D18 cells (Fig 14F). To assess how Kv2. 1 regulates spontaneous contractions, we recorded isometric tension from ∼3 × 10 mm muscle strips of the longitudinal layer of the mouse myometrium, examining the effects of 100 nM ScTx on the activity. Fig 15A and 15B show representative traces of the mechanical activity of D15 and D18 myometrium. In the presence of 100 nM ScTx, contraction amplitude was increased (13. 32±1. 79 to 14. 24±1. 87 mN, D15; 9. 69±1. 79 to 11. 02±1. 87 mN, D18; Fig 15C), as was frequency (0. 74±0. 05 to 0. 89±0. 07 contractions per minute, D15; 0. 82±0. 06 to 0. 99±0. 07 contractions per minute, D18; Fig 15E), whereas contraction half-width was reduced (22. 97±1. 38 to 17. 65±1. 86 s, D15; not significant in D18; Fig 15D).
Laser-capture microdissection [67] and mRNA-sequencing have led to significant advances in gene expression analysis in both single cells and tissues. Application of these techniques provided a complete repertoire of the mRNA population in MSMCs, which we used to construct a list of possible electrogenic entities by combining all possible combinations of ion channel subunits that are potentially expressed in the MSMC. Since many channels exist as various combinations of subunits, we obtain a repertoire of over thirty electrogenic entities that could be present in the plasma membrane of the MSMC. A mathematical model for every individual potential oligomeric channel complex was formulated on the basis of biophysical data in the literature (typically obtained via heterologous expression). We used native cell behavior to drive the individual conductance currents (viz. voltage and calcium time series as detailed in Materials and Methods). The remaining unknowns were the densities in the plasma membrane of the various species. The observed membrane potential time series is linear in these densities, allowing us to calculate the space of all possible solutions of ion channel densities that are consistent with the observed data; the functional redundancy of the system is represented by the basis vectors of this space. To fix the channel densities we considered the most parsimonious solution defined as having the smallest ℓ1-norm subject to non-negativity. Alternative criteria are possible, such as optimising for agreement with mRNA expression levels for the conductances. However, this approach is conditional on accurate quantitative measurements of expression levels, and on the availability of a reliable mapping from transcriptomics to functional proteomics; neither condition is satisfactorily met at the present state of the art (life times and translation efficiencies of mRNA species vary; moreover, the protein may be stored or degraded instead of being translocated to the membrane). Predictions were made using our model based on the conductance repertoire imposed by the criterion of the smallest ℓ1-norm. These predictions were tested against experimental data, obtained in voltage-clamp experiments for isolated currents that were sensitive to inhibition by antagonist application. The simulated values of peak currents in response to voltage steps closely matched the values that were observed experimentally. In the case of certain channels, for example Kir7. 1, the simulated data failed to predict the observed data accurately, which can be attributed to the fact the particular voltage waveform that was used to drive the model does not distinguish well between currents (e. g. between the Kir7. 1 current and the background potassium current). The myometrium has been modeled at the levels of both tissue and the whole organ [5–7]. Bursztyn et al developed an excitation-contraction model of the uterine muscle cells, comprising voltage-gated calcium channels, calcium transporters, and Na+/Ca2+ exchangers [8]. This model explicitly accounts for the processes of myosin light chain phosphorylation and stress production as per the cross-bridge model of Hai and Murphy [9], but did not include the inactivation process of the calcium current and several other ionic currents. By contrast, Rihana et al [10] modeled the macroscopic currents for each of the relevant ionic species based on voltage-clamp experiments; the dimension of the state space of this model was an order of magnitude lower than ours. Nevertheless, they were able to reproduce various phenomena, such as a single AP, or a train of APs with an RMP of −35 mV to −40 mV and a duration of about 133–200 ms. Another macroscopic ionic current model, accounting for fourteen types of current, was presented by Tong et al [4]. Gating kinetics in the present study were largely derived from equations and parameter estimates taken from the literature and as such representative of the state of the art in the field; accordingly, not all models may have been evaluated to the same standard of rigour as advocated by Fink and Noble [68]. Huys et al [69] proposed an approach similar to ours. However, they illustrated their method on simulated data generated by a macroscopic current model, as opposed to the present work in which real-life data have been confronted with a microscopically detailed account of all relevant currents. One distinguishing feature of our approach is that we can chart the functional redundancy of the system at the level of individual species of electrogenic entities. The redundancy maps show which combinations of channels can be substituted for one another. We investigated the potential scope for redundancy of the system of some of the most-studied potassium channels in MSMCs (e. g. BK, Kv2. 1, SK3, hERG). The use of singular-value decomposition (SVD) to determine the directions of substantial variation in parameter space is well-established [70]. A key advantage of SVD over a naive one-factor-at-a-time approach [71], is that the spectrum of singular values expresses directly which linear combinations of parameters strongly affect the read-out of interest. Sher et al. [72] applied the technique to ionic channel modelling and used the SVD of the sensitivity matrix collecting the partial derivatives of the read-out quantities of interest with respect to the model parameters to probe identifiability and sensitivity, whereas the present method takes a single quantity of interest, the voltage, but treats its entire time course as the output and is accordingly time-global. A regression-based method to perform parameter sensitivity analysis in electrophysiological models was proposed by Sobie [73], who derived quasi-equivalent linear models by treating simulation outputs obtained under randomised parameter variation as statistical input data. Saltelli and Annoni [71] emphasize the essential locality of “one-at-a-time exploration” given the lack of information that would be provided by the higher-order derivatives (including cross terms), but this caveat is not applicable here since the model is linear in the vector of parameters of interest (our model is non-linear in the parameters that were fixed from extraneous data but this is not germane to the point made by these latter authors). The role of BK channels in modulating uterine excitability remains controversial: whilst the presence of the protein and an active conductance are not disputed [74–76], the impact of the BK current is not fully understood [40]. Our functional redundancy analysis indicates that the effect of BK can be compensated by small changes in the channel densities of the other channels in the model. This suggests that the effect which BK exerts on the voltage waveform (as observed in the present study) is readily replaceable by the currents carried by other electrogenic entities. This observation is consistent with the findings of Aaronson et al, who demonstrated that application of BK inhibitors has no effect on spontaneously contracting rodent myometrial strips [40]. If BK can be readily substituted by suitable combinations of other channels, one might ask why it should be part of the conductance repertoire at all. It may be relevant that the unitary conductance of the BK channel is unusually large [77], which means that the cell can be silenced by expressing even a minor number of these molecules. Perhaps BK serves the function of an “emergency shut-down. ” We also investigated the potential involvement of Kv2. 1 in shaping the action potential (AP). We investigated the effect of changes in Kv2. 1 channel density on the AP waveform by running free-running simulations. Decreasing the Kv2. 1 channel density from the estimated value to zero increased the amplitude of the initial AP as well as the rate at which steady state is reached during the plateau phase of the AP. The predicted effects on AP amplitude and frequency were confirmed by experiments although the effect on the resting membrane potential (RMP) could not be verified with certitude. It appears that Kv2. 1 channels do play a role in modulating the excitability of the myometrium by inhibiting both contraction frequency and amplitude, while contributing little to the RMP. As regards the hERG channel, our functional redundancy analysis indicated that small changes in the hERG channel density could not be compensated by physiologically realistic adjustments in the densities of the other channels. In addition, the simulations carried out free-running mode suggest that variations in the hERG channel density substantially affect the AP waveform (Fig 13). Taken together, these findings indicate that hERG is a major suppressor of bursting activity in MSMCs, in keeping with Greenwood et al [41], who observed that the ERG-specific blocker dofetilide-induced contractions in quiescent tissues, while ERG-channel blockers were able to induce contractions in tissue strips that failed to develop spontaneous rhythmic activity. In contrast to hERG’s essential role, the SK channels were found to be readily susceptible to compensation. The redundancy analysis indicates that SK2 can be functionally replaced by adjustments in SK3 and SK4 densities. This is in accordance with the close physiological similarities between the members of the SK family. We also demonstrated the utility of the model in exploring various mechanistic hypotheses, addressing three mechanisms that could trigger spontaneous contractions. Simulations provided supporting evidence for the involvement of extracellular ATP and/or calcium ‘sparks’ affecting CaCCs, whereas we found no support for a role of PIP2 acting on potassium channels. Although far from conclusive in their own right, such simulations can help experimentalists to pinpoint promising research avenues. The model can be improved in several ways. In its present form, the dynamics of intracellular calcium are represented in a minimalistic way and a more physiologically faithful model would be desirable. However, this is subject to the quality of the temporospatial resolution of the cytosolic and sarcoplastic-reticular levels of Ca2+. Another mechanism we have not considered here is the feedback from force-generation on the dynamics of membrane potential, which acts via stretch-sensitive ion channels. Furthermore, our model describes a single cell, whereas the in vivo activity of MSMCs is critically dependent on their interconnectivity. Thus, a multi-cell model with realistic network properties would constitute another future improvement. In the present study a spontaneous voltage signal was imposed. It is also possible to impose artificial voltage waveforms (AVWs), for instance composed as a linear combination of B-splines, and compare the current response to such an AVW to the current response of the model. The advantage of this approach would be that the AVW can be optimised to discriminate maximally between two given conductances, for instance by minimising the inner product between the current responses of the respective conductances, where the two are considered as functions of time. A tighter upper bound on the functional redundancy of the system can be obtained in this manner. The functional redundancy we map is strictly relative to the waveform that is imposed. As more behaviors are explored, more “essential roles” are revealed and the dimensionality of the kernel can generically be expected to diminish; in fact, it could come down all the way to zero, which would imply that there is no intrinsic redundancy in the conductance repertoire. Alternatively, the kernel dimensionality settles on a non-zero value “in the limit of arbitrarily many explored behaviors”—this value would then denote the irreducible functional redundancy of the system. It stands to reason that sundry physiological limitations will impose further constraints on the system, but the considerations in the present paper concern only the functional redundancy vis-à-vis the membrane conductance repertoire. Two key advantages of our approach are that it is modular and expandable. The model can readily accommodate new information regarding additional conductances or novel kinetic properties. Moreover, being expressed in terms of microscopic conductances that correspond to particular molecular entities, the model can in principle be related directly to transcriptomics. This has important practical advantages, because transcriptomics data can be acquired in different tissues, and at different time points, accurately surveying all potential molecular entities. Another facet of the expandability of our approach is that the characterisation of functional redundancy, becomes more accurate, that is, comes closer to the irreducible functional redundancy, as a wider range of physiological behavior of a given cell (both spontaneous and evoked) is observed. Our analysis of functional redundancy provides a means to address the well-known difficulty in electrophysiological modelling that electrophysiology at the whole-cell level is insufficient to constrain the densities of all channels. Traditionally, this problem has been solved by postulating macroscopic currents that represent the combined action of several similar conductances, but this approach has the drawback that a direct correspondence to genomic or transcriptomic data is lost. We would even argue that the approach based on macroscopic currents is a false economy since, from an experimental point of view, the acquisition of the macroscopic currents is itself underdetermined, whereas microscopic currents can be unequivocally isolated and thereby fully characterised.
All procedures were conducted within the guidelines of The Declaration of Helsinki and were subject to local ethical approval (REC-05/Q2802/107). Prior to surgery, informed written consent for sample collection was obtained. The complete repertoire of the electrogenic proteins that are potentially expressed in the MSMC was determined on the basis of mRNA expression data, which have been made available as entire raw data via the GEO database (GEO series accession number GSE50599) and, in processed form, as the supplementary material associated with the original publication by Chan et al [13]. We constructed a complete list of every individual oligomeric channel complex that has been attested in the literature and that is consistent with the subunits in the mRNA expression set; the mathematical model, shown in Fig 1A, incorporates 31 time-dependent ionic currents. These currents include outward currents such as the voltage-gated potassium current, the voltage- and calcium-gated potassium current, as well as the calcium-gated potassium current. In addition, the model comprises two inward, depolarising currents attributed to the two voltage-gated calcium channels (L-type and T-type). Chloride fluxes are represented as a calcium-activated chloride current, as well as a background current. A background potassium current was also accounted for. Finally, the model includes pumps and exchangers: the Na+-Ca2+ exchanger (NCX), the plasma membrane calcium ATPase (PMCA), and the Na+-K+ pump (NaK). A separate mathematical model for each oligomeric channel was formulated, either based on models that have already been proposed in the literature or on the basis of the available data. The biophysical and kinetic parameters of several of these entities could in many cases be taken directly in the literature, or else have been obtained by means of least-squares fitting to the experimental data on heterologous expression systems, taken from the literature, as detailed below. The mathematical expressions that describe the individual channel kinetics, together with the associated parameters descriptions, units, and values were previously listed in McCloskey et al [39] and are reproduced here for the sake of clarity in the Supplementary Materials, while a summary of the potential conductance species included in our model is shown in Table 1.
A well-known problem in electrophysiologal modeling is that the parameters of the gating kinetics of the ion channels cannot be uniquely determined from observed behavior at the cellular level. One solution is to employ simplified “macroscopic” currents that mimic the behavior of aggregates of distinct entities at the protein level. The gating parameters of each channel or pump can be determined by studying it in isolation, leaving the general problem of finding the densities at which the channels occur in the plasma membrane. We propose an approach, which we apply to uterine smooth muscle cells, whereby we constrain the list of possible entities by means of transcriptomics and chart the indeterminacy of the problem in terms of the kernel of the corresponding linear transformation. A graphical representation of this kernel visualises the functional redundancy of the system. We show that the role of certain conductances can be fulfilled, or compensated for, by suitable combinations of other conductances; this is not always the case, and such “non-substitutable” conductances can be regarded as functionally non-redundant. Electrogenic entities belonging to the latter category are suitable putative clinical targets.
voltage-gated ion channels
uterus
calcium-activated potassium channels
medicine and health sciences
reproductive system
membrane potential
electrophysiology
neuroscience
ion channels
cellular structures and organelles
proteins
calcium channels
biophysics
cell membranes
physics
biochemistry
myometrium
anatomy
cell biology
physiology
biology and life sciences
potassium channels
physical sciences
voltage-gated calcium channels
neurophysiology
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