Paracrine communication maximizes cellular response fidelity in wound signaling

Population averaging due to paracrine communication can arbitrarily reduce cellular response variability. Yet, variability is ubiquitously observed, suggesting limits to paracrine averaging. It remains unclear whether and how biological systems may be affected by such limits of paracrine signaling. To address this question, we quantify the signal and noise of Ca2+ and ERK spatial gradients in response to an in vitro wound within a novel microfluidics-based device. We find that while paracrine communication reduces gradient noise, it also reduces the gradient magnitude. Accordingly we predict the existence of a maximum gradient signal to noise ratio. Direct in vitro measurement of paracrine communication verifies these predictions and reveals that cells utilize optimal levels of paracrine signaling to maximize the accuracy of gradient-based positional information. Our results demonstrate the limits of population averaging and show the inherent tradeoff in utilizing paracrine communication to regulate cellular response fidelity. DOI: http://dx.doi.org/10.7554/eLife.09652.001


Introduction
Cellular variability is likely a biological trait with significant phenotypic consequences. Technological advances in single-cell measurement methodologies reveal substantial cellular variability. For instance, single-cell quantification of protein concentration variability between cells shows that the concentration of many signaling molecules can vary by~25% (coefficient of variation) (Sigal et al., 2006;Bar-Even et al., 2006;Spencer et al., 2009). Furthermore, a large and rapidly growing body of single-cell transcriptomics experiments further demonstrates that cells homogeneous in 'type' have substantially heterogeneous gene expression patterns (Junker and Van Oudenaarden, 2014). The origin of this cellular variability has been traced to fundamental properties of gene expression. Notably, single-molecule kinetics regulates gene expression and, as a result, is an inherently stochastic process (Sanchez and Golding, 2013).
While the costs and benefits of cellular variability are likely dependent on the specific physiological context, the functional significance of cellular variability suggests that cellular variability magnitude is regulated. Functional analysis of cellular response variability demonstrates that the observed cellular variability affects the core function of signaling networks. Despite a homogenous environment, cells respond in a heterogeneous manner due to biological variability. Response variability potentially degrades transmitted information and decreases downstream effector ability to reliably respond to environmental changes (Selimkhanov et al., 2014;Cheong et al., 2011;Voliotis et al., 2014;Hansen and O'Shea, 2015). The abundance of cellular variability throughout biological processes and the potential consequences of information degradation suggest that biological systems have developed mechanisms to regulate cellular variability. However, cellular variability is not Figure 1. Local averaging using paracrine signaling reduces response variability in a communication distance dependent manner. (A) A hypothesis for local averaging based reduction of response variability using paracrine signaling for ERK activation by P2Y receptors. ATP binds to P2Y receptors to increase cytosolic Ca 2+ levels with high variability between cells despite equivalent ATP dosage per cell (pink shading). EGF release from each cell is proportional to the primary response to ATP (green arrows). Due to diffusion of EGF, the local concentration will be the average of EGF released from nearby cells and subject cells in a local neighborhood to the same level of EGF (blue arrows) to result in similar ERK activation. (B) ATP activates ERK in a paracrine fashion in MCF-10A cells. ERK response to 10 mM ATP addition with (green dashed) and without (green solid) 1 mM of EGFR inhibitor tryphosphitin AG1478. 0 mM ATP addition with (blue dashed) and without (blue solid) 1 mM AG1478 shown as controls. (C) MCF-10A cells were clustered based on their spatial proximity so that cells within a cluster were within a specific maximal communication distance (MCD) and cells in other clusters were farther than the MCD. Cluster denisty was calcuated by dividing the number of cells per cluster by the circular area inhabitated by each cluster. (D) Standard deviation of ERK response per cluster to 10 mM ATP. Standard deviation decreases with increasing cluster density (p-value <0.05, Pearson correlation). (E) Average maximum Ca 2+ response with increasing ATP dosage. The standard deviation, that is noise, of the Ca 2+ response to each ATP dosage is large when compared to the increase in average response with increasing ATP dose, that is signal (SNR = 0.9, gray error bars [standard deviation]). Noise decreases when Ca 2+ response is locally averaged with a PCD of 100 mm (SNR = 20.2, black error bars [standard deviation]).
x-axis shown in log-scale. F. Single-cell maximum Ca 2+ response is locally averaged across an area specified by the PCD to produce a predicted single-cell ERK response. Variability between cells in the predicted ERK response decreases with increasing PCD. Response magnitude of Ca 2+ and ERK response indicated by pink to yellow and blue to yellow colorbars, respectively. G. The SNR of the predicted ERK response from locally averaged Ca 2+ data continually increases with increasing PCD shown for a model with rapid diffusion (green) or limited by the diffusion rates and integration time of paracrine signals (blue, Materials and methods). SNR calculated in same manner as panel E with increasing PCD. Shaded area is SEM (N = 5). DOI: 10.7554/eLife.09652.003 The following figure supplements are available for Figure 1: The initial paracrine signaling pathways that are activated in response to Damage Associated Molecular Patterns (DAMPs) are a good model system for investigating the influence and limits of paracrine communication on cellular response fidelity. Paracrine communication is pervasive during initial wound response. Wound healing begins as soon as the wound occurs and the initial cellular wound response provides the foundation for proper downstream healing. The initial cellular wound response relies on external environmental cues as well as programs inherent to the cell, including DAMPs as primary danger signals (Enyedi and Niethammer, 2015). DAMPs are released from necrotic cells and bind to extracellular receptors on surrounding cells. This binding initiates a signal in the surrounding cells to secrete a secondary set of cytokines and growth factors required to coordinate the wound healing process. Many DAMP signals, such as extracellular ATP, are transient and released in limited quantities. As a result, the initial wound response to such DAMPs shows high cellular variability and low fidelity. Despite the limited fidelity of the initial wound response, the wounded epithelium is able to establish a robust healing response. The complicated and multi-step wound healing process utilizes several paracrine communication mechanisms to share cellular information and coordinate the overall healing program.
Here we use the paracrine release of epidermal growth factor (EGF) ligands initiated by ATP binding to P2Y receptors as a model to investigate the limits of cellular information sharing through paracrine communication to mitigate biochemical noise ( Figure 1A). We show that paracrine communication increases extracellular signal-regulated kinase (ERK) response fidelity using live single-cell quantitative fluorescent imaging of primary Ca 2+ and secondary ERK responses downstream of P2YR and EGFR, respectively. Statistical analysis of the primary response signal-to-noise ratio (SNR) demonstrates that the increase in response fidelity is limited by paracrine communication distance (PCD). To analyze this pathway in the physiological context of wound response we developed a new microfluidics device to monitor the spatial propagation of initial wound response signaling. Our results demonstrate that the interplay between the wound induced spatial signaling gradient and the cellular noise pattern produces an optimal PCD. The optimal PCD balances the benefits of decreased noise from local averaging with the cost of reduced signal of the spatial signaling gradient due to over-averaging. Empirical measurements of the PCD reveal that cellular communication occurs at a distance to maximize cellular response fidelity.

Paracrine signaling reduces response variability
Here we establish that the paracrine activation of ERK by ATP provides a suitable system to investigate signaling response fidelity changes due to paracrine communication. Extracellular ATP binding to P2YR results in EGF family ligand release to bind EGFR and activate ERK response as monitored by ERK activity following ATP addition. In the mammary epithelial cell line MCF-10A, addition of extracellular ATP increases ERK kinase activity in an EGFR dependent manner ( Figure 1B) similar to results reported in other in vitro epithelial models (Yin et al., 2007). ERK, as measured by the genetically encoded FRET sensor EKAREV (Albeck et al., 2013;Komatsu et al., 2011), increases when stimulated with ATP. Inhibiting EGFR with tryphostin AG1478 prevents ERK activation upon ATP addition showing that ERK activation depends on secreted EGF binding to EGFR (Wetzker and Bö hmer, 2003).
With our paracrine communication system established, we next confirmed the influence of paracrine communication on cellular response variability. Under conditions of low cell density we used a spatial clustering analysis to group cells such that the distance between groups effectively constrained communication to cells within groups ( Figure 1C,D, Figure 1-figure supplement 1,2). In the case that cellular coordination is beneficial, we anticipated that the ERK response within groups of higher cellular density, that is cells have increased communication ability, would have reduced response variability than cell clusters with decreased communication ability. ERK response variability within clusters decreases with increasing cluster density indicating increased intercellular communication ability. Increased intercellular communication ability is not observed in the absence of paracrine communication such as with the primary Ca 2+ response to ATP ( Figure 1D  . Together these observations support the hypothesis that paracrine communication decreases cellular variability by increasing cellular coordination.

Computational analysis of variability reduction resulting from paracrine information sharing
Next we developed a computational model that mimics the coordination-effects of paracrine communication (see Materials and methods for details). This computational model quantifies the overall observed benefit of paracrine coordination and predicts the potential reduction in variability. Experimental single-cell dose response data of primary Ca 2+ (prior to paracrine communication) response to ATP is used as an input to predict the secondary ERK response. We quantify cellular response fidelity by using a simple signal-to-noise analysis (SNR). In this analysis the cellular response magnitude of the input ligand (signal) is divided by the cellular response variability (noise). The signal is estimated by calculating the spread between the average cellular Ca 2+ responses from multiple ATP concentrations using multi-well dose response data. Noise is calculated from the average variability between cellular responses to a single input ligand concentration. SNR is simply the ratio of these two estimates ( Figure 1E). To mimic the benefit of paracrine communication our computational model performs a local, spatially weighted average (convolution) of the primary Ca 2+ response to predict the variability of response post paracrine communication (ERK) ( Figure 1F). In short, the convolution averages the signal for every cell with its associated surrounding cells by weighting the surrounding cells based on a Gaussian function parameterized with varying PCDs. The PCD represents how far the paracrine molecule travels from a single-cell to activate its associated surrounding cells. Local spatial averaging provides an upper bound of the possible benefit resulting from cellular communication in conditions where no additional noise exists in the paracrine pathway. This analysis indicates that paracrine averaging using a PCD of 100 mm increases response SNR from 0.9 to 20.2 by decreasing noise, or response variability, of the predicted ERK response ( Figure 1E, gray/black). To investigate the limits of paracrine averaging, we repeated this analysis for multiple PCDs. Interestingly, our analysis estimates that the overall response SNR can increase up to eightyfold at PCDs of 500 mm when paracrine diffusion is not limiting, and up to twenty fivefold when diffusion of the paracrine ligand is limiting ( Figure 1G). More sophisticated statistical measures, such as mutual information, produce similar results ( Figure 1G, Figure 1-figure supplement 5). The large maximal SNR benefit suggests a potentially noise-free ERK response. However, experimental measurements of ERK response fidelity shows substantial ERK variability indicating potential factors that limit the benefit gained from paracrine communication (data not shown).
Cellular response fidelity depends on the extent of paracrine signaling during wound response Extracellular ATP released from necrotic cells act as DAMPs to activate healthy cells proximal to the wound (Cordeiro and Jacinto, 2013). Given this role, the spatial component produced by the ATP concentration gradient and the resulting cellular positional information relative to the wound may be important in the analysis of paracrine communication that occurs over hundreds of microns from the wound. Our previous SNR analysis demonstrating an increasing SNR with increasing communication distance was done based on multi-well experiment data. However, the bolus addition of ATP creates a spatially uniform ligand concentration in the well and does not represent a physiologically relevant spatial component. To examine whether ATP spatial patterns influence the paracrine communication benefit we repeated the SNR analysis using single-cell wound response data.
In order to measure the spatial wound response for epithelial cells, we first developed a convection-free, small-volume wounding device. Scratch-assays, where a monolayer of cells is mechanically wounded using a pipet tip, are traditionally used for epithelial cell wounding (Sholley et al., 1977). Although the scratch-assay is useful for studying cell-migration following wounding, scratch-assays lack the ability to study paracrine signaling. The large volume above the cells and convection caused by the scratch present challenges to examine paracrine signaling due to the dilution and inadvertent mixing of any paracrine molecules released from a cell into the surrounding media. To circumvent these technical issues we developed a microfluidics based wounding device (Figure 2A,B). Our device has two components: an air channel (black) and a cell chamber with a~2.5 mL volume (orange). The ceiling of the cell chamber has a PDMS pillar that, when air pressure is increased in the upper air channel, lowers down on to the cells, thereby wounding the cells in the cell chamber in a highly controlled and reproducible manner ( Figure   Increasing the air pressure in the air channel lowers a pillar in the ceiling of the cell chamber until cells below the pillar are mechanically crushed (middle). The pillar returns to the original height when air pressure is released (bottom). (C) Ca 2+ response visualized by the Fluo-4 Ca 2+ indicator dye over a period of 5 min following a 300 mm diameter wound (black circle). (D) Top: Maximum single-cell (dots) Ca 2+ response to a 300 mm wound. Inset shows maximum Ca 2+ response to 300 mm wound in the presence of the ATP scavenger apyrase. Bottom: Maximum single-cell (dots) of Ca 2+ response to wound according to distance from the wound. (E) Same as D but for maximum ERK response. (F) Top: Cells are binned according to distance from the wound ( Figure 3A) and the average and standard deviation (error bars) are found for each bin. Middle: Coefficient of variation (CV) calcuated by dividing the standard deviation of each bin by the mean of that bin. Bottom: Ca 2+ has higher variability than ERK response for the wound according to the CV of every bin for all wounds (Black bar = average CV, p-value by t-test).  We used our wounding device to monitor Ca 2+ and ERK response to a 300 mm diameter wound using a stable, dual reporter MCF-10A cell line expressing the genetically encoded Ca 2+ indicator RGECO (Akerboom et al., 2013;Zhao et al., 2011) and the EKAREV FRET reporter for ERK (Albeck et al., 2013;Komatsu et al., 2011) (Figure 2D,E; Video 2). We verified the key role of ATP in initial wound response by wounding in the presence of apyrase, an enzyme that rapidly hydrolyzes ATP. Wounding in the presence of apyrase inhibits both Ca 2+ and ERK response ( Figure 2D,E, insets). From each wound we quantified single-cell time traces for over 3000 cells. Notably, the maximum activity per cell shows a larger response in cells closer to the wound compared to cells farther away from the wound for both Ca 2+ and ERK. These response gradients demonstrate the importance of the cellular position to determine the cellular response, or positional information ( Figure 2D,E). We used coefficient of variation (CV) to measure the variability of the post-paracrine ERK response and the pre-paracrine Ca 2+ response in the wound ( Figure 2F). Indeed, the CV for Ca 2+ wound response shows statistically higher variability than ERK wound response indicating that paracrine communication reduces response variability during initial wound response.
We adapted the computational SNR analysis to wound response data to determine the influence of spatial patterns on response fidelity. As opposed to the dose-response data, the wound response data uses the distance of each cell from the wound as the input rather than the concentration of activating ligand ( Figure 3A). Similar to the dose-response data, noise is estimated by averaging the cellular response variability over all distances. The variability between the average response magnitude of each distance constituted the signal ( Figure 3B). Other statistical measures of response fidelity such as mutual information were also adapted for the wound context (Figure 3-figure supplement 1).
The maximum primary Ca 2+ response shows highly variable cellular response when plotted according to distance ( Figure 3C, pink). This variability complicates the ability for a cell to distinguish its respective position to the wound based on its response We again mimicked paracrine communication to predict the post-paracrine ERK response by locally averaging the single-cell Ca 2+ wound response using a Gaussian kernel ( Figure 1F). Locally averaging the cellular Ca 2+ response creates a smoother predicted ERK response pattern versus distance from the wound ( Figure 3C, gray). However, the reduction in variability also decreases the overall response pattern trend. Locally averaging the Ca 2+ signal using increasing PCDs decreases the magnitude of change of the average predicted ERK response between cells closest to the wound and farthest from the wound ( Figure 3D). In other words, the response gradient becomes less obvious when cells are averaged over larger distances. Therefore, although the increase in PCD decreases response noise, the corresponding decrease in signal demonstrates the limit of PCD on the SNR benefit ( Figure 3E). The difference in rates at which the signal and noise decrease results in a maximum SNR at a PCD of 91.0+/-6.3 mm (SEM, N = 5) ( Figure 3F). This peak corresponds to a PCD where the amount of noise is decreased to the lowest amount possible without reducing the response gradient due to 'over-averaging'. The predicted PCD with maximal benefit did not change when we expand the model to consider limitations due to diffusion ( Figure 3F green curve, Materials and methods). Similar analysis using mutual information statistics shows a similar PCD with the maximal mutual information at the distance that showed maximal SNR ( Figure 3E,F, Figure 3-figure supplement 1B). This analysis shows that the benefits from paracrine communication depend on how far a paracrine molecule travels which, in this specific case, has a maximal benefit at~100 mm, or approximately three cell diameters.

Direct measurement of Paracrine Communication Distance
We next empirically measured the PCD in our experimental system to compare to the PCD predicted to maximize the SNR in the wound context. To measure the PCD of ERK activation we first established a co-culture system that allows us to separate the effects of autocrine and paracrine signaling. Our assay utilizes a synthetic GPCR: Designer Receptors Exclusively Activated by Designer Drugs (DREADD). The Gq human muscarinic derived GPCR DREADD is activated by a synthetic small molecule, clozapine-N-oxide (CNO), that has no known endogenous receptors (Armbruster et al., 2007). In addition, DREADD activates the Gq pathway similar to purinergic ATP receptors (Dong et al., 2010). Using a co-culture of DREADD expressing (activated by CNO) and non-expressing cells (not activated by CNO), we can determine which cells release EGF (DREADD expressingred) and which cells accept EGF (non-expressing-gray) ( Figure 4A). Using a synthetic system allows us to directly measure the average communication distance of EGF. CNO addition selectively activates Ca 2+ response in DREADD expressing cells while the surrounding non-expressing cells show no response indicating a lack of paracrine activation of Ca 2+ response in cells ( Figure 4A). Although some systems show that Ca 2+ response can propagate from cell-to-cell through gap junctions (Ross 2012), this does not appear to be the case in MCF-10A cells as non-expressing cells showed no cytosolic Ca 2+ increase upon activation of DREADD cells. Alternatively, ERK response was found in both DREADD expressing and the surrounding non-expressing cells upon CNO addition but was inhibited in both cell types in the presence of the EGFR inhibitor tryphosphitin AG1478, confirming paracrine activation of ERK in the DREADD system ( Figure 4B). The Ca 2+ and ERK responses in the DREADD system suggest that local averaging takes place only at the EGF level between Ca 2+ and ERK response. Additionally, increasing the ratio of DREADD cells to non-expressing cells shows an increasing ERK response magnitude in non-expressing cells, further supporting that paracrine communication locally activates ERK ( Figure 4C). We measured the PCD by monitoring local paracrine ERK activation with our DREADD co-culture assay. We co-cultured DREADD cells at a low concentration compared to non-expressing cells to ensure that neighboring non-expressing cells were activated only by a single DREADD cell. We then analyzed the ERK response of~1500 non-expressing cells neighboring a DREADD cell. (Figure 4D, E). Local ERK activation of non-expressing cells surrounding a DREADD cell show decreasing response with increasing distance from the DREADD cell. ERK response as a function of distance follows a Gaussian fit, consistent with how the concentration of diffusing molecules, like in paracrine signaling, changes over distance ( Figure 4E) (Berg, 1993). The PCD was determined by calculating the spread, or sigma, of this Gaussian curve. According to our fit, the paracrine activation of ERK has

t-test). (C) The average maximum ERK response of non-expressing cells in a given well increases linearly with an increasing percentage of DREADD cells per well. (D) Representative images from a timelapse experiment showing ERK activation (cyan) in non-expressing cells surrounding a single activated DREADD cell (yellow) over a 40 min time period. ERK activation level indicated by black to cyan colorbar.
(E) Average maximum ERK activation in non-expressing cells surrounding single DREADD cell according to distance from the DREADD cell. PCD was calcuated as the spread, or sigma, of the fitted Gaussian curve (dashed line) and measured to be 99.5 mm+/-19.6 mm (blue *) (SEM indicated by error bars, N = 12 DREADD cells). Scale bar represents average length of a single-cell. F. Comparison between calculated optimal PCD per wound (pink, Figure 3F) and experimentally measured PCD found using the DREADD co-culture assay per DREADD cell (blue, Figure 4E) (p-value>0.5, t-test). Horizontal bars represent average. DOI: 10.7554/eLife.09652.019 a communication distance of 99.5+/-19.6 mm (SEM, N = 12 DREADD cells). This empirically measured value is statistically similar to the predicted communication distance value that maximizes the SNR of wound response ( Figure 4F). In other words, the cellular communication distance is tuned to maximize the overall response fidelity during wound response signaling.

Discussion
Multicellular organisms utilize cellular diversity for specialization and division of labor. However, the variability between cells can be detrimental due to the potential loss of response fidelity (Uda et al., 2013;Hansen and O'Shea, 2015;Voliotis et al., 2014;Cheong et al., 2011;Selimkhanov et al., 2014). Paracrine communication can serve to share information between cells to regulate cellular variability. In this study we analyzed the benefits and limitations of paracrine communication based information sharing between cells as a mechanism to control cellular response variability.
We analyzed the limits of paracrine communication on cellular response fidelity in two cases. First we analyzed the response to a spatially uniform ligand. Our analysis reveals that, under these conditions, the magnitude of single-cell response fidelity increases as a function of the PCD with no upper bound. In order for cells to facilitate larger PCDs, cells would need to synthesize larger amounts of paracrine signaling molecules or utilize fast diffusing paracrine ligands like H 2 O 2 (Enyedi and Niethammer, 2015). The increased energy required to synthesize the additional molecules is likely to be minor in comparison to the overall energetic demand of a cell. Therefore cells could potentially take advantage of large PCDs to substantially mitigate biochemical noise. However, a spatially uniform input, while common in cell culture experiments, is likely an inadequate representation of physiological conditions. In the second case we analyzed the response of cells to spatially defined inputs in the form of a mechanical epithelial wound. We analyzed the cellular response to extracellular ATP gradients, a damage associated molecule, following a controlled wounding of an epithelial monolayer in vitro. Similar to developmental systems, an extracellular input ligand conveys positional information in wound response (Dubuis et al., 2013;Sonnemann and Bement, 2011). Analysis of the paracrine communication benefit in our novel quantitative wound response assay with defined spatial perturbations demonstrates that paracrine communication increased cellular response fidelity, but with limitations. Unlike the spatially uniform ligand in the first case, the magnitude of the response fidelity benefit varies with increasing PCD. The maximal increase of cellular response fidelity occurrs at a PCD of~100 mm, or approximately 3 cell diameters. In vivo work measuring ERK propagation using the same EKAREV FRET sensor also showed propagation extending~100 mm (Hiratsuka et al., 2015). Our results demonstrate that the paracrine information sharing benefit depends on the input ligand spatial scale, or PCD. Furthermore, empirical measurements of paracrine communication match the physiologically relevant spatial wound response maximum communication distance.
The process of wound healing is a complex multi-stage program that coordinates the action of multiple cell types over multiple timescales, from minutes to weeks, to address an acute need. The initial steps of wound healing programs propagate information concerning the wound in a manner that is appropriate to the magnitude of damage. Both inflammatory and fibrotic processes, critical steps in wound response signaling, are damaging when they go awry. Therefore, the initial cellular responses and the establishment of signaling gradients are key steps in wound healing. The mechanistic details underlying how tissues robustly match the wound response magnitude to the extent of wound-induced damage remain unknown. Our results demonstrate that intercellular communication during the initial wound response is optimized to increase overall response fidelity and provides the initial evidence that matching the wound response to wound damage is a critical aspect to wound healing programs. Future work is needed to further investigate tissue level response fidelity during wound healing programs.
Paracrine communication increases the fidelity of response at the single-cell level by mitigating biological noise at the single-cell level. Each cell integrates information from its local neighborhood to increase its individual response fidelity. Local averaging at the cellular level is a distinct mechanism compared to the benefit of global averaging at the population level. Without any paracrine communication, the reliability of the average of a cell population response can only increase with the size of the population. This is a consequence of the central limit theorem where the uncertainty of a sample average decreases with sample size. However, this increase in reliability is only true for the population average and not for individual cells in the population. Therefore, in cases where the biologically significant output is the collective action of the population, for example the secretion of a cytokine, intercellular information sharing is not required. However, when biologically significant output requires single-cell action, local information sharing via paracrine communication increases cellular response fidelity. Therefore, whether paracrine communication is required remains context dependent. It is possible that paracrine information sharing is more prevalent in signaling networks that support individual cellular decisions and less prevalent in cases where biologically meaningful outcomes result from population averages.
In cases where paracrine information sharing is used as a method to mitigate biological noise, the breakdown of this system could be detrimental. In vivo studies of ERK response in mammary tumor cells using the same EKAREV FRET sensor utilized here show highly variable ERK response that may lead to the survival and propagation of cancer cells (Kumagai et al., 2014). Although the cause for this heterogeneity is unknown, one possible mechanism may be the breakdown of paracrine communication between cells similar to how our partial inhibition of paracrine communication showed no decrease in ERK variability (Figure 1-figure supplement 4).
The abundance of paracrine communication in mammals, that is the activation of a receptor by a ligand synthesized by another cell, demonstrates the heavy utilization of intercellular communication (Ben-Shlomo et al., 2003;. Paracrine averaging demonstrates how intercellular communication enables cellular collective decision making where the 'wisdom of the crowd' is greater than the individual cell. Theoretical and empirical work in humans and animal collectives has shown that the benefit of collective decision making depends on the size of the group; big crowds are not always better than small crowds (Sasaki et al., 2013;Kao et al., 2014;Hoare et al., 2004;Sueur et al., 2011). Therefore, it is likely that the extent of secretion of each paracrine ligand is adjusted to the level of cellular information sharing to ensure an effective collective decision.
The optimal PCD we identified is not universal. Rather, the optimal distance depends on the specific shape of the spatial pattern of the initial activating ligand and the noise pattern of the primary response. Additionally, propagation patterns of the same activating ligand can depend on the physiological signaling context as demonstrated by differences found during in vivo ERK propagations under wound and normal conditions (Hiratsuka et al., 2015). The effective PCD can be regulated at the cellular level by several possible factors to optimize the benefit of paracrine communication to the specific noise and spatial patterns characteristic to each signaling system (Batsilas et al., 2003;Muratov and Shvartsman, 2003a;2003b). PCD also depends on the effective diffusion coefficient of the secreted molecule, transmitted signal strength (e.g. number of secreted molecules), and receiver cell sensitivity (e.g. receptor K d ). The diffusion coefficients of paracrine signaling molecules can vary by two orders of magnitude, let alone differences in signal strength and receptor sensitivity in individual paracrine signaling pathways (Kreuz et al., 1965;Gregor et al., 2007). Fine-tuning each of these factors provides a possible mechanism for cells to regulate the PCD and thereby the extent a cell locally communicates. The ability to specifically tune PCD raises the possibility that evolutionary pressures can tune paracrine communication to provide the optimal benefit in many other paracrine communication systems.

Materials and methods
Ca 2+ and ERK measurements in MCF-10A cells MCF-10A cells were cultured following established protocols (Debnath et al., 2003). Before plating cells, each surface was first treated with a collagen (Life Technologies, Carlsbad, CA), BSA (New England Biolabs), and fibronectin (Sigma-Aldrich) solution in order for cells to completely adhere, according to established methods. In order to maintain a viable environment, cells were imaged at 32˚C and 5% CO 2 . All EGF (PeproTech) titrations and DREADD experiments were conducted in 96well plates using extracellular hepes buffer (ECB) to reduce background fluorescence (5 mM KCl, 125 mM NaCl, 20 mM Hepes, 1.5 mM MgCl 2 , and 1.5 mM CaCl 2 , pH 7.4). All imaging for wounding was done in MCF-10A assay media (Debnath et al., 2003).

ERK and Ca 2+ activation by DREADD
Cells were plated at a density of 2,000,000 cells/100mm plate and allowed to adhere overnight. Cells were transfected with the Gq-coupled DREADD HA-tagged hM3D with an mCherry tag using a 3:1 ratio of FuGene HD (Promega) to DNA and allowed to incubate overnight (Dong et al., 2010). In order to measure the paracrine signal from a single-cell, non-transfected cells were mixed with DREADD-transfected cells at ratios of 1:0, 1:1, 1:2, 1:5, 1:7, and 0:1 (non-transfected:DREADD) and plated in 96-well plates at a density of 30,000 cells/well. The following day, cells were loaded with 1 mM Hoechst dye for nuclear imaging for 30 min for cell segmentation purposes. 5 mM clozapine-Noxide (CNO) (Enzo Life Sciences) was added to each well to specifically activate DREADD cells. ERK activation was monitored using the EKAREV FRET reporter (Albeck et al., 2013;Komatsu et al., 2011) and Ca 2+ activation was monitored using the Ca 2+ indicator dye Fluo-4 using the published protocol (Invitrogen).

Cell clustering assay and analysis
In order to measure the standard deviation of Ca 2+ and ERK activity within a small group of cells, MCF-10A cells were plated at densities of 1000, 2000, and 3000 cells per well in a 96-well plate, taking advantage of the natural tendency for MCF-10A cells to cluster together. Cells were stimulated with 10 mM or 100 mM ATP and imaged for 5 min every 3 s (Ca 2+ ) or 30 min every minute (ERK).
Standard deviation and average expression of Ca 2+ and ERK were analyzed by grouping cells in to clusters based on the distances between cells and clusters ( Figure 1-figure supplement 1). Following the cluster analysis, the average and standard deviation of Ca 2+ and ERK activation were calculated for each cluster. ERK activation was measured using the ERK FRET reporter (Albeck et al., 2013;Komatsu et al., 2011) and Ca 2+ activation was monitored using the genetically encoded sensor R-GECO (Zhao et al., 2011).

Wounding device design, fabrication, and wounding assay
Master molds for the microfluidics based wounding device were created using silicon wafers and layer-by-layer photolithography using established methods (Ferry et al., 2011). A separate mold for both the air layer and cell layer were made using negative photoresists and masks. Chips were made by pouring uncured polydimethylsiloxane (PDMS) onto each mold, allowing the PDMS to harden, and bonding the layers together and subsequently to a glass slide. Cells were loaded into the devices through the inlet port using a 20G needle. During wounding the outlet port was plugged using tape and the inlet port held a reservoir of media to prevent evaporation in the chamber. Wounding was accomplished by increasing the air pressure in the top layer of the device until the pillar made contact with the bottom of the device after which the air pressure was released to raise the pillar back up. Cells were loaded in to the wounding device at a density of 15,000,000 cells/mL using a 20G needle. Following trypsinization and resuspension, cells were put on ice to prevent aggregation. Two o-rings were attached to the device surrounding both the inlet and outlet ports for media reservoirs. Each o-ring was attached using a thin film of vacuum grease. Wounding devices were kept in an empty pipet box filled with water to prevent media evaporation. Cells were allowed to adhere for 18-24 hr before wounding.

Imaging and image analysis
Imaging was accomplished using a Nikon Plan Apo l 10X/0.45objective with a 0.7x demagnifier and Nikon Eclipse Ti microscope with a sCMOS Zyla camera. All imaging was accomplished using custom automated software written using MATLAB and Micro-Manager (Edelstein et al., 2010). Image analysis was accomplished using a custom MATLAB code published previously (Selimkhanov et al., 2014) and is available through GitHub repository https://github.com/rwollman/CellSegmentation. git.

Model for paracrine communication based on local isotropic diffusion
The paracrine ligand concentration (P) for a cell at position (x,y) observed by (P[x,y]) is the local average of the concentration of ligand released by cells in the local neighborhood ( Figure 1A). We modeled this paracrine ligand local average using a convolution of two functions: S(x,y) that represents the amount of ligand secreted by each cell and D(x,y) that represents the expected diffusion of the paracrine ligand during the timescale of paracrine signal integration: The function S(x,y) was estimated using experimental cellular Ca 2+ response data according to: The function D(x,y) was approximated to follow Gaussian weights with a length-scale we named the Paracrine Communication Distance (PCD): Detailed analysis of how PCD depends on diffusion, number of secreted paracrine molecules, sensitivity of detection, physiological levels of fluid flow (Polacheck et al., 2011), and cellular decoding of time varying paracrine signal are presented in Materials and methods. In cases where biologically relevant integrations times may influence the predicted paracrine communication response, the PCD did not exceed the approximate distance EGF could travel before the first ERK response ( Figure 1G, Figure 3F, Materials and methods). Based on single-cell ERK data to ATP stimulation, this time was found to be~5 min which resulted in a maximum PCD of~300 mm (data not shown) based on EGF diffusion coefficient of 50 mm 2 /s.

Signal to noise analysis
Signal to Noise ratio analysis on Ca 2+ response to ATP titration data was estimated as was done previously (Selimkhanov, 2014). Briefely, the signal S was calculated using: The noise N was calculated by: Where Ca 2+ (t) is the temporal time series of Ca 2+ response measured experimentally. Cells are separated into bins according to either different dosages of ATP added to multiple wells (Figure 1) or different distances from the wound source ( Figure 3). SNR was then simply: SNR = S/N.

Analysis of the effects of diffusion, secretion, and integration time on paracrine communication
In this section we analyze how the Paracrine Communication Distance (PCD), the characteristic length-scale of paracrine communication, depends on factors related to the paracrine signal. Specifically we look into how the PCD depends on the diffusion coefficient D, the number of molecules released from a cell N r , the number of molecules needed for detection N d , and the total integration time T.
To understand how PCD depends on the factors mentioned above, we considered the diffusion of a paracrine ligand from a single cell to its surrounding neighbors. We considered a 2D-like geometry where cylindrical cells, each of height h c and radius , grow in a chamber of total h f height. We simplify the below analysis by approximating the cell monolayer geometry to a series of 'cell cylinders'. The key results of the scaling of PCD and required integration time are similar for other comparable geometries (data not shown). Under these conditions one could write the analytical solution of the diffusion equations: Where Cðr; tÞ is the concentration of paracrine ligand for distance r and time t. For a neighboring cell to respond to this paracrine signal, a critical number of molecules N d needs to reach the volume surrounding the cell. We assume that a cell 'senses' a volume comparable to the volume of a cell itself. For a cylindrical cell of area p 2 and height h c the critical concentration required for cellular response will be: This is simply the required number of molecules divided by the cell volume. Combining equations 6,7 we can solve for the distance and time of where the critical concentration will be reached. Solving for distance we get that The concentration of the paracrine ligand is diluted as it diffuses from the source. Therefore, there is a point in space which is the maximal distance from the source that the critical detection concentration C detect will be reached at some point in time. Distances that are greater than the critical distance will only experience concentrations lower than the critical detection concentration C detect . The existence of such maximum can also be seen by the non-monotonous dependency of r detect on t in equation 8. To find the maximal distance we can simply find the maximum of 1.3 in respect to t. Doing so we get that: We can simplify the analysis by the introduction of two dimensionless variables: 1) S ¼ Nr Nd represents the strength of the signal and is defined as the ratio of released molecules N r and the number of molecules needed to detect the signal N d . 2) h ¼ hc hf represents the fraction of the height of the flow chamber that cells occupy. When we substitute the new variables into Equation 1,4 we get that: Interestingly this shows that the value of PCD does not depend on the diffusion coefficient. Rather, PCD scales as a function of the square root of the strength of the signal S with a multiplicative constant that depends on the specific cell geometry. PCD also depends on cell geometry with the cell radius and the relative height of a cell in the effective environment h. While the analysis above shows that the diffusion coefficient has no influence on the overall PCD, the time required to reach this maximal distance has important biological implications. 'Paracrine averaging' requires cells to integrate the signal. However, the time required for signal integration must be biologically feasible given the cellular response time and diffusion coefficient.
From equation 8 we can identify the time by which the PCD is maximal to be: The integration time grows linearly with signal strength S. This is because the PCD itself scales as a square root of S and the diffusion time grows with the square of the distance. The integration time decreases with increasing diffusion coefficient as expected. Figure 1-figure supplement 7 shows the scaling of the integration time with the diffusion coefficient for a few PCD values.
For diffusion coefficient values of~10-100 mm 2 /s and a PCD of 100 mm (similar to the distance measured in Figure 3) integration times ranged between 0.5 to 5 min. Given that ERK activation is observed only after 5 min post-activation, the required integration time does not pose an issue. However, larger PCDs will require higher diffusion coefficients to allow proper integration of the paracrine signal. Interestingly, H 2 O 2 , another key paracrine signaling molecule critical to initial wound response signaling, has a diffusion coefficient of~2000 mm 2 /s. A larger diffusion coefficient could allow for a much longer PCD with reasonable biological integration times.

Analysis of the effect of fluid flow on paracrine communication
All the analysis above assumed static conditions, that is no fluid mixing or advection of any kind. In this section we analyze the degree to which the principles of paracrine communication are applicable in non-static conditions. Non-static fluid conditions potentially have two effects on mass transport. 1) Non-static fluid conditions can create mixing due to turbulence and 2) Laminar advection can transport secreted molecules away from the secretion source. Since the extracellular environment is characterized by a low Reynold's number there is effectively no turbulent mixing in biologically relevant parameters.
To analyze the relative contribution of advection and diffusive transport we utilize a dimensionless number, the Pé clet number (P), that represents the ratio between the contribution of advection and diffusion: Where v is the interstitial flow rate, D is the diffusion coefficient and L is the characteristic length scale. In our case, the characteristic length scale is the PCD, which depends on the signal strength as described above (equation 10 and Figure 1-figure supplement 6). Therefore, the P number can be expressed as a function of the signal strength S and diffusion coefficient: Graphical representation of this expression is shown in Figure 1-figure supplement 8 where the map of D and S is color coded by the Pé clet number with three highlighted regions: A red region where flow will dominate, a cyan region where diffusion will dominate, and the region in between where both advection and diffusion contribute to paracrine communication.
To gain further insight into the relative contribution of advection and diffusion we looked at the distance molecules will travel via advection for a specific signal strength (S = 1000). As can be seen in Figure 1-figure supplement 8b, for diffusion coefficients of small protein ligands advection will contribute minimally.
When considering positional accuracy of cellular response, an important consideration is that advection can potentially 'shift' the effects of paracrine signaling downstream of the flow. Even if the shift is characterized by low Pé clet number, advection can interfere with positional information accuracy (as analyzed in Figure 3). To estimate the potentially degrading effects of flow we calculate the expected level of positional accuracy error induced by flow. We estimate that the wound induced signaling gradient ( Figure 3C) to be >500 mm. Therefore the effect on positional accuracy will be minimal (<10%) at advection distances up to 50 mm, or for a PCD of 100 mm, a Pé clet number up to 0.5. The isocline of a Pé clet number of 0.5 is shown in Figure 1-figure supplement 8a as a dotted black line. This shows that for paracrine ligands with a diffusion coefficient >40 mm 2 /sec, advection will have little effect on positional accuracy of initial wound response signaling.

