Communication shapes sensory response in multicellular networks
Abstract
Collective sensing by interacting cells is observed in a variety of biological systems, and yet a quantitative understanding of how sensory information is collectively encoded is lacking. Here we investigate the ATPinduced calcium dynamics of monolayers of fibroblast cells that communicate via gap junctions. Combining experiments and stochastic modeling, we find that increasing the ATP stimulus increases the propensity for calcium oscillations despite large celltocell variability. The model further predicts that the oscillation propensity increases not only with the stimulus, but also with the cell density due to increased communication. Experiments confirm this prediction, showing that cell density modulates the collective sensory response. We further implicate cellcell communication by coculturing the fibroblasts with cancer cells, which we show act as “defects” in the communication network, thereby reducing the oscillation propensity. These results suggest that multicellular networks sit at a point in parameter space where cellcell communication has a significant effect on the sensory response, allowing cells to simultaneously respond to a sensory input and to the presence of neighbors.
I Significance Statement
Cells routinely sense and respond to their environment, and they also communicate with each other. Communication is therefore thought to play a crucial role in sensing, but exactly how this occurs remains poorly understood. We study a population of fibroblast cells that responds to a chemical stimulus (ATP) and communicates by molecule exchange. Combining experiments and mathematical modeling, we find that cells exhibit calcium oscillations not only in response to the ATP stimulus, but also in response to an increase in cell density. To confirm that the density dependence is a result of increased cellcell communication, we combine the fibroblasts with cancer cells, which we show have weakened communication properties. The oscillations indeed decrease with the fraction of cancer cells. Our results show that when cells are together, their sensory responses reflect not just the stimulus level, but also the degree of communication within the population.
Ii Introduction
Decoding the cellular response to environmental perturbations, such as chemosensing, photosensing, and mechanosensing has been of central importance in our understanding of living systems. Up to date, most studies of cellular sensation and response have focused on single isolated cells or population averages. An emerging picture from these studies is the set of design principles governing cellular signaling pathways: these pathways are organized into an intertwined, often redundant network, and the network architecture is closely related with the robustness of cellular information processing Lauffenburger2000 (); Darnell2000 (). On the other hand, many examples suggest that collective sensing by many interacting cells may provide another dimension for the cells to process environmental cues Barritt1994 (). Quorum sensing in bacterial colonies Bassler2001 () and cyclic adenosine monophosphate (cAMP) signaling in Dictyostelium discoideum populations Sawai2010 () generate dramatic group behaviors in response to environmental stimuli. Various types of collective sensing have also been observed for higher level multicellular organisms, such as in the olfaction of insects Vosshall2007 () and mammals Smear2011 (), the glucose response in the pancreatic islet benninger2008gap (), and the visual processing of animal retinal ganglion cells Bialek2006 (); Chichilnisky2011 (). These examples suggest a fundamental need to revisit cellular information processing in the context of multicellular sensation and responses, since even weak celltocell interaction may have strong impacts on the states of multicellular network dynamics Bialek2006 (). In particular, we seek to examine how the sensory response of cells in a population differs from that of isolated cells, and whether we can tune between these two extremes by controlling the degree of cellcell communication.
Previously, we have described the spatialtemporal dynamics of collective chemosensing of a mammalian cell model system Sun2012 (); Sun2013b (). In this model system, high density mouse fibroblast cells (NIH 3T3) form a monolayer that allows nearestneighbor communications through gap junctions Gilula1996 (). When extracellular ATP is delivered to the cell monolayer, P2 receptors on the cell membrane trigger release of the second messenger IP3 (Inositol 1,4,5trisphosphate) Irvine1989 (). IP3 molecules, when binding to their receptors at the membrane of endoplasmic reticulum (ER), open ion channels to allow influx of Ca from the ER to the cytoplasm. The fast increase in cytosolic calcium concentration is followed by a slow decrease when Ca are pumped back to the ER Snyderc1990 (). This simple picture, however, is complicated by the nonlinear feedback between Ca and ion channel opening probability, which leads to rich dynamic behaviors such as cytosolic calcium oscillations Putney2011 (). In the situation of collective ATP sensing, we have found that gap junction communications dominate intercellular interactions Sun2012 (). Furthermore, these shortrange interactions propagate and turn the cell monolayer into a percolating network Sun2013b (). These characteristics make the system ideal for studying how sensory responses are modulated by extensive communication in multicellular networks.
Here we use this model system to examine how cellcell communication affects collective chemosensing. Combining experiments with stochastic modeling, we find that cells robustly encode the ATP stimulus strength in terms of their propensity for calcium oscillations, despite significant celltocell variability. The model further predicts that the oscillation propensity depends not only on the stimulus, but also on the density of cells, and that denser monolayers have narrower distributions of oscillation frequencies. We confirm both predictions experimentally. To verify that the mechanism behind the density dependence is the modulation of cellcell communication, we introduce cancer cells (MDAMB231) into the fibroblast cell monolayer. As we show, MDAMB231 cells act as “defects” in the multicellular network, as they have distinct calcium dynamics when compared with the fibroblasts due to reduced gap junction communication Kanno1966 (); Laird2006 (); Jiang2014 (). We find that the oscillation propensity of the fibroblasts decreases as the fraction of cancer cells increases, confirming that the sensory response is directly affected by the cellcell communication. Our findings indicate that cells’ collective response to a sensory stimulus is distinct from their individual responses, and that in this case the response simultaneously encodes both the stimulus strength and the degree of communication within the population.
Iii Results
In order to study the sensory responses of a multicellular network, we use single channel microfluidic devices and deliver ATP solutions to monolayers of fibroblast (NIH 3T3) cells (see Appendix A). The stimuli solutions of ATP concentrations varying from 0 M to 200 M are pumped though the channel at rate of 50 l/min, ensuring minimal flow perturbation to the cells and fast delivery of ATP across the field of view. During this time we record the calcium dynamics of each individual cell loaded with fluorescent calcium indicator (Fluo4, Life Technologies) at 4 frames/sec (Hamamatsu Flash 2.8, Hamamatsu). Our analysis uses 400 seconds of data from each recording selected such that ATP arrival is at approximately 50 seconds.
We modulate the degree of communication in two ways. First, we vary the cell density. Smaller cell densities correspond to larger celltocell distances, which we expect to reduce the probability of forming gap junctions. Second, we coculture the fibroblasts with breast cancer (MDAMB231) cells in the flow channel (see Appendix A). As we later show, MDAMB231 cells have reduced communication properties and therefore act as defects in the multicellular network. To distinguish the two cell types, MDAMB231 cells are prelabeled with red fluorescent dye (Cell Tracker Red CMTPX, Life Technologies). Varying cell density and the fraction of cancer cells allows us to control the architecture of the multicellular network over a wide range.
Figure 1A shows the composite image of a highdensity cell monolayer with cocultured fibroblast and cancer cells. In this example, MDAMB231 cells make up a fraction 15% of the total population which has a total cell density of 2500 cells/mm. At this density, each cell has an average of 6 nearest neighbors from which extensive gap junction communication is expected. After identifying cell centers from the composite image (see Appendix A), we compute the timedependent average fluorescent intensity near the center of each cell which represents the instantaneous intracellular calcium concentration at the single cell level.
iii.1 Collective response to ATP stimuli
Typical responses of cells to excitation by ATP are shown in Fig. 1B. We see that, on average, higher concentrations of ATP trigger larger increases in calcium levels. Celltocell variations are significant; for example, response times as well as subsequent calcium dynamics of individual cells vary dramatically amongst cells in the same multicellular network. In many cells, the initial calcium increase is followed by transient calcium oscillations. We quantify the oscillation propensity by computing the fraction of nonoscillating cells using a peakfinding algorithm (see Appendix B). We see in Fig. 1C that higher concentrations of ATP cause a larger percentage of cells to oscillate, and thus a smaller .
The period of the oscillation is characterized by the interspike interval (ISI), which has been proposed to dynamically encode information about the stimuli Tang1995 (); Falcke2014 (). To investigate the characteristics of ISI in the context collective chemosensing, we systematically study the statistics of the ISI from 30,000 cells. Figure 1D shows the histogram (event counts) of ISI values, normalized by the number of cells, of a typical experiment where the ATP concentration is 50 M. We see that the distribution is broad, which underscores the high degree of celltocell variability in the responses. Figure 1E summarizes the distribution at each ATP concentration using a boxandwhisker plot. We see that there is no significant dependence of the ISI on the ATP concentration. This is at odds with a familiar property of calcium oscillations, termed frequency encoding, in which the oscillation frequency (or ISI) depends on the strength of the stimulus Putney2011 (); Tang1995 (); Woods1986 (); Meyer1991 (). However, we will see in the next section that the lack of a dependence here is likely due to the high degree of celltocell variability.
Finally, we characterize the spatial correlations of the ISI within the monolayer by computing the crosscorrelation function as a function of topological distance between cells. In particular, for each experiment, we first identify all oscillatory cells and compute the average ISI for each cell . We then define , and , where is the topological distance between cells and . Figure 1F shows that falls off immediately for . This is surprising, since the cells experience identical ATP stimuli, and one might hypothesize that communication between cells would result in the ISI values for nearby cells being correlated. However, as described next, evidence from mathematical modeling suggests that this correlation is removed by the celltocell variability.
iii.2 Stochastic modeling of the collective response
To obtain a mechanistic understanding of the experimental observations, we turn to mathematical modeling. We develop a stochastic model of collective calcium signaling based on the work of Tang and Othmer Tang1995 (); Othmer1993 (). Their model captures the ATPinduced release of IP3, the IP3triggered opening of calcium channels in the ER, and the nonlinear dependence of the opening probability on the calcium concentration, as schematically illustrated in Fig. 2A. The model neglects more complex features of calcium signaling observed in some cell types, such as upstream IP3 oscillations Stryer1988 (); Hofer2006 () and spatial correlations among channels Champeil1999 (); Falcke2004 (). The model predicts that at a critical ATP concentration, the calcium dynamics transition from nonoscillating to oscillating. However, it was previously only analyzed deterministically for a single cell Tang1995 (); Othmer1993 (). Therefore we extend it to include both intrinsic noise and cellcell communication via calcium exchange (see Appendix C). We also explicitly include the dynamics of IP3, which has a constant degradation rate and a production rate that we take as proportional to the ATP concentration. We simulate the dynamics using the Gillespie algorithm Gillespie1977 (), and we vary the density of cells on a square grid, which modulates the degree of communication.
Figure 2B shows the dependence of on , where is the model analog of the experimentally controlled ATP concentration. Consistent with the experimental findings in Fig. 1C, we see that decreases with . In the model, the decrease is due to the fact that intrinsic noise broadens the transition from the nonoscillating to the oscillating regime. Thus, instead of a sharp transition from to as predicted deterministically, the transition occurs gradually over the range of shown in Fig. 2B. Figure 2C shows the dependence of the ISI on in the model (see the green box plots). We see that the ISI decreases with , which is expected since frequency encoding is a component of the TangOthmer model Tang1995 (); Othmer1993 (). Yet, this property is not consistent with the experimental observation in Fig. 1E, where the ISI shows no clear dependence on ATP concentration. Further, Fig. 2D shows the dependence of the correlation function on the topological distance in the model (see the green curve). We see that decreases gradually with , indicating nonzero spatial correlations in the ISI. This feature is again inconsistent with the experimental findings seen in Fig. 1F.
Motivated by the high level of celltocell variability evident in Figs. 1B and 1D, we hypothesize that celltocell variability is responsible for these discrepancies between the model and the experiments. Indeed, inspecting the ISI histogram from the model reveals a very narrow distribution of ISI values, as seen in Fig. 2E (green bars), which is in contrast to the broad distribution observed experimentally in Fig. 1D. To incorporate celltocell variability, we allow the model parameters to vary from cell to cell. Lacking information about the susceptibility of particular parameters to variation, we allow all model parameters to vary by the same maximum fold change (i.e., parameters are sampled uniform randomly in log space). The maximum fold change is found by equating the variance of the resulting ISI distribution with that from the experiments, which yields the value (see Appendix C). As seen in Fig. 2E (blue bars), the resulting ISI distribution is consistent with that observed in Fig. 2D, both in width and in shape.
We see in Fig. 2C (blue box plots) that including celltocell variability in the model severely weakens the decrease of the ISI with , therefore agreeing with the experimental results shown in Fig. 1E (with variability, is defined as the mean of the values sampled for each cell). We also see in Fig. 2D (blue curve) that variability removes the correlations for , which is consistent with the immediate falloff observed experimentally in Fig. 1F. Importantly, even with variability, the decrease of with seen in Fig. 2B persists, as demonstrated in Fig. 2F. This decrease remains consistent with the experimental observation in Fig. 1C. Indeed, variability significantly broadens the range of values over which the transition occurs, as expected (compare Fig. 2F to Fig. 2B), which is consistent with the broad range over which the transition occurs experimentally (Fig. 1C).
iii.3 Effects of communication on the sensory response
Having validated the model, we now use it to make predictions about the effect of cellcell communication on collective calcium dynamics. Communication in the model is controlled by cell density, with higher density leading to more celltocell contacts and thus a higher degree of communication. Therefore we first investigate the dependence of the oscillation propensity on the cell density. Fig. 3A shows as a function of both cell density and the ATPinduced IP3 production rate . We see that the fraction of nonoscillating cells transitions from to as a function of and that there is also a dependence of on . Not only has intrinsic noise and celltocell variability broadened the transition as a function of , but evidently cellcell communication has introduced a dependence of the transition on cell density. The result is that at low , is everywhere large and is independent of cell density (Fig. 3B and left box in Fig. 3A). However, at intermediate , is a decreasing function of cell density (Fig. 3C and right box in Fig. 3A). In this regime, increasing the degree of communication causes more cells to exhibit oscillatory calcium dynamics (thus decreasing ), even with a fixed sensory stimulus . At large (beyond the range of Fig. 3A), we have checked that the nonoscillating fraction is driven to low values as expected, and the density dependence of diminishes.
The prediction in Fig. 3C is striking because it implies that cellcell communication causes more cells to oscillate, even while celltocell variability causes their ISI values to be spatially uncorrelated (Fig. 2D). Therefore, we wondered whether communication would have an effect on the width of the ISI distribution in this regime. The width, or more generally the amount of uncertainty contained in the ISI distribution, is characterized by the entropy. For a continuous variable , the entropy becomes the differential entropy, defined as , where is the probability density. As seen in Fig. 3E, the entropy of the ISI distribution increases with . This indicates that as communication decreases the nonoscillating fraction, it also narrows the distribution of ISI values.
We now test these predictions in our experimental system. To test our predictions about how the nonoscillating fraction of cells should depend on cell density, we measure as a function of for various ATP concentrations. We see in Fig. 4A that with no ATP, is large at both low and high densities, and there is no statistically significant correlation between and . Then, we see in Fig. 4B that at intermediate ATP concentrations ( M), significantly decreases with (see Appendix B for pairwise statistical comparison of particular density values). Finally, we see in Fig. 4C that at large ATP concentration (200 M), is small at both low and high densities, and again there is no statistically significant correlation between and . These results confirm the predictions in Fig. 3.
To test the prediction that the entropy of the ISI distribution increases with the nonoscillating fraction of cells, we measure as function of . As seen in Fig. 4D, increases with , consistent with the prediction in Fig. 3D. This implies that cell density regulates the spectrum of the ISI response. In particular, it suggests that increasing the degree of cellcell communication narrows the distribution of ISI, making the ISI values less variable across the cell population. We have also checked that the entropy of the distribution of crosscorrelation values for nearest neighbors’ entire calcium trajectories Sun2012 (); Sun2013b () decreases as a function of cell density (see Appendix D). is not only a statistical characterization of the collective cellular dynamics, but also may provide robust channel for cellular information encoding Bialek2006 (); Meister1996 (); Meister2003 (). Together, these results imply that cellcell communication has a significant effect on the collective sensory response. This finding is especially striking given the strong effects of celltocell variability (Fig. 1E and F). We conclude that the effects of communication observed here persist in spite of extensive celltocell variability.
iii.4 Effect of cancer cell “defects”
We have seen that increasing cell density increases the propensity of cells to oscillate in response to an ATP stimulus. This behavior is consistent with our model, which predicts that the mechanism is through increased cellcell communication. However, it could be in the experiments that increasing the cell density introduces other effects beyond increased gapjunction communication, such as mechanical coupling between cells or coupling to the substrate Janmey2011 (). To modulate the communication directly, we vary the fraction of cancer cells with which the fibroblasts are cocultured, while keeping the density of all cells fixed. Since cancer cells are known to have reduced gap junction communication Kanno1966 (); Laird2006 (); Jiang2014 (), we expect the fraction of nonoscillating cells to have the opposite dependence on that it does on cell density (Fig. 4B).
First we investigate whether MDAMB231 cells indeed have reduced communication in our system. Figure 5A shows several examples of singlecell calcium dynamics for NIH 3T3 and MDAMB231 cells in a typical experiment. We see that both cell types exhibit immediate increases in cytosolic calcium levels at the arrival of ATP, but cancer cells typically show long relaxation times while fibroblast cells tend more often to exhibit oscillations after stimulation. These qualitative features are maintained across all ATP concentrations. Figure 5B shows a comparison of the intercellular diffusion coefficients in the two cell types, obtained from a fluorescence recovery after photobleaching (FRAP) analysis Didelon2008 () (see Appendix A). We see in Fig. 5B that gap junctionmediated diffusion between MDAMB231 cells is significantly weaker than between NIH 3T3 cells, consistent with previous reports Kanno1966 (); Laird2006 (); Jiang2014 (). It is therefore evident that MDAMB231 cells can be treated as communication defects in the cocultured multicellular network. Indeed, Fig. 5C shows the spatial distribution of these defects in the monolayer. In Fig. 5C, the mean ISI for each cell is shown in color, with nonoscillating cells in black. We see that cancer cells, labeled by white circles, are more likely to be nonoscillating, which is consistent with the qualitative characteristics shown in Fig. 5A. We have further quantified the distinction between the two cells types in Appendix B, where we show using the distributions of ISI values that oscillatory events are at least five times less likely to occur for the MDAMB231 cells.
Having established that the presence of cancer cells reduces the degree of cellcell communication in the monolayer, we now vary the fraction of cancer cells and measure the oscillation propensity of the remaining fibroblasts. Fig. 5D shows the nonoscillating fraction of fibroblasts (blue bars) as a function of the cancer cell fraction for a typical experiment at fixed cell density ( cells/mm). We see that significantly increases with (see Appendix B for pairwise statistical comparison of particular values). We also see that for all cells (both fibroblasts and cancer cells, red bars) significantly increases with , and that, as expected, is larger for all cells than for just fibroblasts. These findings imply that reduced cellcell communication decreases the propensity for calcium oscillations, which is consistent with the effects of varying cell density (Fig. 4B). Finally, we also investigate the effect of cancer cells on the entropy of the ISI distribution. As shown in Appendix B, is higher for cells that are surrounded by a large number of cancer cells, and lower for cells with pure fibroblast neighbors. In the latter case, also increases as the number of nearest neighbors decreases. These findings imply that reduced cellcell communication increases the entropy of the ISI values, even at the local level of a cell’s microenvironment, which is consistent with the effects seen in Fig. 4D. Taken together, we conclude that the calcium dynamics of individual cells are strongly regulated by the degree of gap junction communication inside the cell monolayer.
Iv Discussion
We have characterized the collective calcium dynamics of multicellular networks with varying degrees of cellcell communication when they respond to extracellular ATP. These networks consist of cocultured fibroblast (NIH 3T3) and cancer (MDAMB231) cells. We have revealed their properties using a combination of controlled experiments and stochastic modeling. We have found that increasing the ATP stimulus increases the propensity for cells to exhibit calcium oscillations, which is expected at the single cell level. However, we have also found that increasing the cell density alone, while keeping the stimulus fixed, has a similar effect, revealing a purely collective component to the sensory response. Modeling suggests that this effect is due to an increased degree of molecular communication between cells. In line with this prediction, we have found that increasing the fraction of cancer cells in the monolayer reduces the oscillation propensity, as cancer cells act as defects in the communication network.
Typical physiological concentrations of extracellular ATP are tens of nM and larger schwiebert2003extracellular (); falzoni2013detecting (); trabanelli2012extracellular (), but hundreds of M has been associated with disease pellegatti2008increased (). This suggests that healthy tissues in vivo are restricted to the regimes of Fig. 4A and B, and not C ( M). In the regime of Fig. 4A, the calcium dynamics encode the stimulus strength. In the regime of Fig. 4B, the calcium dynamics encode two pieces of information simultaneously: the stimulus strength and the cell density. One may wonder why a system would evolve to encode two pieces of information in the same quantity. It is tempting to speculate that there is a benefit to having the sensory response of the system be a function of the population density, as in quorum sensing Bassler2001 () or socalled dynamical quorum sensing based on collective oscillations mehta2010approaching (); de2007dynamical (); taylor2009dynamical (). The particular benefit of such a strategy for these cells is unclear. In a similar vein, one may wonder whether it is possible for cells to deconvolve the two pieces of information from the calcium signal using a downstream signaling network. Such “multiplexing” has been shown to be possible with simple biochemical networks deRonde2011 (), although the ways in which dynamic information is stored in, and extracted from, cellular signals is a topic of ongoing research Cheong2011 (); Wollman2014 ().
Our results suggest that the dependence of the calcium response on both sensory and collective parameters persists despite significant celltocell variability. Indeed, we have found that certain measures are robust to variability, such as the oscillation propensity and the entropy of the ISI distribution, while others are not, such as spatial correlations in the ISI and its dependence on the ATP input (frequency encoding). This implies that traditional measures of information encoding in calcium dynamics, such as frequency encoding Putney2011 (); Tang1995 (); Hofer2006 (); Woods1986 (); Meyer1991 (), may have to be rethought in contexts where celltocell variability is pronounced. It is becoming increasingly understood that variability is common in cell populations, and recent examples suggest it may even be beneficial. For example, recent studies in a related system (NFB oscillations in fibroblast populations) also found a large degree of celltocell variability Tay2010 () and demonstrated that this variability allows entrainment of the population to a wider range of inputs Tay2015 ().
In our model, the transition from the nonoscillatory to the oscillatory regime occurs due to a saddlenode bifurcation, a critical point in parameter space where the number of dynamical fixed points changes (see Appendix C). This transition is broadened by intrinsic noise and celltocell variability into a critical “region”, while cellcell communication causes the oscillation propensity to depend on cell density within this region (Fig. 3A). Nonetheless, the underlying mechanism remains the critical dynamical transition. Our finding that this region is broad, and our suggestion that it may be of some functional use for the system, resonates with recent studies that have argued that biological systems are poised near critical points in their parameter space Mora2011 (); Krotov2014 (); Hidalgo2014 (). However, the connection between dynamical criticality, as in our model, and criticality in manybody statistical systems, remains to be fully explored.
Gap junctional communications exist among many types of cells. Therefore, our results may have farreaching implications for other biological model systems such as neuronal networks or cardiovascular systems. It will be interesting to explore the extent to which distinctions in the calcium dynamics in these systems originate from differences in their degrees of cellcell communication.
V Materials and Methods
v.1 Fabrication of Microfluidic Device
See Appendix A for detailed device fabrication and characteristics.
v.2 Cell Culture and Sample Preparation
NIH 3T3 and MDAMB231 cells were cultured in standard growth mediums (Dulbeccoâs modified Eagle medium (DMEM) supplemented with 10% bovine calf serum and1% penicillin and DMEM supplemented with 10% fetal bovine serum, 1% penicillin, and 1% nonessential amino acids respectively). To prepare samples, cells were detached from culture dishes using TrypLE Select (Life Technologies) and suspended in growth mediums before pipetted into the microfluidics devices and allowed to form monolayers. If MDAMB231 cells were the dominant species (a fraction greater than 50% of all cells), they were first allowed to attach the glass bottom of the microfluidics devices. Red fluorescent tag (CellTracker, Life Technologies) was then applied and subsequently washed with growth medium. Finally, NIH 3T3 cells were injected into the device so that the desired cell density (1000 cells/mm) was reached. If NIH 3T3 cells were the dominant species, they were allowed to attach the glass bottom of the microfluidics devices first. MDAMB231 cells already loaded with CellTracker were then injected into the devices to reach the desired cell density. After incubating the microfluidics devices containing cell monolayers overnight, fluorescent calcium indicator was applied (Fluo4, Life Technologies) making the samples ready for imaging.
v.3 Fluorescence Imaging and Image Analysis
Fluorescence was detected using a inverted microscope (Leica DMI 6000B) coupled with a Hamamatsu Flash 2.8 camera. Movies were taken at a frame rate of 4 frame/sec with a 20x oil immersion objective. Image analysis and data processing were performed in MATLAB. (Details in SI).
v.4 Stochastic model of collective calcium dynamics
We extend the model in Tang1995 (); Othmer1993 () to include intrinsic noise and cellcell communication. We simulate the ensuing calcium dynamics of cells of various densities on a 3by3 grid (and 7by7 grid for Fig. 2D) using the Gillespie algorithm Gillespie1977 () with tauleaping Gillespie2001 (); Rathinam2003 (); Cao2006 (); Cao2007 ()). We include celltocell variability by sampling all model parameters from distributions that are uniform in log space, which varies foldchange up to a maximal value. For further details, see Appendix C.
Acknowledgements.
Discussions with Prof. D. Roundy and Prof. G. Schneider are gratefully acknowledged. The project is supported by NSF Grant PHY1400968. This work is also supported by a grant from the Simons Foundation (#376198 to A.M.)Appendix A Additional information of experimental setup
a.1 Fabrication and characteristics of microfluidic device
The organic elastomer polydimethylsiloxane (PDMS, Sylgard 184, DowCorning) used to create the microfluidic devices is comprised of a two part mixture  a base and curing agent  that is mixed in a 10:1 ratio, degassed, and poured over a stainless steel mold before curing at 65C overnight. Once cured, the microfluidic devices are cut from the mold, inlet/outlet holes are punched, and the device is affixed to a No. 1.5 coverslip via corona treatment. Figures 6AC provide schematics of the design of the device as well as inner dimensions of the flow chamber.
To characterize the stimuli profile in the flow chamber, we use fluorescein (Sigma Aldrich) to replace ATP as our probe and record the fluorescent intensity in the same flow condition as in the ATP sensing experiments. Figure 6D shows the normalized fluorescent intensity as a function of time from two devices. We estimate the time takes to fill the field of view with chemical stimuli is about 1 second.
a.2 Obtaining individual cell calcium dynamics
Application of CellTracker Red CMPTX (LifeTechnologies) to breast cancer cells and the calcium binding dye Alexa Fluo4 (LifeTechnologies) to both breast cancer cells and fibroblast cells within the microfluidic device allows for the differentiation of the two cell types. To do so, successive images are taken of the two fluorophores and first analyzed in ImageJ (imagej.nih.gov). Overlaying the two images generates a composite image that determines the network architecture, as is shown in Figure 1A of the main text. The centroids for each individual cell are manually determined and a dot is placed at the center of each cell as shown in Fig. 7. Using the dotted image as a mask, intensity is averaged over the 40 pixels comprising the dot for each cell and for all acquired frames in the experiment.
a.3 Quantification of gap junctional intercellular communication
In order to quantify the gap junction communications for the two cell types we have used (NIH 3T3 and MDAMB231), we perform FRAP (fluorescent recovery after photobleaching) experiments to measure the effective diffusion coefficients in a cell monolayer for each cell type. (some experimental details, eg dye used frame rates etc). As shown in Fig. 8B, at each time point , we represent the average concentration of fluorescent dye in the bleached area as and nonbleached area as . The concentration of the dye is proportional to the recorded intensity as , where is a constant. We approximate the time evolution of and based on the following assumptions: (a) there is an overall bleaching rate while imaging, affecting both and . (b) get an influx due to a concentration gradient approximately proportionally to , where is the effective diffusion coefficient. (c) the cell monolayer has a uniform thickness of . We have
(1)  
(2)  
The above equations can be simplified and we have
(3) 
In our experiment, we take image at every 30 secs, therefore is set to be 1 minute. The recorded time series , from all experiments are plugged into the above equation to obtain a distribution of the effective diffusion coefficient as shown in Fig. 8D.
Appendix B Additional information of InterSpike Interval (ISI) analysis
b.1 InterSpike Interval Determination
As calcium spike trains are unsuitable for Fourier analysis due to their varying amplitudes and irregular periods, the interspike interval (ISI) has been used as a suitable substitute for determining how the period of calcium oscillation affects information encoding within and amongst cells in the network. To determine the ISI’s from a single cell’s calcium dynamics , we first smooth the time series to remove outlying data points Garcia2010 (), to obtain . We then find all local peaks in using the following criteria: a local peak must be higher than its neighboring valleys (local minimums) on both directions by a certain value, which we refer to as the sensitivity parameter; the first and last frames are excluded as peaks or valleys. Third, we obtain ISI’s from the distances of peaks . Finally, to obtain the true ISI’s corresponding to calcium oscillations, we only keep fall in the range of 5 seconds and 150 seconds. This final step excludes about 5% of ISI’s.
Figure 9A and 9B provide two examples of the peak finding algorithm’s ability to determine peak locations for a fibroblast cell calcium profile (A) and cancer cell calcium profile (B). The sensitivity parameter used for all experiments and in Figure 9A and 9B is the standard deviation of the normalized intensity as determined empirically. Figure 9C and 9D show the algorithms diminished ability to determine peaks and to detect spurious peaks when the sensitivity parameter is set to and respectively.
b.2 Spatial map of ISI in multicellular network
We have shown a map of ISI in Fig. 5C of the main text. To generate such a map, we first calculate the ISI’s corresponding to each individual cell as shown in the above. We then calculate the average ISI of each cell. Fig. 10 shows two more examples of the ISI spatial map. Figure 10A corresponding to a map from a experiment with 20M ATP concentration and Fig. 10B corresponding to a 50M ATP concentration. Black regions within the maps indicate cells that did not oscillate in the given experiment and a red circle indicates the location of a cancer cell. It is apparent from the figure that cancer cells are more likely to be nonoscillatory than fibroblast cells. This is also shown in Fig. 11, which shows ISI distributions of both cell types for a typical experiment.
b.3 Entropy of ISI distributions
Interestingly, the network architecture also regulates the spectrum of ISI. This can be demonstrated from the perspective of each individual cell’s microenvironment. We have grouped NIH 3T3 cells based on their number of nearest neighbors (), as well as the fraction of cancer cells in their neighbors. For each group at a fixed ATP concentration (50 M), we have computed the differential entropy of interspike intervals . The differential entropy is computed using the knearest neighbor method, as described in Appendix D. As shown in Fig. 12, is higher for cells that are surrounded by a large number of cancer cells, and lower for cells with pure fibroblast neighbors. In the latter case, also increases as the number of nearest neighbors decreases. These results show that the spectrum of ISI is regulated by cellcell communications both locally in the immediate microenvironment of a cell and globally throughout the whole defective multicellular network.
b.4 Statistics of , the fraction of nonoscillating cells
As shown in Figs. 4 and 5 of the main text, the fraction of nonoscillating cells systematically depends on the network architecture. Fig. 13 shows pairwise statistical comparison of the results. The results are obtained by ANOVA followed with Fisher’s least significant difference procedure, implemented in MATLAB.
Appendix C Theoretical Model
c.1 Deterministic Model
We model the dynamics of calcium release following Tang and Othmer Tang1995 (). Their model is deterministic and valid for a single cell. First we describe the behavior of their model. Then we generalize it to include stochastic effects. Next, we generalize it further to include cellcell communication. Finally, we add celltocell variability. This procedure allows us to systematically explore the effects of stochasticity, communication, and variability on the calcium dynamics.
ATP triggers the production of IP within a cell. IP binds to receptors on the endoplasmic reticulum (ER). The ER is an internal store for calcium ions. When the receptors are active, they act as channels that allow passage of calcium between the ER and the cytoplasm of the cell. Receptor activation requires binding of both IP and cytoplasmic calcium, such that calcium triggers its own release. However, an excess of calcium deactivates receptors, closing the channels. These positive and negative feedbacks, respectively, lead to oscillations and other rich dynamics.
The model of Tang and Othmer is based upon the following reactions,
(4)  
(5) 
Eqn. 4 contains the reactions describing calcium transport between the ER and the cytoplasm, while Eqn. 5 contains the reactions describing the sequential (un)binding of molecules to receptors. In Eqn. 4, the first reaction describes leakage between the calcium in the ER () and calcium in the cytoplasm (); is the ratio of the ER volume to the cytoplasm volume . The second reaction describes transport of calcium through the receptor channels: receptors with both IP3 and calcium bound () are active channels. The last reaction describes the pumping of calcium from the cytoplasm to the ER; denotes the nonlinear pumping propensity. In Eqn. 5, the first reaction describes the binding of IP3 to bare receptors (), creating the complex ; is the concentration of IP3, which is treated as a parameter in their model (we later generalize this feature, making a dynamic variable and its production rate the control parameter, since this is more consistent with the experimental setup). The second reaction describes the subsequent binding of cytoplasmic calcium to the complex, which activates the complex (). The last reaction describes the further binding of cytoplasmic calcium to another site on the complex, which deactivates the complex (). Only receptors in the state serve as active channels.
The deterministic model describes the dynamics of the mean numbers of molecules. The dynamics follow directly from Eqns. 4 and 5. Introducing and as the numbers of calcium molecules in the cytoplasm and ER, respectively, and as the numbers of receptors in state , the means (denoted with an overbar) obey
(7)  
(8)  
(9)  
(10)  
(11) 
The pumping propensity is a Hill function with coefficient and parameters and .
Two quantities are conserved in the model: the total number of calcium molecules and the total number of receptors ,
(12)  
(13) 
The conservation relations reduce the number of dynamical equations from six to four. Specifically, we eliminate and . Furthermore, we assume that the number of calcium molecules is much larger than the number of receptors, . In this limit, introducing as the concentration of calcium in the cytoplasm (both unbound and bound to receptors) and as the fraction of receptors in each state, Eqns. C.111 become
(14)  
(15)  
(16)  
(17) 
where , , , , , , and is the average calcium concentration in both the cytoplasm and the ER. Eqns. 1417 reproduce Eqn. 4 in Tang1995 (). The parameter values used in Tang1995 (), which we also use here, are given in Table 1.
Deterministic parameters  Stochastic parameters 

m  
M  
s  s 
s  s 
Ms  s 
M  
Ms  s 
Ms  s 
Ms  s 
s  s 
s  s 
s  s 
(varied)  
s  
s 
Tang and Othmer recognize that the (un)binding of calcium to the deactivating site is slow compared to the other (un)binding reactions, and . This permits the quasisteadystate approximations and , which reduce Eqns. 1417 to two equations. In dimensionless form they read
(18)  
(19) 
where , , , , , , , , , and .
Eqns. 18 and 19 exhibit four dynamic regimes. As shown in Fig. 14, the system can be (A) monostable, with low cytoplasmic calcium, (B) excitable, (C) oscillatory, or (D) monostable, with high cytoplasmic calcium. In the last regime (D), near the boundary, oscillations can also be supported for certain initial conditions. All four regimes are accessible by tuning the IP concentration (Fig. 14E). Transitions between regimes occur at the critical concentrations , , and . Of particular interest is the transition from the excitable to the oscillatory regime (), since these are the dynamics exhibited by cells in the experiments: transient pulsing or sustained oscillations.
c.2 Stochastic Model
Stochastic dynamics are simulated using an adaptive tauleaping method Gillespie2001 (); Rathinam2003 (); Cao2006 (); Cao2007 (), which is a computationallyefficient approximation of the Gillespie algorithm Gillespie1977 (). The simulation accounts for the facts that molecule numbers are integervalued and that reactions fire at random, exponentiallydistributed times. The reaction propensities follow directly from Eqns. 45, and have the same forms as the individual terms in Eqns. C.111 (without the overbars). The stochastic model parameters are related to the deterministic model parameters by the cytoplasmic volume and the total number of receptors ; see Table 1. Together these factors determine the noise in the system via the molecule numbers and . Since relative fluctuations decrease with molecule number, large corresponds to low noise, while small corresponds to high noise. Typical trajectories are shown in Fig. 15.
We estimate the cell volume from the experiments. An upper bound is obtained from the measurements of the cell density , since the cell area cannot be larger than . The maximal density mm gives the sharpest upper bound, for an area of m, or a length scale of m. Since this is an upper bound, we assume a length scale range of m. For computational efficiency, we perform simulations assuming m, and taking a monolayer width of m, this leads to a cell volume of m. Absent an experimental estimate for the number of receptors , we take .
In the stochastic model, we no longer treat (the concentration of IP3) as a parameter; instead, we treat the number of IP3 molecules as a stochastic variable. The control parameter is then the production rate of IP3, , which is assumed to be proportional to the ATP concentration in the experiments. IP3 is degraded at a rate , meaning that the stochastic reactions simulated are Eqns. 45, with the first reaction of Eqn. 5 replaced with
(20) 
The degradation rate of IP3 has been measured to be s Wang1995 (), and we vary in the range s such that the mean IP3 level is molecules per cell, or M.
c.3 CelltoCell Communication
We introduce communication into the system by allowing calcium molecules to diffuse between neighboring cells, simulating communication via gap junctions in the experiments. Doing so introduces new reactions into our model:
(21) 
where index cell lattice sites, and is the hopping rate of calcium ions from cell to cell. We estimate from the experimentally measured diffusion coefficient m/s (see Fig. 1D of the main text). Since in our model, ions are hopping between cells on a twodimensional lattice, corresponds to the expected time to diffuse a celltocell distance . Using m as above, we find s, and therefore we use s.
Our calculations are done on a square lattice, where each lattice site is either empty or contains one cell, and therefore each cell can have up to four neighbors with which to communicate. Density is varied by changing the number of empty lattice sites. For each chosen density, we sample over individual realizations, in which cell locations are assigned randomly. Thus, the statistics encompass many possible spatial distributions of cells for each density.
c.4 Effects of noise and communication
Once we add noise, we can no longer distinguish the dynamic regimes (monostable, excitable, oscillatory) using the deterministic criterion, namely the number and stability of the fixed points (Fig. 14E). Instead, we define a criterion based on the trajectories themselves, similar to the procedure for the experimental trajectories. We simulate the system for a given time , where s as in the experiments. The criterion is then:

If the calcium level spikes zero times in , the dynamics are monostable (nonoscillatory),

If the calcium level spikes one time in , the dynamics are excitable (nonoscillatory), and

If the calcium level spikes two or more times in , the dynamics are oscillatory.
When analyzing trajectories to determine ISI values, excitations separated by less than s are excluded, as are excitations separated by more than s, just as in the experimental analysis. We validate the above criterion using a complementary approach: we also compute the power spectrum of the trajectories. For periodic trajectories, the power spectrum exhibits a peak at the oscillation frequency. As shown in Fig. 15BD, trajectories that are periodic, and thus pass the above criterion, indeed given peaked power spectra as expected.
We find that noise broadens the boundaries between dynamical regimes (Fig. 2B of the main text), and that noise shifts the boundaries between regimes (Fig. 15A, compare the bars in the middle and bottom rows). In particular, the oscillatory regime is expanded at both ends. At low values of the control parameter, noise promotes oscillations where the deterministic dynamics are excitable. At high values of the control parameter, noise promotes oscillations where the deterministic dynamics are monostable. Both effects have been observed in stochastic models of similar excitable systems, and are likely due to the ability of noise to cause repeated excitations and prevent damping of oscillations, respectively Mugler2016 ().
c.5 CelltoCell Variability
The model predicts frequency encoding, but experimental measurements show no trend in ISI values as a function of ATP concentration (Fig. 1E). We thus introduce celltocell variability into the model as a possible explanation for the lack of observed frequency encoding.
Variability is introduced into our model by allowing parameters for each cell to be drawn from a distribution. All parameters are allowed to vary from their baseline values listed in Table 1. Parameters are sampled from uniform distributions in log space such that loglog,log, where is the mean value for the parameter , and is a parameter that changes the strength of variability. Therefore each parameter is sampled from a distribution such that .
Fig. 16 shows the ratio of the standard deviation in mean ISI values over their mean as a function of . As is increased, this ratio increases until saturation occurs at . Also shown is box plot of the experimental values of this ratio, which was obtained from twelve experiments with similar density. By the saturating value , the theoretical values lie within the experimental range; therefore we use the value to obtain the results that include variability in the main text.
Appendix D Nearest neighbor crosscorrelation analysis and computation of differential entropy
d.1 Definition
We have previously demonstrated the effectiveness of using nearest neighbor crosscorrelation functions () to characterize the calcium dynamics of collective chemosensing Sun2012 (). As a result, we first study how network architecture may affect the spectrum of . For any cell , we represent its calcium dynamics (fluorescent intensity profile) as , and define as:
(22) 
To evaluate , we numerically differentiate the response curve using the fivepoint stencil method; is the time average. Note that is the variance of , which normalizes to be dimensionless. Based on Delauney triangulation of the multicellular network, topological distance can be defined for each cell pair where for nearest neighbors. The mean nearest neighbor crosscorrelation function is obtained by averaging over all nearest neighbor pairs .
d.2 Statistics of nearest neighbor pairs
[ATP]  10 M  20 M  50 M  100 M  Total 

# of cells  7,055  4,245  13,078  5,062  30000 
# of nearest neighbor pairs  40,596  24,462  75,092  32,222  172,372 
By using Delauney triangulation coupled with the manually determined cell centers (see Fig. 7 for example), we are able to determine the connectedness of the monolayer network. In particular, we focus on the communication between nearest neighbor cells to reflect proper intercellular communication which correspond to a topological distance of in the triangulation. The total number of nearest neighbor pairs and cells analyzed for each ATP concentration are given in table 2.
d.3 of collective calcium dynamics
We compute the spectrum of using differential entropy, a scalar value that characterizes the randomness of a set of observables. The differential entropy is defined as , where is the probability density of a continuous variable . Higher values of differential entropy are associated with broader distributions; this allows us to directly compare the spectrum of nearest neighbor crosscorrelations. To avoid biasing the entropy calculation due to binning of data, the differential entropy is determined by using the nearest neighbor algorithm selim (); Loft (), as described below. As shown in Fig. 17, increasing cell density decreases the differential entropy of , suggesting that higher levels of network connectivity will suppress the fluctuations in the nearest neighbor crosscorrelations.
d.4 Evaluations of differential entropies
In order to characterize the spectrum of nearest neighbor crosscorrelations and inter spike intervals of the collective calcium dynamics, we have implemented a knearest neighbor method to calculate the differential entropy of a set of random variables. We first make benchmark tests using random variables drawn from Gaussian distributions. Since most experimental datasets contain more than 1500 elements, we compare the differential entropy calculated from 1500 Gaussian random variables with respect to the analytic results. As shown in Fig. 18A, numerical calculations agree with exact values very well. The performance of knearest neighbor method depends on the sample size. To evaluate the dependence, we employ knearest neighbor method on random variables drawn from standard normal distribution with varying sample sizes. As shown in Fig. 18B, the relative error is around 7% for sample size equals to 100, and quickly decay to less than 1% when sample size is greater than 500.
In calculating the differential entropy using knearest neighbor method, every element in a dataset contribute to the final result. In order to estimate how much fluctuates as a result of resampling, we would like to calculate from many datasets each of which contains distinct elements drawn from the same probability distribution. However, the distribution of experiment observables are unknown. In order to perform resampling, we employ statistical bootstrap Duval1993 (). As a particular example, we take a randomly selected set of nearest neighbor crosscorrelations and calculated the differential entropy of 1000 bootstrap resamplings. Fig. 18C shows the probability distribution of the results. As can be seen, the fluctuation is small, and is typically less than 2%.
References
 (1) Lauffenburger, DA (2000) Cell signaling pathways as control modules: Complexity for simplicity? Proc. Natl. Acad. Sci. 97:5031–5033.
 (2) Lodish, H et al. (2000) Molecular Cell Biology, 4th edition (W. H. Freeman, New York, United States).
 (3) Barritt, G (1994) Communications Within Animal Cells (Oxford Science Publications, Oxford, UK).
 (4) Miller, MB, Bassler, BL (2001) Quorum sensing in bacteria. Annu. Rev. Microbiol. 55:165–199.
 (5) Gregor, T, Fujimoto, K, Masaki, N, Sawai, S (2010) The onset of collective behavior in social amoebae. Science 328:1021–1025.
 (6) Jones, WD, Cayirlioglu, P, Kadow, IG, Vosshall, LB (2007) Two chemosensory receptors together mediate carbon dioxide detection in drosophila. Nature 445:86–90.
 (7) Smear, M, Shusterman, R, O’Connor, R, Bozza, T, Rinberg, D (2011) Perception of sniff phase in mouse olfaction. Nature 479:397–400.
 (8) Benninger, RK, Zhang, M, Head, WS, Satin, LS, Piston, DW (2008) Gap junction coupling and calcium waves in the pancreatic islet. Biophysical journal 95:5048–5061.
 (9) Schneidman, E, II, MJB, Segev, R, Bialek, W (2006) Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440:1007–1012.
 (10) Greschner, M et al. (2011) Correlated firing among major ganglion cell types in primate retina. J. Physiol. 589:75–86.
 (11) Sun, B, Lembong, J, Normand, V, Rogers, M, Stone, HA (2012) The spatialtemporal dynamics of collective chemosensing. Proc. Nat. Aca. Sci 109:7759–7764.
 (12) Sun, B, Doclos, G, Stone, HA (2013) Network characterization of collective chemosensing. Phys. Rev. Lett. 110:158103.
 (13) Kumar, NM, Gilula, NB (1996) The gap junction communication channel. Cell 84:381–388.
 (14) Berridge, MJ, Irvine, RF (1989) Inositol phosphates and cell signalling. Nature 341:197–205.
 (15) CD, CDF, Huganir, HR, Snyderc, SH (1990) Calcium flux mediated by purified inositol 1,4,5trisphosphate receptor in reconstituted lipid vesicles is allosterically regulated by adenine nucleotides. Proc. Natl. Acad. Sci. 87:2147–2151.
 (16) Dupont, G, Combettes, L, Bird, GS, Putney, JW (2011) Calcium oscillations. Cold Spring Harb Perspect Biol. 3:a004226.
 (17) Loewenstein, WR, Kanno, Y (1966) Intercellular communication and the control of tissue growth: lack of communication between cancer cells. Nature 209:1248–1249.
 (18) McLachlan, E, Shao, Q, Wang, HL, Langlois, S, Laird, DW (2006) Connexins act as tumor suppressors in threedimensional mammary cell organoids by regulating differentiation and angiogenesis. Cancer Res. 66:9886–9894.
 (19) Zhou, JZ, Jiang, JX (2014) Gap junction and hemichannelindependent actions of connexins on cell and tissue functions: An update. FEBS Lett. 588:1186–1192.
 (20) Tang, Y, Othmer, HG (1995) Frequency encoding in excitable systems with applications to calcium oscillations. Proceedings of the National Academy of Sciences 92:7869–7873.
 (21) Thurley, K et al. (2014) Reliable encoding of stimulus intensities within random sequences of intracellular Ca spikes. Science Signaling 7:ra59.
 (22) Woods, N, Cuthbertson, K, Cobbold, P (1985) Repetitive transient rises in cytoplasmic free calcium in hormonestimulated hepatocytes. Nature 319:600–602.
 (23) Meyer, T, Stryer, L (1991) Calcium spiking. Annual review of biophysics and biophysical chemistry 20:153–174.
 (24) Othmer, HG, Tang, Y (1993) in Experimental and theoretical advances in biological pattern formation, eds Othmer, HG, Maini, PK, Murray, JD (Springer US).
 (25) Meyer, T, Stryer, L (1988) Molecular model for receptorstimulated calcium spiking. Proc. Natl. Acad. Sci. 85:5051–5055.
 (26) Politi, A, Gaspers, LD, Thomas, AP, Höfer, T (2006) Models of ip3 and ca2+ oscillations: Frequency encoding and identification of underlying feedbacks. Biophys. J. 90:3120–3133.
 (27) Swillens, S, Dupont, G, Combettes, L, Champeil, P (1999) From calcium blips to calcium puffs: Theoretical analysis of the requirements for interchannel communication. Proc. Natl. Acad. Sci. 96:13750–13755.
 (28) Falcke, M (2004) Reading the patterns in living cells the physics of ca2 signaling. Adv. Physics 53:255–440.
 (29) Gillespie, DT (1977) Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry 81:2340–2361.
 (30) Meister, M (1996) Multineuronal codes in retinal signaling. Proc. Natl. Acad. Sci. 93:609–614.
 (31) Schnizer, MJ, Meister, M (2003) Multineuronal firing patterns in the signal from eye to brain. Neuron 37:499–511.
 (32) Janmey, PA, Miller, RT (2011) Mechanisms of mechanical signaling in development and disease. Journal of cell science 124:9–18.
 (33) Abbaci, M, BarberiHeyob, M, Blondel, W, Guillemin, F, Didelon, J (2008) Advantages and limitations of commonly used methods to assay the molecular permeability of gap junctional intercellular communication. BioTechniques 45:33–62.
 (34) Schwiebert, EM, Zsembery, A (2003) Extracellular atp as a signaling molecule for epithelial cells. Biochimica et Biophysica Acta (BBA)Biomembranes 1615:7–32.
 (35) Falzoni, S, Donvito, G, Di Virgilio, F (2013) Detecting adenosine triphosphate in the pericellular space. Interface Focus 3:20120101.
 (36) Trabanelli, S et al. (2012) Extracellular atp exerts opposite effects on activated and regulatory cd4+ t cells via purinergic p2 receptor activation. The Journal of Immunology 189:1303–1310.
 (37) Pellegatti, P et al. (2008) Increased level of extracellular atp at tumor sites: in vivo imaging with plasma membrane luciferase. PloS one 3:e2599.
 (38) Mehta, P, Gregor, T (2010) Approaching the molecular origins of collective dynamics in oscillating cell populations. Current opinion in genetics & development 20:574–580.
 (39) De Monte, S, d’Ovidio, F, Dano, S, Sorensen, PG (2007) Dynamical quorum sensing: Population density encoded in cellular dynamics. Proceedings of the National Academy of Sciences 104:18377–18381.
 (40) Taylor, AF, Tinsley, MR, Wang, F, Huang, Z, Showalter, K (2009) Dynamical quorum sensing and synchronization in large populations of chemical oscillators. Science 323:614–617.
 (41) de Ronde, W, Tostevin, F, Ten Wolde, PR (2011) Multiplexing biochemical signals. Physical review letters 107:048101.
 (42) Cheong, R, Rhee, A, Wang, CJ, Nemenman, I, Levchenko, A (2011) Information transduction capacity of noisy biochemical signaling networks. science 334:354–358.
 (43) Selimkhanov, J et al. (19652014) Accurate information transmission through dynamic biochemical signaling networks. Science 346:1370–1373.
 (44) Tay, S et al. (2010) Singlecell nf[kgr] b dynamics reveal digital activation and analogue information processing. Nature 466:267–271.
 (45) Kellogg, RA, Tay, S (2015) Noise facilitates transcriptional control under dynamic inputs. Cell 160:381–392.
 (46) Mora, T, Bialek, W (2011) Are biological systems poised at criticality? Journal of Statistical Physics 144:268–302.
 (47) Krotov, D, Dubuis, JO, Gregor, T, Bialek, W (2014) Morphogenesis at criticality. Proceedings of the National Academy of Sciences 111:3683–3688.
 (48) Hidalgo, J et al. (2014) Informationbased fitness and the emergence of criticality in living systems. Proceedings of the National Academy of Sciences 111:10095–10100.
 (49) Gillespie, DT (2001) Approximate accelerated stochastic simulation of chemically reacting systems. The Journal of Chemical Physics 115:1716–1733.
 (50) Rathinam, M, Petzold, LR, Cao, Y, Gillespie, DT (2003) Stiffness in stochastic chemically reacting systems: The implicit tauleaping method. The Journal of Chemical Physics 119:12784–12794.
 (51) Cao, Y, Gillespie, DT, Petzold, LR (2006) Efficient step size selection for the tauleaping simulation method. The Journal of chemical physics 124:044109.
 (52) Cao, Y, Gillespie, DT, Petzold, LR (2007) Adaptive explicitimplicit tauleaping method with automatic tau selection. The Journal of chemical physics 126:224101.
 (53) Garcia, D (2010) Robust smoothing of gridded data in one and higher dimensions with missing values. Computational statistics & data analysis 54:1167–1178.
 (54) Wang, S, Alousi, AA, Thompson, SH (1995) The lifetime of inositol 1, 4, 5trisphosphate in single cells. The Journal of general physiology 105:149–171.
 (55) Mugler, A et al. (2016) Noise expands the response range of the bacillus subtilis competence circuit. PLoS Computational Biology Accepted.
 (56) Selimkhanov, J et al. (2014) Accurate information transmission through dynamic biochemical signaling networks. Science 346:1370–1373.
 (57) Loftsgaarden, DO, Quesenberry, CP et al. (1965) A nonparametric estimate of a multivariate density function. The Annals of Mathematical Statistics 36:1049–1051.
 (58) Mooney, CZ, Duval, RD, Duval, R (1993) Bootstrapping: A nonparametric approach to statistical inference (Sage) No. 9495.