Cellular Interrogation: Exploiting Cell-to-Cell Variability to Discriminate Regulatory Mechanisms in Oscillatory Signalling

The molecular complexity within a cell may be seen as an evolutionary response to the external complexity of the cell’s environment. This suggests that the external environment may be harnessed to interrogate the cell’s internal molecular architecture. Cells, however, are not only nonlinear and non-stationary, but also exhibit heterogeneous responses within a clonal, isogenic population. In effect, each cell undertakes its own experiment. Here, we develop a method of cellular interrogation using programmable microfluidic devices which exploits the additional information present in cell-to-cell variation, without requiring model parameters to be fitted to data. We focussed on Ca2+ signalling in response to hormone stimulation, which exhibits oscillatory spiking in many cell types and chose eight models of Ca2+ signalling networks which exhibit similar behaviour in simulation. We developed a nonlinear frequency analysis for non-stationary responses, which could classify models into groups under parameter variation, but found that this question alone was unable to distinguish critical feedback loops. We further developed a nonlinear amplitude analysis and found that the combination of both questions ruled out six of the models as inconsistent with the experimentally-observed dynamics and heterogeneity. The two models that survived the double interrogation were mathematically different but schematically identical and yielded the same unexpected predictions that we confirmed experimentally. Further analysis showed that subtle mathematical details can markedly influence non-stationary responses under parameter variation, emphasising the difficulty of finding a “correct” model. By developing questions for the pathway being studied, and designing more versatile microfluidics, cellular interrogation holds promise as a systematic strategy that can complement direct intervention by genetics or pharmacology.

ImageJ macro used for segmentation of cells after imaging. 28 Image segmentation and Ca 2+ trajectory generation. 29 Heterogeneity control for CIAD imaging. 29 Pseudo-code of the spike filtering and pattern search algorithms. 30 Flow chart and pseudo-code for choosing random ICs and PVs and determining the natural period. 31 Flow chart and pseudo-code for pattern counting. 32 Absence of photobleaching. 33 Flow chart and pseudo-code for nonlinear amplitude analysis. 34 Effect of sample size on nonlinear frequency analysis under uniform sampling 35 Effect of sample size on nonlinear frequency analysis under lognormal sampling 36 Effect of sample size on nonlinear amplitude analysis under uniform sampling 37 Effect of sample size on nonlinear amplitude analysis under lognormal sampling 38

General comments about the models
Eight models of oscillatory Ca 2+ spiking were selected from the literature to cover a range of behaviours, as described in the Paper. The dynamical variables used in the various models are listed in Table 1 below. The differential equations for the models are listed in the sub-sections that follow, using the abbreviations species variable units cytoplasmic calcium cc µM ER calcium cer µM IP 3 p µM fraction of active IP 3 receptors n u.l. fraction of receptors in state x x u.l. fraction of receptors in state y y u.l. fraction of receptors in state y 2 y 2 = x + y u.l.  Figure 3. These equations follow the model descriptions in the original papers. Also listed in the accompanying table are the reference parameter values used in the original paper, at which the model exhibits limit cycle oscillations in response to step stimulation.
In the original papers, an input stimulation by hormone is regarded as varying one of the parameter values, so that the corresponding term in the equations is effectively taken to aggregate the effects of the hormone receptor and any downstream signalling between the receptor and the mechanism in which the parameter appears. We did not want to modify the original models and thereby make it difficult to compare our work to previous studies of calcium oscillation. We therefore chose to use the same approach. The input parameter for each model is flagged with an asterisk, * , in the corresponding table of reference parameter values.

Meyer-Stryer-MST
This is from the 1988 paper by Meyer and Stryer, "Molecular model for receptor-stimulated calcium spiking" [4].

Goldbeter-Dupont-Berridge-GDB
This is from the 1990 paper by Goldbeter, Dupont and Berridge, "Minimal model for signal-induced Ca2+ oscillations and for their frequency encoding through protein phosphorylation" [2].
variable init. cond. cer 0.25 µM cc 1.5 µM Table 3: GDB reference parameter values (left) and initial conditions (right) from [2] This is from the 1993 paper by Atri, Amundson, Clapham and Sneyd, "A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte" [1], with modifications for the class I/II distinction from the 2006 paper by Sneyd, Tsaneva-Atanasova, Reznikov, Bai, Sanderson and Yule, "A method for determining the dependence of calcium oscillations on inositol trisphosphate oscillations" [6].

Atri et al class 2-AT2
This is also from the 1993 paper by Atri, Amundson, Clapham and Sneyd, "A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte" [1], with modifications for the class I/II distinction from the 2006 paper by Sneyd, Tsaneva-Atanasova, Reznikov, Bai, Sanderson and Yule, "A method for determining the dependence of calcium oscillations on inositol trisphosphate oscillations" [6].
3.915194 µM cc 0.083406 µM µ 0 and b to J IP 3 Table 5: AT2 reference parameter values (left) and initial conditions (right) from [6] This is from the 1994 paper by Li and Rinzel, "Equations for InsP 3 receptor-mediated [Ca 2+ ] i oscillations derived from a detailed kinetic model: a Hodgkin-Huxley like formalism" [3], with modifications for the class I/II distinction from the 2006 paper by Sneyd, Tsaneva-Atanasova, Reznikov, Bai, Sanderson and Yule, "A method for determining the dependence of calcium oscillations on inositol trisphosphate oscillations" [6].

Li-Rinzel class 2-LR2
This is also from the 1994 paper by Li and Rinzel, "Equations for InsP 3 receptor-mediated [Ca 2+ ] i oscillations derived from a detailed kinetic model: a Hodgkin-Huxley like formalism" [3], with modifications for the class I/II distinction from the 2006 paper by Sneyd, Tsaneva-Atanasova, Reznikov, Bai, Sanderson and Yule, "A method for determining the dependence of calcium oscillations on inositol trisphosphate oscillations" [6]. 1.8 Sneyd-LeBeau two and three state-SB2 and SB3 These are from the 2000 paper by Sneyd, LeBeau and Yule, "Traveling waves of calcium in pancreatic acinar cells: model construction and bifurcation analysis" [5].
2 state: variable init. cond.  1.9 Hybrid AT1-based and LR1-based models The 11 hybrid models are listed below, following the numbering scheme in the paper, with the AT1-based hybrid given first and the corrresponding LR1-based hybrid given second. The specific modification is explained and the modified term is marked with an asterisk.

Hybrid 1
-AT1: replaced the hyperbolic function in J cc by a Hill function of coefficient 2.
replaced the Hill function in J cc by a hyperbolic one.

Hybrid 2
Replaced hyperbolic dependence of J IP 3 on cc by a cubic function as in LR1.
replace the cubic dependence of J IP 3 on cc by a linear one.

Hybrid 3
Did the same thing as in Hybrid 2 but for IP 3 , using the term µ(p) 3 replace the cubic dependence of J IP 3 on p by a linear one.

Hybrid 4
Replaced the linear dependence of J IP 3 on n by a cubic one.
replace the cubic dependence of J IP 3 on n by a linear one.
version: combine the changes in Hybrids 2, 3 and 4.

Hybrid 6
-AT1 version: starting with Hybrid 5, set b = 0 and µ 0 = 0 Hybrid 7 -AT1 version: assumed a hyperbolic function for n inf (cc) Hybrid 10 -AT1 version: set b = 0 and µ 0 = 0  The row in the middle, coloured in red, corresponds to the sample of 20,000 sets used in Paper Figure 5. The four rows above the middle were obtained by random sub-sampling of 5000 sets from among those 20,000 sets. The four rows below the middle were obtained by independently generating further samples of 20,000 sets.