Analysis of the effect of cellular decoding schemes on paracrine communication
The analysis in the previous two sections assumes that the concentration of the paracrine ligand decreases over increasing distance from the source of secretion according to a Gaussian fit where the diffusion length-scale represents the PCD. The cellular response to a paracrine ligand depends on cellular decoding of the temporal paracrine concentration profile a cell observes. As both the temporal profile of the secreted paracrine molecule and the temporal cellular decoding are unknown, we consider the simple assumption of a Gaussian profile reasonable. To quantitatively test this assumption we compared the Gaussian profile to an alternative model that could be addressed analytically. In the alternative model, we assume that all paracrine molecules are released at T = 0 and that cellular decoding of the paracrine signal is simple temporal averaging. Under these assumptions one can write an expression of the temporal average of the paracrine concentration at a distance r from the source as: Where all symbols follow equation 6 and ei represent the exponential integral: eiðxÞ ¼ ð x À¥ e t t dt Comparison of the two models can be seen in Figure 1-figure supplement 9. Overall, the two models generate very similar Paracrine Averaging Weights (the effect of cellular decoding of paracrine signal). There is a small discrepancy between the two models at very low distances (<50 mm). However, this discrepancy is most likely a result of the assumption in the alternative model that all the paracrine molecules are released at once. Under the more realistic assumption where paracrine molecule release duration is not much smaller than the time to diffuse 50 mm (12.5 s at D = 50 mm 2 / s) we anticipate that the similarity between these two profiles will further increase.