The ability to sense the environment is heterogeneously distributed in cell populations
Abstract
Channel capacity of signaling networks quantifies their fidelity in sensing extracellular inputs. Low estimates of channel capacities for several mammalian signaling networks suggest that cells can barely detect the presence/absence of environmental signals. However, given the extensive heterogeneity and temporal stability of cell state variables, we hypothesize that the sensing ability itself may depend on the state of the cells. In this work, we present an informationtheoretic framework to quantify the distribution of sensing abilities from singlecell data. Using data on two mammalian pathways, we show that sensing abilities are widely distributed in the population and most cells achieve better resolution of inputs compared to an ‘average cell’. We verify these predictions using livecell imaging data on the IGFR/FoxO pathway. Importantly, we identify cell state variables that correlate with cells’ sensing abilities. This informationtheoretic framework will significantly improve our understanding of how cells sense in their environment.
eLife assessment
In this valuable article, the authors use an existing theoretical framework relying on information theory and maximum entropy inference in order to quantify how much information single cells can carry, taking into account their internal state. They reanalyze experimental data in this light. Despite some limitations of the data, the study convincingly highlights the difference between singlecell and population channel capacities. This result should be of interest to the quantitative biology community as it contributes to explaining why channel capacities are apparently low in cells.
https://doi.org/10.7554/eLife.87747.3.sa0Introduction
In cell populations, there is a significant overlap in responses to environmental stimuli of differing strengths. This raises a fundamental question (Levchenko and Nemenman, 2014): do signaling networks in cells relay accurate information about their environment to take appropriate action? And if not, where along the signal transduction pathway is the information lost (Rhee et al., 2012; ten Wolde et al., 2016)? Mutual information (MI) quantifies the information content in an intracellular output about extracellular inputs. For an input $u$ (e.g., concentration of a ligand) and an output $x$ (e.g., intracellular species concentration, Cheong et al., 2011; Selimkhanov et al., 2014; Suderman et al., 2017) or a cellular phenotype (Varennes et al., 2019; Mattingly et al., 2021), the MI is defined as (Cover and Thomas, 2006)
Experimental singlecell methods such as flow cytometry (Wu and Singh, 2012), immunofluorescence (Wu and Singh, 2012), mass spectrometry (Slavov, 2021), or livecell imaging (Lemon and McDole, 2020) allow us to estimate response histograms $p\left(xu\right)$ across several inputs. Using these distributions, we can estimate the maximum of the MI (the channel capacity, CC) by optimizing Equation 1 over input distributions $p\left(u\right)$. The CC quantifies fidelity of signal transduction (Levchenko and Nemenman, 2014; Rhee et al., 2012). For example, a CC of 1 bit implies that that the cells can at best distinguish between two input levels (e.g., presence versus absence), with higher CCs indicating that cells can resolve multiple input states. Importantly, CC can be used to identify bottlenecks in signaling (Rhee et al., 2012; Hansen and O’Shea, 2015).
The CC has been estimated for input–output relationships in several mammalian signaling networks (Cheong et al., 2011; Selimkhanov et al., 2014; Suderman et al., 2017; Billing et al., 2019; Garner et al., 2016; Frick et al., 2017). When the output is defined as levels of a single protein a fixed time, the CC was found to be surprisingly low, ∼1–1.5 bits. These estimates have been improved by considering multidimensional outputs (Selimkhanov et al., 2014) or timevarying inputs (Lee et al., 2021). While these modifications led to somewhat higher CC estimates, the overall conclusion that cells know little about their environment remains well established. In contrast, significantly higher CC estimates have been found when the output at the level of cell population averages is considered (Suderman et al., 2017), suggesting that the only way to overcome lowsensing fidelity is population average response.
These previous calculations estimated onechannel capacity for all cells in a population, implicitly assuming that individual cells have similar sensing abilities. However, we now know that cells exhibit extensive heterogeneity in cell state variables (Trapnell, 2015; Wollman, 2018) such as abundances of key signaling proteins (Meyer et al., 2012; Frankel et al., 2014), mRNA abundances (Chen et al., 2015), and chromatin conformation (Nagano et al., 2013) and accessibility (Buenrostro et al., 2015). Notably, the time scale of fluctuations in these variables can be significantly slower than relevant signaling time scales (Flusberg et al., 2013; Spencer et al., 2009), sometimes extending across multiple generations (Pleška et al., 2021). Heterogeneity in cell state variables leads to a heterogeneity in response to extracellular cues, including chemotherapeutic drugs (Flusberg et al., 2013; Spencer et al., 2009; Emert et al., 2021), mitogens (Matson and Cook, 2017), hormones (Norris et al., 2021), chemotactic signals (Varennes et al., 2019; Pleška et al., 2021; Moon et al., 2021; Waite et al., 2016), and other electrical and chemical stimuli (Clark and Campbell, 2019; Yao et al., 2016). Therefore, we hypothesize that that the ability to sense the environment varies from cell to cell in populations in a cell statedependent manner.
There is no conceptual framework to estimate the variation in sensing abilities in cell populations and its dependence on cell state variables. To that end, we introduce an informationtheoretic quantity CeeMI: Cell statedependent Mutual Information. We show that using typically collected singlecell data and computational models of signaling networks, we can estimate the distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ of singlecell signaling fidelities (singlecell mutual information values). We also show that we can identify cell state variables that make some cells better and others worse at sensing their environment.
Using an illustrative example, we show that in heterogeneous cell populations, estimates of mutual information that average over cell states can be significantly lower than the mutual information of signaling networks in typical cells. Next, using previously collected experimental data, we estimate ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ for two important mammalian signaling pathways (Lewis Cantley and Sever, 2013): the epidermal growth factor (EGF)/EGF receptor pathway and the insulinlike growth factor (IGF)/Forkhead Box protein O (FoxO) pathway. We show that while the cell stateagnostic CC estimates for both pathways are ∼1 bit, most individual cells are predicted to be significantly better at resolving different inputs. Using livecell imaging data for the IGF/FoxO pathway, we show that our estimate of variability in sensing abilities matches closely with experimental estimates. Importantly, for this pathway, we also verify our prediction that specific cell state variables dictate cells’ sensing abilities. Finally, using a simple receptor/ligand model, we show how the time scales cell state dynamics affects cells’ individuality and therefore sensing abilities. We believe that CeeMI will be an important tool to explore variability in cellular sensing abilities and in identifying what makes some cells better and others worse in sensing their environment.
Results
Conditional mutual information models single cells as cell statedependent channels
Consider a cell population where cells are characterized by state variables $\mathit{\theta}$ (Figure 1). These include abundances of signaling proteins and enzymes, epigenetic states, etc. We assume that cell states are temporally stable, that is, $\mathit{\theta}$ remains constant over a time scale that is longer than typical fluctuations in environmental inputs. Later, using a toy simple model, we explicitly study the effect of cell state dynamics on cells’ ability to sense their environment.
We assume that cell state variables are distributed in the population according to a distribution $p\left(\mathit{\theta}\right)$ . If $x$ denotes an output of choice (e.g., phosphorylation levels of one or more protein(s) at one or more time point(s)) and $u$ denotes the input (e.g., ligand concentration), the experimentally measured response distribution $p\left(xu\right)$ can be decomposed as
where $p\left(xu,\mathit{\theta}\right)$ is the distribution of the output $x$ conditioned on the input $u$ and cell state variables $\mathit{\theta}$. We note that in most cases the same cell cannot be probed multiple times. Consequently, $p\left(xu,\mathit{\theta}\right)$ may not be experimentally accessible. However, it is conceptually well defined and mathematically accessible when interactions among molecules are specified (Feinberg and Levchenko, 2023).
Using $p\left(xu,\mathit{\theta}\right)$ , we can define the cell statedependent mutual information $I\left(\mathit{\theta}\right)$ for a fixed input distribution $p\left(u\right)$ analogously to Equation 1:
$I\left(\mathit{\theta}\right)$ quantifies individual cells’ ability to sense their environment as a function of the cell state parameters $\mathit{\theta}$. The distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ of singlecell sensing abilities is
where $\delta \left(\cdot \right)$ is the Dirac delta function. We can also compute the joint distribution between the singlecell mutual information and any cell state variable of interest $\chi $ (e.g., abundance of cell surface receptors):
where $\chi \left(\mathit{\theta}\right)$ is the value of the biochemical parameter when cell state variables are fixed at $\mathit{\theta}$. The distributions in Equation 5 quantify the interdependency between a cell’s signaling fidelity $I\left(\mathit{\theta}\right)$ and cellspecific biochemical parameters $\chi \left(\mathit{\theta}\right)$. As we will see below, the distributions in Equations 4 and 5 can be experimentally verified when appropriate measurements are available. Finally, we define the population average of the cell statedependent mutual information:
In information theory, ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ is known as the conditional mutual information (Cover and Thomas, 2006) between the input $u$ and the output $x$ conditioned on $\mathit{\theta}$. If cell state variables remain constant and if distribution over cell states is independent of the input distribution, it can be shown that ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}\ge I$ (Section 1 of Materials and Methods). It then follows that at least some cells in the population have better sensing abilities compared to the cell stateagnostic mutual information (Equation 1).
Since ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ depends on the input distribution $p\left(u\right),$ we can find an optimal input distribution $p\left(u\right)$ that maximizes ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ (Section 1 of Materials and Methods). Going forward, unless an input distribution is specified, the distributions ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ and ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I,\chi \right)$ are discussed in the context of this optimal input distribution.
Maximum entropy inference can estimate ${\mathit{p}}_{\mathbf{C}\mathbf{e}\mathbf{e}\mathbf{M}\mathbf{I}}\left(\mathit{I}\right)$
Estimation of ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ requires decomposing the experimentally observed response $p\left(xu\right)$ into cellspecific output distributions $p\left(xu,\mathit{\theta}\right)$ and the distribution of cell state variables $p\left(\mathit{\theta}\right)$ (Equations 3 and 4). This problem is difficult to solve given that neither $p\left(xu,\mathit{\theta}\right)$ nor $p\left(\mathit{\theta}\right)$ are accessible in typical experiments. However, for many signaling networks, stochastic biochemical models can be constructed to obtain the cellspecific output distribution $p\left(xu,\mathit{\theta}\right)$ . Here, $\mathit{\theta}$ comprise protein concentrations and effective rates of enzymatic reactions and serve as a proxy for cell state variables. Given the experimentally measured cell stateaveraged response $p\left(xu\right)$ and the modelpredicted cellspecific output distribution $p\left(xu,\mathit{\theta}\right)$ , we need a computational method to infer $p\left(\mathit{\theta}\right)$ (see Equation 2). The problem of inference of parameter heterogeneity is a nontrivial inverse problem (Dixit et al., 2020). While there are several proposed methods to solve this problem (reviewed in Loos and Hasenauer, 2019), most cannot efficiently infer $p\left(\mathit{\theta}\right)$ for signaling networks with even a moderate ($n\sim 10$) number of parameters. Here, we use our previously develop maximum entropybased approach to infer $p\left(\mathit{\theta}\right)$ (Dixit et al., 2020). This way, we can estimate ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ using experimentally obtained cell stateagnostic response $p\left(xu\right)$ and a stochastic biochemical model $p\left(xu,\mathit{\theta}\right)$ of the underlying signaling network.
An ‘average cell’ can discern much less than a typical cell about the environment
To illustrate the effect of heterogeneity of cell state variables on the cell stateagnostic estimate of mutual information (which we call $I}_{\mathrm{C}\mathrm{S}\mathrm{A}$ from now onward), we consider a simple stochastic biochemical network of a receptorligand system. Details of the model and the calculations presented below can be found in Section 2 of Materials and Methods. Briefly, extracellular ligand $L$ binds to cell surface receptors $R.$ Steadystate levels of the ligandbound receptor $B$ are considered the output. The signaling network is parameterized by several cell state variables $\mathit{\theta}$ such as receptor levels, rates of binding/unbinding, etc. For simplicity, we assume that only one variable, steadystate receptor level ${R}_{0}$ in the absence of the ligand, varies in the population. Calculations with variability in other parameters are presented in Section 2 of Materials and Methods.
In this population, cells’ response $p\left(BL,{R}_{0}\right)$ is distributed as a Poisson distribution whose mean is proportional to the cell state variable ${R}_{0}$ (Section 2 of Materials and Methods). That is, when all other parameters are fixed, a higher ${R}_{0}$ leads to lower noise (coefficient of variation). To calculate cell statedependent mutual information (Equation 3), we assume that $p\left(L\right)$ is a gamma distribution. As expected, $I\left({R}_{0}\right)$ (Equation 3) between the output $B$ and the input $L$ increases monotonically with ${R}_{0}$ (inset in Figure 2A). Moreover, given that ${R}_{0}$ varies in the population (also assumed to be a gamma distribution), the sensing ability varies in the population as well (Figure 2A). Notably, the average ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ of $I\left({R}_{0}\right)$ remains relatively robust to variation in ${R}_{0}$ . At the same time, the traditional estimate ${I}_{\mathrm{C}\mathrm{S}\mathrm{A}}$ , which is the mutual information between the input $L$ and the cell stateagnostic population response $p\left(BL\right)=\int p\left(BL,{R}_{0}\right)p\left({R}_{0}\right)d{R}_{0}$ (response of the ‘average’ cell, Equations 1 and 2), decreases as the population heterogeneity in ${R}_{0}$ increases. Importantly, ${I}_{\mathrm{C}\mathrm{S}\mathrm{A}}$ is significantly lower than the sensing ability of most cells in the population (Figure 2A). This is because the overlap in the population response distributions is significantly larger than that in singlecell response distributions (Figure 2B). This simple example illustrates that the traditional mutual information estimates may severely underestimate cells’ ability to resolve inputs, especially when cell state variables are heterogeneously distributed. Moreover, it is crucial to explicitly account for heterogeneity in cell state variables when estimating fidelity of cellular communication channels.
Experimental verification of ${\mathit{p}}_{\mathbf{C}\mathbf{e}\mathbf{e}\mathbf{M}\mathbf{I}}\left(\mathit{I}\right)$ using growth factor signaling networks
In cell populations, state variables that govern signaling dynamics such as protein levels (receptors, kinases, dephosphatases, etc.) (Meyer et al., 2012; Frankel et al., 2014), as well as effective rate constants such as endocytosis rates (Kallenberger et al., 2017), ligand binding rates (Chung et al., 1997), etc., vary from cell to cell. Therefore, we expect the downstream phenotype of environmental sensing to be widely distributed as well. To experimentally verify the computational prediction of the distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ of sensing abilities, we need a system that allows us to approximate the cell statespecific response distribution $p\left(xu,\mathit{\theta}\right)$ . The IGF/FoxO pathway (Lewis Cantley and Sever, 2013) is an ideal system for these explorations for several reasons. First, following IGF stimulation, the transcription factor FoxO is pulled out of the nucleus. GFPtagged variant of FoxO can be used to detect the dynamics of nuclear levels of FoxO at the singlecell level (Gross et al., 2019). Second, nuclear FoxO levels reach an approximate steady state within 30–45 min of constant IGF stimulation, with FoxO levels decreasing with increasing IGF dose (Gross et al., 2019). As a result, an approximate cell statespecific distribution $p\left(xu,\mathit{\theta}\right)$ of steadystate levels of FoxO can be obtained by stimulating the same cell with successive IGF levels. Finally, the biochemical interactions in the IGF/FoxO are well studied (Figure 2C), allowing us to build a stochastic biochemical model based on previous computational efforts (Wimmer et al., 2014; Lyashenko et al., 2020) that fits the singlecell data accurately. Another system where ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ can in principle be verified is the EGF/EGFR pathway (Figure 2C). Here too, abundance of cell surface EGFR levels can be tracked in live cells following EGF stimulation (Fortian and Sorkin, 2014), allowing us to obtain cell statespecific response distribution $p\left(xu,\mathit{\theta}\right)$ . Below, we show estimates of ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ for both pathways and an experimental verification of our estimates for the IGF/FoxO pathway where livecell imaging data was previously collected (Gross et al., 2019).
The details of the calculations presented below can be found in Materials and Methods (Figure 2—figure supplement 7; Figure 2—figure supplement 6; Figure 2—figure supplement 5; Figure 2—figure supplement 3; Figure 2—figure supplement 2; Figure 2—figure supplement 1; Figure 2—figure supplements 4). Briefly, we first constructed stochastic biochemical models for the two pathways based on known interactions (Figure 2C) and previous models (Wimmer et al., 2014; Lyashenko et al., 2020). The output for the EGFR pathway was defined as the steadystate levels of cell surface EGF receptors, and the output for the IGF/FoxO pathway was defined as the steadystate nuclear levels of the transcription factor FoxO. Using previously collected singlecell data on the two pathways (Dixit et al., 2020; Gross et al., 2019; Lyashenko et al., 2020) and our maximum entropybased framework (Dixit et al., 2020), we estimated the distribution over parameters $p\left(\mathit{\theta}\right)$ for the two pathways (Section 3 of Materials and Methods). Using these inferred distributions, and the modelpredicted cell statespecific response distribution $p\left(xu,\mathit{\theta}\right)$ , we could compute ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ for any specified input distribution $p\left(u\right)$. We choose the support of the input distribution as the ligand concentrations used in the experimental setup. The estimates shown in Figure 2D and E show ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ corresponding to the input distribution $p\left(u\right)$ that maximizes the conditional mutual information ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ (see Equation 6). This input distribution is shown in the inset of Figure 2D and E.
Similar to the illustrative example (Figure 2A), there is a wide distribution of singlecell sensing fidelities in real populations (Figure 2D and E). Moreover, most cells are better sensors compared to the ‘average cell’, a cell whose response $p\left(xu\right)$ is averaged over cell state variability which was estimated to be ∼1 bit for both pathways. The cellular signaling fidelity skews toward the upper limit of 2 bits, which corresponds to the logarithm of the number of inputs used in the experiment. Indeed, as seen in the insets of Figure 2D and E and , the input distribution corresponding to the maximum of the cell stateagnostic mutual information ${I}_{\mathrm{C}\mathrm{S}\mathrm{A}}$ is concentrated on the lowest and the highest input for both pathways, indicating that cells may be able to detect only two input levels. In contrast, the input distribution corresponding to the maximum of ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ is close to uniform, suggesting that individual cells can in fact resolve different ligand levels.
To verify our computational estimate of ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ for the IGF/FoxO pathway, we reanalyzed previously collected data wherein several cells were stimulated using successive IGF levels (Gross et al., 2019). The details of our calculations can be found in Section 4 of Materials and methods. Briefly, the cells reached an approximate steady state within 60 min of each stimulation and nuclear FoxO levels measured between 60 and 90 min were used to approximate an experimental cell statespecific response distribution $p\left(xu,\mathit{\theta}\right)$ . The distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ was then obtained by maximizing the average mutual information ${I}_{\mathrm{C}\mathrm{e}\mathrm{e}}$ (averaged over all cells) with respect to the input distribution. As seen in Figure 2D and E, the experimentally evaluated singlecell signaling fidelities match closely with our computational estimates. Moreover, as predicted using our computational analysis, individual cells in the population were significantly better at sensing than what is implied by the maximum of ${I}_{\mathrm{C}\mathrm{S}\mathrm{A}}.$ Indeed, the distribution of steadystate FoxO levels was found to be well resolved at the singlecell level as well (Figure 2E). Livecell imaging data for the EGFR pathway was not available, and we leave it to future studies to validate our predictions.
Our calculations show that real cell populations comprise cells that have differing sensing fidelities, individual cells are significantly better at sensing their environment than what traditional estimates would indicate, and importantly, the CeeMI approach can accurately estimate the variation in signaling performances using readily collected experimental data and stochastic biochemical models. Notably, the variability in cell states and therefore the heterogeneity in sensing abilities are likely to be stable over time; the same cell’s FoxO responses to the same input were found to have significantly less variation compared to the variation within the population (Figure 2—figure supplement 7).
CeeMI identifies biochemical parameters that determine cells’ ability to sense their environment
Like other phenotypes (Frankel et al., 2014) and signaling dynamics (Meyer et al., 2012; Norris et al., 2021), we expect that cells’ ability to sense their environment depends on their state variables. For example, cells with faster endocytosis rates (Kallenberger et al., 2017) may integrate environmental fluctuations with higher accuracy (Lyashenko et al., 2020). Or cells with higher receptor numbers (Frankel et al., 2014) will lead to a lowered relative noise in bound ligand concentration.
To systematically identify the variables that differentiate between cells’ ability to sense the environment, we quantify the joint distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I,\chi \right)$ of singlecell signaling fidelity and biochemical state variables that take part in the signaling network. To test whether we can identify variables that are predictive of cellular fidelity, we estimated the joint distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I,\chi \right)$ for two variables that were experimentally accessible, response range of nuclear FoxO (Figure 3, left, see inset) and total nuclear FoxO levels prior to IGF stimulation (Figure 3, right, see inset). In both figures, the shaded regions show computational estimates of the joint probability distributions, and the red circles represent real cells. The green and the cyan trend lines represent computational and experimental binned averages.
One may expect that higher total nuclear FoxO levels could result in lower noise (coefficient of variation) and therefore better sensing abilities. However, we find that total nuclear FoxO levels only weakly correlate with cell statedependent mutual information (correlation coefficient $r=0.16$ for computational estimates and $r=0.04$ for experimental data). In contrast, cell statedependent mutual information depended strongly on the dynamic range of the response (correlation coefficient $r=0.53$ for computational estimates and $r=0.29$ for experimental data). Importantly, the model captured the observation that cells with a small response range had a variable sensing abilities while cells with a large response range all performed well in resolving extracellular IGF levels. Surprisingly, the total nuclear FoxO levels only weakly correlated with the cellspecific mutual information. In Section 5 of Materials and Methods (Figure 3—figure supplement 1), we show the predicted joint distributions ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I,\chi \right)$ for several other biochemical variables that can potentially govern single cells’ response to extracellular IGF stimuli. This way, CeeMI can be used to systematically identify cell state variables that differentiate between good sensors and bad sensors of the environmental stimuli.
Time scale of stochastic dynamics of cell states dictates the divergence between statespecific and stateagnostic sensing abilities
A key limitation of the presented analysis is the assumption that cell state variables remain approximately constant over the time scale of typical environmental fluctuations. While many cell state variables change slowly, remaining roughly constant over the life spans of cells and beyond (Flusberg et al., 2013; Spencer et al., 2009), state changes may occur within the life span of a cell as well (MacLean et al., 2018). These dynamical changes can be accommodated easily in our calculations. Here, instead of fixing the cell state variables $\mathit{\theta},$ we can treat them as initial conditions and propagate them stochastically with their own dynamics. In the limit of very fast dynamics where individual cells rapidly transition cell states, we expect that the individuality of cells in a population vanishes and consequently cell statespecific mutual information for each cell will agree with traditional cell stateagnostic estimates of the channel capacity. In contrast, if the cell state dynamics are slow, our framework highlights differences between cells in the population.
To elucidate the role of cell state dynamics, we built a simple model of ligand–receptor system (Section 6 of Materials and Methods, Figure 4). Briefly, the model included production and degradation of the receptor mRNA, translation of the receptor protein from the mRNA, and ligand binding to the receptor (Figure 4A). The number of ligandfree receptors and the number mRNA molecules were together considered the state variables. We tuned the mRNA dynamics by keeping the average mRNA number constant while simultaneously changing mRNA production and degradation rates. The time scale of mRNA dynamics is denoted by $\tau .$
When the mRNA dynamics are slower than translation, ligand binding, and receptor degradation, individual cells are effectively frozen around a particular receptor number (Figure 4B). The distribution of the number of ligandbound receptors after a shortterm ligand exposure reflects this frozen state. As a result, the cell stateconditioned mutual information depends strongly the cell state, with higher mutual information for cells with higher receptor numbers. In contrast, as the mRNA dynamics become faster, cells quickly lose the memory of cell state conditioning. As a result, the distribution of the number of ligandbound receptors after a shortterm ligand exposure reflects an ergodic averaging over the underlying mRNA dynamics. This averaging effectively washes out differences between cells (Figure 4C).
This example shows that when cells states change at time scales slower than signaling dynamics, cells in a population can be identified by their cell state and their ability to sense the environment can differ from one another.
Discussion
Cell populations are characterized by heterogeneity in cell state variables (Trapnell, 2015; Wollman, 2018) that is responsible for important phenotypic differences and selective advantage, for example, in sensitivity to drugs (Flusberg et al., 2013; Spencer et al., 2009; Emert et al., 2021), response to chemotactic signals (Mattingly and Emonet, 2022), and following proliferation cues (Matson and Cook, 2017). Therefore, it is reasonable to expect that cells’ ability to sense their environment depends on their state and is therefore variable across cells in a population. To quantify this heterogeneity, here, we developed a novel informationtheoretic perspective that allowed us to quantify the distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I\right)$ of singlecell sensing abilities using easily measurable singlecell data and stochastic models of signaling networks. We also quantified the joint distribution ${p}_{\mathrm{C}\mathrm{e}\mathrm{e}\mathrm{M}\mathrm{I}}\left(I,\chi \right)$ of cellspecific sensing ability and biochemical cell state variable. Importantly, using two growth factor pathways, we showed that individual cells in real cell populations were much better at sensing their environment compared to what is implied by the traditional estimate of channel capacities of signaling networks. Typical singlecell data are time stamped and do not give information about dynamics in the single cell. Moreover, these data do not provide information about dependence of celltocell variability on cell state variables. Therefore, computational modeling of cellular trajectories from singlecell data and stochastic models is essential in deciphering mechanistic origins of heterogeneity in cellular phenotypes.
The approach presented here will be useful in identifying bottlenecks in signal transduction. Many cellular phenotypes such as chemotaxis (Varennes et al., 2019; Mattingly et al., 2021) and cell proliferation (Gross and Rotwein, 2016) exhibit a weak correlation between cellular outputs (e.g., directional alignment with chemical gradients) with the input (e.g., gradient strength), resulting in a low channel capacity even for individual cells. In such cases, it is important to understand where exactly along the information transduction pathway is the information about the gradient is lost. If traditional calculations are pursued, for example, for movement of mammalian cells under growth factor gradients (Varennes et al., 2019; Moon et al., 2021), one may conclude that the information loss likely happens right at the receptor level (Figure 2D). In contrast, CeeMI will allow us to disentangle the effect of cell state heterogeneity and noisy cellular response to precisely pinpoint intracellular signaling nodes that are responsible for signal corruption.
How do we contrast our results with previous low estimates of channel capacities? There is extensive population heterogeneity in cell state variables and this heterogeneity often is stable over time. Nonetheless, this heterogeneity arises from stochasticity in underlying processes including dynamics of epigenetic transitions and gene expression. When one specifies state variables (Wollman, 2018) $\mathit{\theta}$, one effectively conditions on the stochastic trajectory of the underlying dynamics to a specific narrow window (Figure 4). This extra conditioning narrows the cells’ response distributions, revealing the heterogeneous nature of cell populations (Wollman, 2018).
In summary, we showed that like other phenotypes the ability to sense the environment is itself heterogeneously distributed in a cell population. Moreover, we also showed that when conditioned on cell state variables, mammalian cells appear to be significantly better at sensing their environment than what traditional mutual information calculation suggests. Finally, we showed that we could identify cell state variables that made some cells better sensors compared to others. We believe that CeeMI will be an important framework in quantifying fidelity of input/output relationships in heterogeneous cell populations.
Materials and methods
Scripts used to generate mutual information results as well the experimental data used to train our models can be obtained at https://github.com/adgoetz186/Cell_signalling_information (copy archived at Goetz, 2024). Script for maximum entropy inference of cellular parameters can be obtained at https://github.com/hodaakl/MaxEnt (copy archived at Akl, 2024).
The description of the methods is organized as follows: Section 1 describes basic information theoretical concepts and our numerical approach to estimate mutual information and channel capacity. Section 2 describes the in silico toy model. Section 3 describes the maximum entropy approach. Section 4 describes numerical methods for quantification of mutual information and channel capacity using livecell imagining data. Section 5 shows the joint distribution of mutual information and biochemical parameters. Section 6 describes the methods to examine the role of cell state dynamics on sensing ability.
1. Information theory primer
Here, we give a brief description of information theoretic quantities used in the article. The readers are referred to a textbook (Cover and Thomas, 2006) for a detailed discussion.
A communication channel (e.g., a signaling network) is an input/output relationship between two random variables: the input $U$ (say, a ligand concentration) and the response $R$ (e.g., levels of some intracellular protein). The mutual information (MI) between $U$ and $R$ is the reduction in the uncertainty about $U$ due to the access of the outcome of $R$. The MI is defined as
where the choice of two in the base of the logarithm gives the value of information in bits. The summation is replaced by integrals when considering continuous random variables. In the context of cell signaling, we will use the notation $I}_{CSA$ (instead of I in Equation S1) to denote the cell stateagnostic mutual information, which is the mutual information between the response distribution $p\left(ru\right)$ and the input distribution $p\left(u\right)$ .
We also consider a more nuanced situation where instead of a single channel we have a family of channels whose states are characterized by a random variable $\mathrm{\Theta}$. The channel state $\mathit{\theta}\in \mathrm{\Theta}$ uniquely determines the probabilistic response relationship $p\left(ru,\mathit{\theta}\right)$ between $U$ and $R$ when $\mathrm{\Theta}$ is fixed. Similar to Equation S1, we can define the MI between $U$ and $R$ conditioned for a specific realization $\mathit{\theta}\in \mathrm{\Theta}$
The average of $I\left(\mathit{\theta}\right)$ over $p\left(\mathit{\theta}\right)$ is traditionally known as the conditional mutual information (Cover and Thomas, 2006). We define the Cell statedependent mutual information $I}_{Cee$
The channel capacity (CC) is a measure of the optimal performance of the signaling network with respect to the distribution of $U$. It can be defined for both $I$ and $I}_{Cee$ as follows:
The difference between $I$ and $I}_{Cee$ is called the interaction information $I\left(U;R;\mathrm{\Theta}\right)$. The interaction information can be simplified as
where $I\left(U;\mathrm{\Theta}\right)=0$ follows from the statistical independence of the input signal $U$ and channel state $\mathrm{\Theta}$. Equation S5 shows that ${I}_{Cee}>I$ as long as $I\left(U;\mathrm{\Theta}\right)=0$. We note that this may not be true in general. However, it is true in the context of cellular signaling networks where the inputs are chosen by the experimentalists while the cell states are an inherent property of the cell population. It is trivial to extend this argument to the channel capacities of both terms.
Numerical estimation of mutual information
Request a detailed protocolEvaluating the mutual information between an input and an output (Equation S1 and S2) requires numerical integration over the input and the output distribution. We limit the input distributions to a finite support, specifically to the ligand concentrations that were used in the experimental setup. These are $L=\left[0,0.0078,0.01,0.03,0.06,0.125,0.25,0.5,1,100\right]$ ng/mL for the EGF/EGFR pathway and $L=\left[0,\text{}17.5,\text{}37.5,\text{}125\right]$ pM for the IGF/FoxO pathway.
The numerical integration in Equation S1 and S2 requires summing over probabilities of all possible responses and inputs. When the response distributions are approximated using histograms (either from experimental data or from Markov chain Monte Carlo simulations of a model), the summation can be error prone. To avoid this, we assume that the responses are distributed according to a gamma distribution (see Figure 2—figure supplement 4). The gamma distribution was chosen here and in several other places because it has been shown to accurately approximate real distributions of protein/mRNA abundances (Taniguchi et al., 2010).
To speed up the calculations, inspired by previous approaches (Rhee et al., 2012), we use a binning strategy. Specifically, we bin the response distribution using a constant bin width. The bin width was chosen to be equal to 5% of the smallest interquartile range across response distributions corresponding to all considered inputs. The binning procedure requires truncating the response distributions to a finite support. We ensured that our binning captured at least 99.95% of the entire mass of the distribution. The same strategy was used to obtain $I}_{CSA$ in Equation S1 and $I\left(\mathit{\theta}\right)$ in Equation S2. The only difference in computing $I\left(\mathit{\theta}\right)$ compared to $I}_{CSA$ was that the response distribution $p\left(ru,\mathit{\theta}\right)$ was obtained computationally using a stochastic differential equation model with network parameters fixed at $\mathit{\theta}$. Samples from the joint distribution ${p}_{CeeMI}\left(I,\chi \right)$ where $\chi$ is a cell state variable were obtained by sampling cell state variables $\mathit{\theta}$ from $p\left(\mathit{\theta}\right)$ (see Section 3) and simultaneously evaluating $I\left(\mathit{\theta}\right)$ and $\chi \left(\mathit{\theta}\right)$.
Once the mutual information could be obtained numerically for a given input distribution $p\left(u\right)$, the corresponding channel capacity, the maximum of the mutual information over all possible input distributions, can be obtained by solving the following optimization problem:
The optimization problem was solved via a trust region constrained algorithm (Byrd et al., 1999) using the SciPy optimization library (Virtanen et al., 2020). When the distribution over cell state variables $p\left(\mathit{\theta}\right)$ is available, we can estimate $I}_{Cee$ as the average $\u27e8I\left(\mathit{\theta}\right)\u27e9}_{\mathit{\theta}$ (Equation S3). The optimum value of $I}_{Cee$ over all input distributions serves as a cell statedependent analogue to the channel capacity of $I}_{CSA$ and is obtained by solve a similar optimization problem:
The procedure needed to evaluate singlecell mutual information values using livecell imaging data on the IGF/FoxO pathway is described in Section 4.
2. Description of the in silico toy model
The in silico cell receptor/ligand system comprises two components (unbound receptors and bound receptors) and five reactions (Figure 2—figure supplement 1).
In the model, receptors are constantly shuttled to the membrane and removed from the membrane and degraded. The extracellular ligand (concentration denoted by $L$) binds to cell surface receptors. We assume that the ligand concentration is kept constant in the environment. Steadystate abundance of ligandbound receptor (denoted by $B$) is taken as the output of the system. At steady state, the mean field equations for the average species level are
Solving for steady state, at the singlecell level, the mean number of bound receptors is given by
In Equation S10, ${k}_{\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{d}}$ is the ligand binding rate, ${k}_{\mathrm{u}\mathrm{n}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{d}}$ is the ligand unbinding rate, ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ is the degradation rate, and ${R}_{0}$ is the average value of the receptor level in the absence of the ligand. The bound receptor levels are Poisson distributed with a mean given by Equation S10 (Figure 2—figure supplement 2).
Using this toy network, we created two in silico cell populations. In both populations, we fixed $k}_{\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{d}}=1\phantom{\rule{thinmathspace}{0ex}}se{c}^{1}a.u{.}^{1$ and $k}_{\mathrm{u}\mathrm{n}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{d}}=10\phantom{\rule{thinmathspace}{0ex}}se{c}^{1$ . In the first population, every parameter was kept constant across cells except for the cell surface receptor levels ${R}_{0}$ . In the second population, every parameter was kept constant across cells except for the receptor degradation rate ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ . In the first population (when ${R}_{0}$ was varied), we fixed $k}_{\mathrm{d}\mathrm{e}\mathrm{g}}=5\phantom{\rule{thinmathspace}{0ex}}se{c}^{1$ . In the second population (when ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ was varied), we fixed $R}_{0}=50\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{s}/\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l$. The parameter that varied across cells was assumed to be distributed according to a gamma distribution. For the two populations, we kept fixed mean value of the variable parameter to be $\u27e8{R}_{0}\u27e9=500\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{s}/\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$ and $\u27e8{k}_{\mathrm{d}\mathrm{e}\mathrm{g}}\u27e9=5\phantom{\rule{thinmathspace}{0ex}}se{c}^{1}$ . We varied coefficients of variation for both populations between $CV={10}^{1.5}$ and $CV={10}^{0.5}.$
Obtaining ${\mathit{I}}_{\mathit{C}\mathit{S}\mathit{A}}$ for the toy network
Request a detailed protocolEquation 1 shows that calculation of MI in a cell stateagnostic manner requires the knowledge of the cell stateaveraged response distribution $p\left(B\rightL)$ and the distribution of inputs $p\left(L\right)$. In our calculations, we restrict $p\left(L\right)$ to be a discrete version of the gamma distribution obtained by equal percentile binning of a gamma distribution. The mean of the gamma distribution was taken to be 10, the coefficient of variation 1, and the number of bins was 25. The discretization step was not essential but was taken to simplify the calculations by making all variables discrete.
We assumed that $p\left(B\rightL)$, which is obtained by averaging over the cell state variables, was distributed as a negative binomial distribution (Figure 2—figure supplement 2). We estimated its first two moments by numerically averaging the first two moments of the singlecell response distribution $p\left(B\rightL,\mathit{\theta})$ (Poisson distribution with mean given by Equation S10) according to the population distribution $p\left(\mathit{\theta}\right)$ of the variable parameter(s) $\mathit{\theta}$ ($\mathit{\theta}\equiv {R}_{0}$ for the first cell population and $\mathit{\theta}\equiv {k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ for the second cell population). As mentioned above, $p\left(\mathit{\theta}\right)$ was assumed to be a gamma distribution whose coefficient was systematically varied. Using the first two moments, we inferred the parameters for the negative binomial distribution.
Once the populationlevel response $p\left(B\rightL)$ is obtained and the input distribution $p\left(L\right)$ is fixed, we can calculate ${I}_{CSA}$ using Equation 1. The values of ${I}_{CSA}$ for the population with variable ${R}_{0}$ are given in Figure 2C while those for the population with variable ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ are given in Figure 2—figure supplement 3.
Obtaining ${\mathit{p}}_{\mathit{C}\mathit{e}\mathit{e}\mathit{M}\mathit{I}}\left(\mathit{I}\right)$ and ${\mathit{I}}_{\mathit{C}\mathit{e}\mathit{e}}$ for the toy receptor network
Request a detailed protocolAs indicated by Equations 4 and 6, calculation of the distribution of cell statedependent mutual information values requires the cell statespecific response distribution $p\left(B\rightL,\mathit{\theta})$ and the distribution $p\left(\mathit{\theta}\right)$ of cell states. In the toy model, $p\left(B\rightL,\mathit{\theta})$ is modeled as a Poisson distribution with a mean given by Equation S8 and $p\left(\mathit{\theta}\right)$ is assumed to be gamma distributed in the variable parameter.
Using $p\left(B\rightL,\mathit{\theta})$ and gamma distributed input distribution $p\left(L\right)$, we obtain cell statespecific mutual information $I\left(\mathit{\theta}\right)$ using Equation S2. We find ${p}_{CeeMI}\left(I\right)$ by sampling multiple values of $\mathit{\theta}$ (${R}_{0}$ for the first population and ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ for the second population). ${I}_{Cee}$ is simply the average of ${p}_{CeeMI}\left(I\right)$. ${I}_{Cee}$ and ${p}_{CeeMI}\left(I\right)$ for the population with variable ${R}_{0}$ are given in Figure 2C. Figure 2—figure supplement 3 shows the same for the population with variable ${k}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ .
3. Maximum entropy inference of cell state variability
Using experimentally collected singlecell data on heterogeneity in protein abundances, we estimate the distribution over cell state variables (biochemical parameters) using the maximum entropy approach. Below, we first describe the data that was used in our analysis. Next, we briefly discuss the maximum entropy approach.
EGF/EGFR pathway, data, and model
Request a detailed protocolWe used previously collected singlecell data on cell surface EGFR levels (Dixit et al., 2020; Lyashenko et al., 2020). Briefly, MCF10A cells were stimulated with 10 different extracellular EGF levels ranging between 0 ng/mL and 100 ng/mL. Cell surface EGFR levels were measured in ∼7000 cells for each ligand concentration after 3 hr of continuous EGF stimulation. The data was measured in arbitrary fluorescence units. To convert the data to the units of number of receptors per cell, we used a mean number of cell surface receptors $R=2.5\times {10}^{5}$ for MCF10A cells (Shi et al., 2016). The population mean of the experimentally measured steadystate receptor count distribution in the absence of the ligand was matched to this number. Receptor count data at every other ligand concentration was scaled appropriately.
Using a previously validated model (Dixit et al., 2020; Lyashenko et al., 2020), we constructed a simplified model of the EGF/EGFR pathway. Specifically, we incorporated ligand binding to receptor, receptor activation, and preferential endocytosis of activated receptors. To keep our model simple, we did not incorporate receptor dimerization and oligomerization following ligand binding. Notably, evidence suggests that oligomers may be preformed in the absence of the ligand as well, which could make them effective monomers. Finally, we assumed that the extracellular ligand concentration was kept constant. The model was represented by the following reaction network.
In Equations S11–S14, $R$ is the level of free receptor, $B$ is the level of ligandbound receptor, and $P$ is the level of phosphorylated receptor. The total number of receptors is given by ${R}_{T}=R+B+P.$ Cell state variables $\mathit{\theta}$ comprised $\mathit{\theta}\equiv \left\{{k}_{prod},\phantom{\rule{thinmathspace}{0ex}}{k}_{bind},\phantom{\rule{thinmathspace}{0ex}}{k}_{unbind},\phantom{\rule{thinmathspace}{0ex}}{k}_{p},\phantom{\rule{thinmathspace}{0ex}}{k}_{dp},\phantom{\rule{thinmathspace}{0ex}}{k}_{deg},\phantom{\rule{thinmathspace}{0ex}}{k}_{deg}^{\ast}\right\}.$
IGF/FoxO pathway, data, and model
Request a detailed protocolWe used previously collected singlecell data on nuclear FoxO levels following continuous IGF stimulation in HeLa cells (Gross et al., 2019). Briefly, cells were continuously stimulated with IGF, and nuclear FoxO levels were measured using GFPtagged FoxO using livecell imaging every 3 min for 90 min. To convert the arbitrary fluorescence units to units of copies of FoxO per cell, we first removed the background fluorescence intensity. Then, we used the previously estimated total FoxO levels in HeLa cells (BekkerJensen et al., 2017) (∼710 molecules/cell) and the nucleartocytoplasmic ratio in the absence of stimulation (Wimmer et al., 2014) (2/3 of the total in the nucleus). There was a small disagreement in mean nuclear FoxO levels in the absence of stimulation across different experiments. To remove this artifact and to start all experiments with the average nuclear FoxO levels, we offset individual experiments such that the mean nuclear FoxO was identical across all experiments.
Similar to the EGF/EGFR pathway, using a previously validated model (Wimmer et al., 2014), we constructed a simplified model of the IGF/FoxO pathway. Specifically, we incorporated ligand binding to IGF receptor, receptor activation, activation of Akt, and Aktdriven phosphorylation of FoxO. Phosphorylated FoxO was prohibited from entering the nucleus. Finally, we assumed that the extracellular ligand concentration was kept constant. The model was represented by the following reaction network.
In Equations S15–S21, $R$ is the level of free IGF receptor, $B$ is the level of ligandbound receptor, and $P$ is the level of phosphorylated receptor. $pAkt$ is phosphorylated Akt, $pFox{O}_{c}$ is cytoplasmic phosphorylated FoxO, $Fox{O}_{c}$ is cytoplasmic unphosphorylated FoxO, and $Fox{O}_{n}$ is nuclear unphosphorylated FoxO. Cell state variables $\mathit{\theta}$ comprised $\mathit{\theta}\equiv \left\{{k}_{prod},\phantom{\rule{thinmathspace}{0ex}}{k}_{bind},\phantom{\rule{thinmathspace}{0ex}}{k}_{unbind},\phantom{\rule{thinmathspace}{0ex}}{k}_{p},\phantom{\rule{thinmathspace}{0ex}}{k}_{dp},\phantom{\rule{thinmathspace}{0ex}}{k}_{deg},\phantom{\rule{thinmathspace}{0ex}}{k}_{deg}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{\left[Akt\right]}_{total},\phantom{\rule{thinmathspace}{0ex}}{k}_{ap},\phantom{\rule{thinmathspace}{0ex}}{k}_{adp},\phantom{\rule{thinmathspace}{0ex}}{k}_{in},\phantom{\rule{thinmathspace}{0ex}}{k}_{ef},\phantom{\rule{thinmathspace}{0ex}}{k}_{fp},\phantom{\rule{thinmathspace}{0ex}}{k}_{fdp},\phantom{\rule{thinmathspace}{0ex}}{\left[FoxO\right]}_{total}\right\}.$
Inference of model parameters
Request a detailed protocolWe assume that cells in a population can be assigned a cellspecific state denoted by a state vector $\mathit{\theta}$ that comprises biochemical parameters relevant to the modeled signaling network. The population variability in cell state parameters is represented by the joint probability density $p\left(\mathit{\theta}\right).$ Typically, $p\left(\mathit{\theta}\right)$ is not experimentally accessible. Therefore, we infer it using a previously developed technique called MEDIRIAN (maximum entropybased framework for inference of heterogeneity in dynamics of signaling networks) (Dixit et al., 2020). MERIDIAN infers the maximum entropy distribution $p\left(\mathit{\theta}\right)$ that reproduces a set of averages computed from experimental singlecell measurements.
MERIDIAN requires a mechanistic model of the signaling network that can predict cell’s response to extracellular perturbation (e.g., ligand) and userspecified population averages computed from experimental data. For the EGF/EGFR and IGF/FoxO networks, we used stochastic biochemical models described by Equations S11–S21, respectively. We use the moment closure approximation (Gillespie, 2009) to approximate the singlecell distributions using the first two moments (see below). The differential equations for the pathways can be found on the GitHub.
Now, we briefly describe the MERIDIAN approach. The entropy of any distribution $p\left(\mathit{\theta}\right)$ is given by
In MERIDIAN, we find $p\left(\mathit{\theta}\right)$ that maximizes $S$ while requiring it to reproduce a set of average constraints evaluated using experimental data. Following Dixit et al., 2020, entropy maximized $p\left(\mathit{\theta}\right)$ is given by the Gibbs–Boltzmann distribution
where ${\lambda}_{m}$ are the Lagrange multiplier corresponding to the mth constraint, ${\psi}_{m}\left(\mathit{\theta}\right)$ is the quantity whose average is constrained, and a $\mathrm{\Omega}$ is the normalization constant.
Given a set of constraints (see below), the Lagrange multipliers can be numerically tuned such that the predictions from the distribution $p\left(\mathit{\theta}\right)$ match their experimental value. We optimize the Lagrange multipliers using gradientbased search (Equation S24) using the ADAM algorithm (Kingma, 2014). The gradients for minimizing a Lagrangian cost function are given by Dixit et al., 2020.
In Equation S24, $\u27e8{\psi}_{m}\left(\mathit{\theta}\right)\u27e9}_{\mathit{\theta}$ denotes the ensemble average computed using $p\left(\mathit{\theta}\right)$ and ${R}_{m}$ are the corresponding measurements. We stop the iterative procedure when the mean absolute relative error $\frac{1}{M}\sum _{m}\frac{\left{R}_{m}{\u27e8{\psi}_{m}\left(\mathit{\theta}\right)\u27e9}_{\mathit{\theta}}\right}{{R}_{m}}$ reaches a predefined value.
The predictions from the Max Ent model depend on the choice of the experimental constraints. To ensure that the constraints represent the entire range of singlecell behaviors, we opt for percentile constraints. Specifically, for experimentally collected singlecell data for a given condition (ligand dose, time point, etc.), we first approximate the singlecell histogram using a gamma distribution. Then, we identify abundances that represent 10th–90th percentiles of this distribution. The fraction of cells belonging to each of these percentile windows is exactly 10%. These become our experimental constraints ${R}_{em}=0.1$. Here, ‘e’ denotes the experimental condition (ligand dose, time point, etc.) and $m\in \left[\mathrm{1,10}\right]$ denotes the percentile window.
The corresponding model predictions of the fraction of cells in a given percentile window $\u27e8{\psi}_{em}\left(\mathit{\theta}\right)\u27e9}_{\mathit{\theta}$ are given by
where
In Equation S26, ${p}_{e}\left(r\mathit{\theta}\right)$ is the modelpredicted singlecell distribution of responses (surface EGFR levels or nuclear FoxO levels) for experimental condition $e$. The integration bounds ${l}_{em}$ and ${u}_{em}$ represent the lower and upper bounds of the mth percentile window. The distribution ${p}_{e}\left(r\mathit{\theta}\right)$ in principle can be approximated by several runs of an explicit simulation using Gillespie’s algorithm (Billing et al., 2019). This may prove to be computationally expensive, especially when sampling through multiple parameter sets $\mathit{\theta}$. Therefore, we resort to moment closure techniques (Gillespie, 2009) to approximate ${p}_{e}\left(r\mathit{\theta}\right)$ as a gamma distribution. We used a previously developed package called MOCA (moment closure analysis) (Schnoerr et al., 2015). We used a Gaussian moment closure to obtain the first and the second moments of the distributions. The moment closure approximation was quite accurate compared to the explicit Gillespie simulation and allowed us to rapidly predict single cell response distributions without performing multiple calculations (Figure 2—figure supplement 4).
The averages in Equation S25 cannot be computed analytically. We therefore resort to Markov chain Monte Carlo techniques to approximate them. Briefly, for a fixed set of Lagrange multipliers, we start 150 parallel MCMC chains in the parameter space with a starting point chosen randomly from the previous iteration. Each step in the MCMC calculation attempted to change between 1 and 5 parameters. The step size for the change was a uniform random number whose maximum was 10% of the parameter bounds for individual parameters. Each MCMC chain was run for approximately for $\sim {10}^{4}$ steps for the EGF/EGFR pathway and $\sim 2.5\times {10}^{5}$ for the IGF/FoxO pathway. The first 1000 steps were discarded, and parameter sets were stored every 50th step after that.
We used ADAM to optimize the Lagrange multipliers (Figure 2—figure supplement 5). The hyperparameters for ADAM were as follows: the exponential decay rate for the first moment estimates was set to 0.8. The exponential decay rate for the secondmoment estimates was set to 0.999, the step size was 0.1. This procedure led to a decrease in the relative error. Using the final set of Lagrange multipliers, we sampled parameter sets that represented an in silico cell population. This population was used for further analysis (see below).
4. Numerical quantification of mutual information and channel capacity using livecell imaging data
Request a detailed protocolHere, we describe how we estimated cell statespecific mutual information from livecell imaging data. The IGF/FoxO pathway reaches an approximate steady state within 30–45 min after the introduction of the IGF ligand. We used data collected on HeLa cells where cells were treated with subsequent doses of IGF every 90 min to approximate the singlecell response to multiple IGF levels. We approximated the singlecell steadystate nuclear FoxO distributions as gamma distributions and estimated their means and variances from data collected between 60 and 90 min of IGF stimulation (Figure 2—figure supplement 6). Using these first two moments, we approximated the singlecell response distribution $p\left(r\rightu,\mathit{\theta})$ as a gamma distribution. The mutual information at the singlecell level and the population average cell statespecific mutual information ${I}_{Cee}$ were obtained using these response distributions. For ${I}_{CSA}$ , the first two moments of the singlecell response distributions were averaged to obtain the populationlevel means and variances. This, along with the assumption that the responses were gamma distributions, provided the cell stateagnostic response distribution $p\left(r\rightu)$ from which ${I}_{CSA}$ was obtained.
To verify that cell states are indeed conserved at the time scale of the experiment, we reanalyzed data generated by Gross et al., 2019 wherein cells were perturbed with IGF (37.5 pM), followed by a washout which allowed the cells to reach prestimulation nuclear FoxO levels, followed by a reperturbation with the same amount of IGF. Nuclear FoxO response was measured at the singlecell level after 90 min with IGF exposure both these times. Since the response $x$ to the same input $u$ was measured twice in the same cell (${x}_{1}$ and ${x}_{2}$), we could evaluate the intrinsic variability in response at the singlecell level. We then compared this intrinsic variability to the extrinsic cell statedependent variability in the population.
To do so, we computed for each cell ${\delta =x}_{1}{x}_{2}$ the difference between the two responses. Figure 2—figure supplement 7 shows the histogram $p\left(\delta \right)$ as computed from the data (pink) and the same computed from the model that was trained on the singlecell data (blue). We also computed ${p(\delta}_{0})$, which represented the difference between responses of two different cells from the same population, shown for both data and the model.
As shown in Figure 2—figure supplement 7, the distribution $p\left(\delta \right)$ is significantly narrower than ${p(\delta}_{0})$, suggesting that intracellular variability is significantly smaller than acrosspopulation variability and that cells’ response to the same stimuli is quite conserved, especially when compared to responses in randomly picked pairs of cells. This shows that cell states and the corresponding response to extracellular perturbations are conserved, at least at the time scale of the experiment. Therefore, our estimates of celltocell variability signaling fidelity are stable and reliable.
5. Modelpredicted joint distributions ${p}_{CeeMI}(I,\chi )$ for several biochemical parameters
Request a detailed protocol6. Examining the role of cell state dynamics on sensing ability
Request a detailed protocolTo elucidate the role of cell state dynamics, we built a simple model of ligand–receptor system (Figure 4). We define the system as follows:
Here all cells are assigned identical kinetic parameters. Instead, the state of the system is defined by the abundances of each molecule prior to ligand dose. Thus, each cell may be described by two state variables, $\overrightarrow{\theta}=\{{R}_{0},mRN{A}_{0}\}$, where $mRN{A}_{0}$ is the initial mRNA levels and ${R}_{0}$ is the initial ligandfree receptor count. We tune the mRNA dynamics by changing mRNA production and degradation rates while keeping the average mRNA copy number constant. A relaxation of conditional responses requires both states have sufficiently high turnover, while we will modulate mRNA turnover, we will simply select a sufficiently high receptor turnover rate as described below.
For our simulations, we set $\frac{kmprod}{kmdeg}=5\phantom{\rule{thinmathspace}{0ex}}a.u.$, $k}_{prod}=50\phantom{\rule{thinmathspace}{0ex}}a.u.\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{s}^{1$ , ${k}_{deg}=0.5\phantom{\rule{thinmathspace}{0ex}}{s}^{1},$ and $k}_{unbind}=1\phantom{\rule{thinmathspace}{0ex}}{s}^{1$ . We also set ${k}_{\mathrm{m}\mathrm{d}\mathrm{e}\mathrm{g}}={\tau}^{1}$, where $\tau $ takes on five values, equally distributed in the log space, between 10^{2} and 10^{4} s. We select two distinct cells to study, ${\overrightarrow{\theta}}_{a}=\left\{\mathrm{300,3}\right\}$ and ${\overrightarrow{\theta}}_{b}=\left\{\mathrm{700,7}\right\}$. For each cell and mRNA turnover time scale of interest, we introduce 20 different doses equally distributed in log scale, leading to values of $L{k}_{bind}$ ranging from ${10}^{2}$ to ${10}^{3}\phantom{\rule{thinmathspace}{0ex}}\left({s}^{1}\right)$. We ran 2500 Gillespie simulations for each condition to obtain approximate distributions of the response. We then fit these samples to gamma distributions to obtain conditional responses from which the channel capacity mutual information can be calculated using an approach similar to the one described in Section 1.
Data availability
The current manuscript is a computational study of previously collected data. No data have been generated for this manuscript. Scripts used to generate mutual information results as well the experimental data used to train our models can be obtained at: https://github.com/adgoetz186/Cell_signalling_information (copy archived at Goetz, 2024). Script for maximum entropy inference of cellular parameters can be obtained at: https://github.com/hodaakl/MaxEnt (copy archived at Akl, 2024).
References

An interior point algorithm for largescale nonlinear programmingSIAM Journal on Optimization 9:877–900.https://doi.org/10.1137/S1052623497325107

Heterogeneity of epidermal growth factor binding kinetics on individual cellsBiophysical Journal 73:1089–1102.https://doi.org/10.1016/S00063495(97)781414

Diverse relaxation rates exist among rat cardiomyocytes isolated from a single myocardial regionThe Journal of Physiology 597:711–722.https://doi.org/10.1113/JP276718

Variability within rare cell states enables multiple paths toward drug resistanceNature Biotechnology 39:865–876.https://doi.org/10.1038/s41587021008373

Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypesMolecular Biology of the Cell 24:2186–2200.https://doi.org/10.1091/mbc.E12100737

Momentclosure approximations for massaction modelsIET Systems Biology 3:52–58.https://doi.org/10.1049/ietsyb:20070031

SoftwareCell signalling information, version swh:1:rev:7f892fe7c3ca612f3a0cbf7d8d1643d1598179b7Software Heritage.

Unraveling growth factor signaling and cell cycle progression in individual fibroblastsThe Journal of Biological Chemistry 291:14628–14638.https://doi.org/10.1074/jbc.M116.734194

Correlated receptor transport processes buffer singlecell heterogeneityPLOS Computational Biology 13:e1005779.https://doi.org/10.1371/journal.pcbi.1005779

Livecell imaging in the era of too many microscopesCurrent Opinion in Cell Biology 66:34–42.https://doi.org/10.1016/j.ceb.2020.04.008

Cellular noise and information transmissionCurrent Opinion in Biotechnology 28:156–164.https://doi.org/10.1016/j.copbio.2014.05.002

BookJeremy Thorner Signal Transduction: Principles, Pathways, and Processes (1st ed)New York: Cold Spring Harbor Laboratory Press.

Mathematical modeling of variability in intracellular signalingCurrent Opinion in Systems Biology 16:17–24.https://doi.org/10.1016/j.coisb.2019.10.020

Exploring intermediate cell states through the lens of single cellsCurrent Opinion in Systems Biology 9:32–41.https://doi.org/10.1016/j.coisb.2018.02.009

Cell cycle proliferation decisions: the impact of single cell analysesThe FEBS Journal 284:362–375.https://doi.org/10.1111/febs.13898

Escherichia coli chemotaxis is information limitedNature Physics 17:1426–1431.https://doi.org/10.1038/s41567021013803

Comparison of different momentclosure approximations for stochastic chemical kineticsThe Journal of Chemical Physics 143:18.https://doi.org/10.1063/1.4934990

Singlecell protein analysis by mass spectrometryCurrent Opinion in Chemical Biology 60:1–9.https://doi.org/10.1016/j.cbpa.2020.04.018

Fundamental limits to cellular sensingJournal of Statistical Physics 162:1395–1424.https://doi.org/10.1007/s1095501514405

Defining cell types and states with singlecell genomicsGenome Research 25:1491–1498.https://doi.org/10.1101/gr.190595.115

Physical constraints on accuracy and persistence during breast cancer cell chemotaxisPLOS Computational Biology 15:e1006961.https://doi.org/10.1371/journal.pcbi.1006961

Nongenetic diversity modulates population performanceMolecular Systems Biology 12:895.https://doi.org/10.15252/msb.20167044

Mathematical modeling reveals modulation of both nuclear influx and efflux of Foxo1 by the IGFI/PI3K/Akt pathway in skeletal muscle fibersAmerican Journal of Physiology. Cell Physiology 306:C570–C584.https://doi.org/10.1152/ajpcell.00338.2013

Robustness, accuracy, and cell state heterogeneity in biological systemsCurrent Opinion in Systems Biology 8:46–50.https://doi.org/10.1016/j.coisb.2017.11.009

Singlecell protein analysisCurrent Opinion in Biotechnology 23:83–88.https://doi.org/10.1016/j.copbio.2011.11.023

Distinct cellular states determine calcium signaling responseMolecular Systems Biology 12:894.https://doi.org/10.15252/msb.20167137
Peer review
Reviewer #2 (Public Review):
In this paper the authors present an existing information theoretic framework to assess the ability of single cells to encode external signals sensed through membrane receptors. The main point is to distinguish actual noise in the signaling pathway from cellcell variability, which could be due to differences in their phenotypic state, and to formalize this difference using information theory. After correcting for this cellular variability, the authors find that cells may encode more information than one would estimate from ignoring it, which is expected. The authors show this using simple models of different complexities, and also by analyzing an imaging dataset of the IGF/FoxO pathway.
I am only partially satisfied by the authors response. To me, the main question that was unanswered, while being at the core of the claim of the paper, was the magnitude of withincell variability across repetitions of the stimulus.
This can only be done on the IGF/FoxO system because, as the authors acknowledge, the EGF/EGFR system does not have any data to support any claim about singlecell information that's not heavily informed by models, which assume by construction that this variability is small, naturally leading the desired conclusion.
The authors now measure withincell, acrossrepetition variability (delta_0) for IGF/FoxO, but:
 they compare it to celltocell variability, finding that it's smaller. That's good and that supports the main claim of the paper that single cells are more precise than a mean cell. However they don't show it in the paper, but only in the response.
 they also don't compare it to withincell, withinstimulation variability across time. But this latter variability is what they (wrongly) used to estimate information, and still do in this revision. However I think this is approximated by the blue "simulation" violin plot in Reviewer Figure 2. The true variability is clearly larger than previously assumed. So it's strange that they conclude that "our estimates of celltocell variability signaling fidelity are stable and reliable."
 they don't use this delta_0 variability to revise their estimate of the information accordingly.
 since variability is small compared to the differences between distinct stimulations, of which there are only 4, all information quantities they get are around 2 bits, which is not approaching the information capacity but merely a statement that the number of tested doses is small.
I strongly recommend that the authors actually report the figure they provided as Reviewer Figure 2 in the manuscript. In addition, they should not claim that the withincell variability (captured by the variability across distinct presentations of the stimulus) is well captured by their initial estimate (based on the variance within a single presentation of the stimulus).
https://doi.org/10.7554/eLife.87747.3.sa1Reviewer #3 (Public Review):
Goetz, Akl and Dixit investigated the heterogeneity in the fidelity of sensing the environment by individual cells in a population using computational modeling and analysis of experimental data for two important and wellstudied mammalian signaling pathways: (insulinlike growth factor) IGF/FoxO and (epidermal growth factor) EFG/EFGR mammalian pathways. They quantified this heterogeneity using the conditional mutual information between the input (eg. level of IGF) and output (eg. level of FoxO in the nucleus), conditioned on the "state" variables which characterize the signaling pathway (such as abundances of key proteins, reaction rates, etc.) First, using a toy stochastic model of a receptorligand system  which constitutes the first step of both signaling pathways  they constructed the population average of the mutual information conditioned on the number of receptors and maximized over the input distribution and showed that it is always greater than or equal to the usual or "cell state agnostic" channel capacity. They constructed the probability distribution of cell state dependent mutual information for the two pathways, demonstrating agreement with experimental data in the case of the IGF/FoxO pathway using previously published data. Finally, for the IGF/FoxO pathway, they found the joint distribution of the cell state dependent mutual information and two experimentally accessible state variables: the response range of FoxO and total nuclear FoxO level prior to IGF stimulation. In both cases, the data approximately follow the contour lines of the joint distribution. Interestingly, high nuclear FoxO levels, and therefore lower associated noise in the number of output readout molecules, is not correlated with higher cell state dependent mutual information, as one might expect. This paper contributes to the vibrant body of work on information theoretic characterization of biochemical signaling pathways, using the distribution of cell state dependent mutual information as a metric to highlight the importance of heterogeneity in cell populations. The authors suggest that this metric can be used to infer "bottlenecks" in information transfer in signaling networks, where certain cell state variables have a lower joint distribution with the cell state dependent mutual information.
The utility of a metric based on the conditional mutual information to quantify fidelity of sensing and its heterogeneity (distribution) in a cell population is supported in the comparison with data. Some aspects of the analysis and claims in the main body of the paper and SI need to be clarified and extended.
Remaining Comments:
 I think Review Figure 2 which is currently in the SI would improve the main body of the paper if moved there. In that case, the discussion of this figure in the main text would have to address more than it currently does, namely "the same cell's FoxO responses to the same input were found to have significantly less variation compared to the variation within the population".
https://doi.org/10.7554/eLife.87747.3.sa2Author response
The following is the authors’ response to the original reviews.
Public Reviews:
Reviewer #1 (Public Review):
The manuscript by Goetz et al. takes a new perspective on sensory information processing in cells. In contrast to previous studies, which have used population data to build a response distribution and which estimate sensory information at about 1 bit, this work defines sensory information at the single cell level. To do so, the authors take two approaches. First, they estimate single cells' response distributions to various input levels from timeseries data directly. Second, they infer these singlecell response distributions from the population data by assuming a biochemical model and extracting the cells' parameters with a maximumentropy approach. In either case, they find, for two experimental examples, that singlecell sensory information is much higher than 1 bit, and that the reduction to 1 bit at the population level is due to the fact that cells' response functions are so different from each other. Finally, the authors identify examples of measurable cell properties that do or do not correlate with singlecell sensory information.
The work brings an important and distinct new insight to a research direction that generated strong interest about a decade ago: measuring sensory information in cells and understanding why it is so low. The manuscript is clear, the results are compelling, and the conclusions are well supported by the findings. Several contributions should be of interest to the quantitative biology community (e.g., the demonstration that single cells' sensory information is considerably larger than previously implied, and the approach of inferring singlecell data from population data with the help of a model and a maximumentropy assumption).
We thank the reviewer for the excellent summary of our research.
Reviewer #2 (Public Review):
In this paper the authors present an existing information theoretic framework to assess the ability of single cells to encode external signals sensed through membrane receptors.
The main point is to distinguish actual noise in the signaling pathway from cellcell variability, which could be due to differences in their phenotypic state, and to formalize this difference using information theory.
After correcting for this cellular variability, the authors find that cells may encode more information than one would estimate from ignoring it, which is expected. The authors show this using simple models of different complexities, and also by analyzing an imaging dataset of the IGF/FoxO pathway.
The implications of the work are limited because the analysed data is not rich enough to draw clear conclusions. Specifically,
the authors do not distinguish what could be methodological noise inherent to microscopy techniques (segmentation etc), and actual intrinsic cell state. It's not clear that cellcell variability in the analyzed dataset is not just a constant offset or normalization factor. Other authors (e.g. Gregor et al Cell 130, 153164) have recentered and renormalized their data before further analysis, which is more or less equivalent to the idea of the conditional information in the sense that it aims to correct for this experimental noise.
We thank the reviewer for the comment. However, we do not believe our analysis is a consequence of normalization artifacts. Prior to modeling the single cell data, we removed welldependent background fluorescence. This should take care of technical variation related to overall offsets in the data. We agree with the reviewer that background subtraction may not fully account for technical variability. For example, some of the celltocell variability may potentially be ascribed to issues such as incorrect segmentation. Unfortunately, however, attempting to remove this technical variability through cellspecific normalization as suggested by the reviewer1 will diminish to a very large extent the true biological effects related to extensivity (cell size, total protein abundance). We note that these effects are a direct function of cell statevariables (see for example CohenSaidon et al.2 who use cellstate specific normalization to improve signaling fidelity). Therefore, an increase in mutual information after normalization does not only reflect removal of technical noise but also accounts for effect of cell state variables.
Nonetheless, as the reviewer suggested, we performed a cellspecific normalization wherein the mean nuclear FoxO levels in each cell (in the absence of IGF) were normalized to one. Then, for each ligand concentration, we collated FoxO response across all cells and computed the channel capacity corresponding to cellstate agnostic mutual information ICSA. As expected, ICSA increases from ∼0.9 bits to ∼1.3 bits when cellspecific normalization was performed (Author response image 1). However, this value is significantly lower than the average ∼1.95 of cellstate specific mutual information ⟨ICee⟩. Finally, we note that the cell specific normalization does not change the calculations of channel capacity at the single cell level as these calculations do not depend on linear transformations of the data (centering and normalization). Therefore, we do not think that our analysis of experimental data suffers from artifacts related to microscopy.
in the experiment, each condition is shown only once and sequentially. This means that the reproducibility of the response upon repeated exposures in a single cell was not tested, casting doubt on the estimate of the response fidelity (estimated as the variance over time in a single response).
The reviewer raises an excellent question about persistence of cell states. To verify that cell states are indeed conserved at the time scale of the experiment, we reanalyzed data generated by Gross et al.3 wherein cells were perturbed with IGF (37.5 pM), followed by a washout which allowed the cells to reach prestimulation nuclear FoxO levels, followed by a reperturbation with the same amount of IGF. Nuclear FoxO response was measured at the single cell level after 90 minutes with IGF exposure both these times. Since the response x to the same input u was measured twice in the same cell (x1 and x2), we could evaluate the intrinsic variability in response at the single cell level. We then compared this intrinsic variability to the extrinsic cellstate dependent variability in the population.
To do so, we computed for each cell δ=x1x2 the difference between the two responses. reviewer Figure 2 show the histogram p(δ) as computed from the data (pink) and the same computed from the model that was trained on the single cell data (blue). We also computed p(δ0) which represented the difference between responses of two different cells both from the data and from the model.
As we see in Author response image 2, the distribution p(δ) is significantly narrower than p(δ0) suggesting that intracellular variability is significantly smaller than acrosspopulation variability and that cells’ response to the same stimuli are quite conserved, especially when compared to responses in randomly picked pairs of cells. This shows that cell states and the corresponding response to extracellular perturbations are conserved, at least at the time scale of the experiment. Therefore, our estimates of celltocell variability signaling fidelity are stable and reliable. We have now incorporated this discussion in the manuscript (lines 275281).
another dataset on the EGF/EGFR pathway is analyzed, but no conclusion can be drawn from it because singlecell information cannot be directly estimated from it. The authors instead use a maximumentropy Ansatz, which cannot be validated for lack of data.
We thank the reviewer for this comment. We agree with the reviewer that we have not verified our predictions for the EGF/EGFR pathway. That study was meant to show the potential generality of our analysis. We look forward to validating our predictions for the EGF/EGFR pathway in future studies.
Reviewer #3 (Public Review):
Goetz, Akl and Dixit investigated the heterogeneity in the fidelity of sensing the environment by individual cells in a population using computational modeling and analysis of experimental data for two important and wellstudied mammalian signaling pathways: (insulinlike growth factor) IGF/FoxO and (epidermal growth factor) EFG/EFGR mammalian pathways. They quantified this heterogeneity using the conditional mutual information between the input (eg. level of IGF) and output (eg. level of FoxO in the nucleus), conditioned on the "state" variables which characterize the signaling pathway (such as abundances of key proteins, reaction rates, etc.) First, using a toy stochastic model of a receptorligand system  which constitutes the first step of both signaling pathways  they constructed the population average of the mutual information conditioned on the number of receptors and maximized over the input distribution and showed that it is always greater than or equal to the usual or "cell state agnostic" channel capacity. They constructed the probability distribution of cell state dependent mutual information for the two pathways, demonstrating agreement with experimental data in the case of the IGF/FoxO pathway using previously published data. Finally, for the IGF/FoxO pathway, they found the joint distribution of the cell state dependent mutual information and two experimentally accessible state variables: the response range of FoxO and total nuclear FoxO level prior to IGF stimulation. In both cases, the data approximately follow the contour lines of the joint distribution. Interestingly, high nuclear FoxO levels, and therefore lower associated noise in the number of output readout molecules, is not correlated with higher cell state dependent mutual information, as one might expect. This paper contributes to the vibrant body of work on information theoretic characterization of biochemical signaling pathways, using the distribution of cell state dependent mutual information as a metric to highlight the importance of heterogeneity in cell populations. The authors suggest that this metric can be used to infer "bottlenecks" in information transfer in signaling networks, where certain cell state variables have a lower joint distribution with the cell state dependent mutual information.
The utility of a metric based on the conditional mutual information to quantify fidelity of sensing and its heterogeneity (distribution) in a cell population is supported in the comparison with data. Some aspects of the analysis and claims in the main body of the paper and SI need to be clarified and extended.
1. The authors use their previously published (Ref. 32) maximumentropy based method to extract the probability distribution of cell state variables, which is needed to construct their main result, namely p_CeeMI (I). The salient features of their method, and how it compares with other similar methods of parameter inference should be summarized in the section with this title. In SI 3.3, the Lagrangian, L, and Rm should be defined.
We thank the reviewer for the comment and apologize for the omission. We have now rewritten the manuscript to include references to previous reviews of works that infer probability distributions4 of cell state variables (lines 156168). Notably, as we argued in our previous work5, no current method can efficiently estimate the joint distribution over parameters that is consistent with measured single cell data and models of signaling networks. Therefore, we could not use multiple approaches to infer parameter distributions. We have now expanded our discussion of the method in the supplementary information sections.
1. Throughout the text, the authors refer to "low" and "high" values of the channel capacity. For example, a value of 11.5 bits is claimed to be "low". The authors need to clarify the context in which this value is low: In some physically realistic cases, the signaling network may need to simply distinguish between the present or absence of a ligand, in which case this value would not be low.
We agree with the reviewer that small values of channel capacities might be sufficient for cells to carry out some tasks, in which case a low channel capacity does not necessarily indicate a network not performing its task. Indeed, how much information is needed for a specific task is a related but distinct question from how much information is provided though a signaling network. Both questions are essential to understand a cell's signaling behavior, with the former being far less easy to answer in a way which is generalizable. In contrast, the latter can be quantitatively answered using the analysis presented in our manuscript.
1. Related to (2), the authors should comment on why in Fig. 3A, I_Cee=3. Importantly, where does the fact that the network is able to distinguish between 23 ligand levels come from? Is this related to the choice (and binning) of the input ligand distribution (described in the SI)?
We thank the reviewer for the comment. The network can distinguish between all inputs used in the in silico experiment precisely because the noise at the cellular level is small enough that there is negligible overlap between single cell response distributions. Indeed, the mutual information will not increase with the number of equally spaced inputs in a sublinear manner, especially when the input number is very high.
1. The authors should justify the choice of the gamma distribution in a number of cases (eg. distribution of ligand, distribution cell state parameters, such as number of receptors, receptor degradation rate, etc.).
We thank the reviewer for the comment. We note that previous works in protein abundances and gene expression levels (e.g. see6) have reported distributions with positive skews that can be fit well with gamma distributions or lognormal distributions. Moreover, many stochastic models of protein abundance levels and signaling networks are also known to result in abundances that are distributed according to a negative binomial distribution, the discrete counterpart of gamma distribution. Therefore, we chose Gamma distributions in our study. We have now clarified this point in the Supplementary Information. At the same time, gamma distribution only serves as a regularization for the finite data and in principle, our analysis and conclusion do not depend on choice of gamma distribution for abundances of proteins, ligands, and cell parameters.
1. Referring to SI Section 2, it is stated that the probability of the response (receptor binding occupancy) conditioned on the input ligand concentration and number of receptors is a Poisson distribution. Indeed this is nicely demonstrated in Fig. S2. Therefore it is the coefficient of variation (std/mean) that decreases with increasing R0, not the noise (which is strictly the standard deviation) as stated in the paper.
We thank the reviewer of the comment. We have now corrected our text.
1. In addition to explicitly stating what the input (IGF level) and the output (nuclear GFPtagged FoxO level) are, it would be helpful if it is also stated what is the vector of state variables, theta, corresponding to the schematic diagram in Fig. 2C.
We thank the reviewer of the comment. We have now corrected our text in the supplementary material as well as the main text (Figure 2 caption).
1. Related to Fig. 2C, the statement in the caption: "Phosphorylated Akt leads to phosphorylation of FoxO which effectively shuttles it out of the nucleus." needs clarification: From the figure, it appears that pFoxO does not cross the nuclear membrane, in which case it would be less confusing to say that phosphorylation prevents reentry of FoxO into the nucleus.
We thank the reviewer of the comment. We have now corrected our text (Figure 2 caption).
1. The explanations for Fig. 2D, E and insets are sparse and therefore not clear. The authors should expand on what is meant by model and experimental I(theta). What is CC input dose? Also in Fig. 2E, the overlap between the blue and pink histograms means that the value of the blue histogram for the final bin  and therefore agreement or lack thereof with the experimental result  is not visible. Also, the significance of the values 3.25 bits and 3 bits in these plots should be discussed in connection with the input distributions.
We thank the reviewer of the comment. We have now corrected our text (Figure 2 caption and lines 249251).
1. While the joint distribution of the cell state dependent mutual information and various biochemical parameters is given in Fig. S7, there is no explanation of what these results mean, either in the SI or main text. Related to this, while a central claim of the work is that establishing this joint distribution will allow determination of cell state variables that differentiate between high and low fidelity sensing, this claim would be stronger with more discussion of Figs. 3 and S7. The related central claim that cell state dependent mutual information leads to higher fidelity sensing at the population level would be made stronger if it can be demonstrated that in the limit of rapidly varying cell state variables, the I_CSA is retrieved.
We thank the reviewer for this excellent comment. We have now added more discussion about interpreting the correlation between cell state variables and cellstate specific mutual information (lines 294306). We also appreciate the suggestion about a toy model calculation to show that dynamics of cell state variables affects cell state specific mutual information. We have now performed a simple calculation to show how dynamics of cell state variables affects cells’ sensing ability (lines 325363). Specifically, we constructed a model of a receptor binding to the ligand wherein the receptor levels themselves changed over time through a slow process of gene expression (Author response image 3, main text Figure 4). In this model, the timescales of fluctuations of ligandfree receptors on the cell surface can be tuned by speeding up/slowing down the degradation rate of the corresponding mRNA while keeping the total amount of steady state mRNA constant. As shown in Author response image 3, the dependence of cellspecific mutual information on cell state variable diminishes when the time scale of change of cell state variables is fast.
Reviewer #1 (Recommendations For The Authors):
My major concerns are mainly conceptual, as described below. With proper attention to these concerns, I feel that this manuscript could be a good candidate for the eLife community.
Major concerns:
1. The manuscript convincingly demonstrates that cells good sensors after all, and that heterogeneity makes their inputoutput functions different from each other. This raises the question of what happens downstream of sensing. For singlecelled organisms, where it may be natural to define behavioral consequences at the singlecell level, it may very well be relevant that singlecell information is high, even if cells respond differently to the environment. But for cells in multicellular organisms, like those studied here, I imagine that most behavioral consequences of sensing occur at the multicellular level. Thus, many cells' responses are combined into a larger response. Because their responses are different, their highinformation individual responses may combine into a lowinformation collective response. In fact, one could argue that a decent indicator of the fidelity of this collective response is indeed the populationlevel information measure estimated in previous works. Thus, a fundamental question that the authors must address is: what is the ultimate utility of reliable, but heterogeneous, responses for a multicellular system? This question has an important bearing for the relevance of their findings.
We thank the reviewer for this thoughtprovoking comment. We agree that the fidelity with which cells sense their environment, especially those in multicellular organisms, may not always need to be very high. We speculate that when the biological function of a collection of cells can be expressed as an average over the response of individual cells; highinformation but heterogeneous cells can be considered equivalent to lowinformation homogeneous cells. An example of such a function is population differentiation to maintain relative proportions of different cell types in a tissue or producing a certain amount of extracellular enzyme.
In contrast, we believe that when the biological function involves collective action, spatial patterning, or temporal memory, the difference between reliable but heterogeneous population and unreliable homogeneous population will become significant. We plan to explore this topic in future studies.
1. The authors demonstrate that the agreement is good between their inference approach and the direct estimation of response distributions from singlecell time series data. In fact, the agreement is so good that it raises the question of why one would need the inference approach at all. Is it because singlecell time series data is not always available? Is that why the authors used it for one example and not the other? The validation is an asset, but I imagine that the inference approach is complicated and may make assumptions that are not always true. Thus, its utility and appropriate use must be clarified.
We thank the reviewer for the comment. As the reviewer correctly pointed out, live cell imaging data is not always available and has limited scope. Specifically, optical resolution limits measurements of multiple targets. Moreover, typical live cell measurements measure total abundance or localization and not posttranslational modification (phosphorylation, methylation, etc.) which are crucial to signaling dynamics. The most readily available single cell data such those measured using single cell RNA sequencing, immunofluorescence, or flow cytometry are necessarily snapshots. Therefore, computational models that can connect underlying signaling networks to snapshot data become essential when imputing single cell trajectories. In addition, the modeling also allows us to identify network parameters that correlate most strongly with cellular heterogeneity. We have now clarified this point in the manuscript (lines 366380).
Minor comments:
1. I would point out that the maximum values in the singlecell mutual information distributions (Fig 2D and E) correspond to log2 of the number of inputs levels, corresponding to perfect distinguishability of each of the equallyweighted input states. It is clear that many of the mutual information values cluster toward this maximum, and it would help readers to point out why.
We thank the reviewer for the comment. We have now included a discussion about the skew in the distribution in the text (lines 251260).
1. Line 216 references Fig 2C for the EGF/EGFR pathway, but Fig 2C shows the FoxO pathway. In fact, I did not see a schematic of the EGF/EGFR pathway. It may be helpful to include one, and for completeness perhaps also one for the toy model, and organize the figures accordingly.
We thank the reviewer for the comment. We did not include three separate schematics because the schematics of the EGF/EGFR model and the toy model are subsets of the schematic of the IGF/FoxO model. We have now clarified this point in the manuscript (Figure 2 caption).
Reviewer #2 (Recommendations For The Authors):
the simple model of Fig. 2A would gain from a small cartoon explaining the model and its parameters.
We thank the reviewer for the comment. We did not include a schematic for the toy model as it is a subset of the schematic of the IGF/FoxO model. The schematic of the toy model is included in the supplementary information.
L should be called u, and B should be called x, to be consistent with the rest of the notations in the paper.
We have decided to keep the notation originally presented in the manuscript.
legend of 2E and D should be clarified. "CC input dose" is cryptic. The x axis is the input dose, the y axis is its distribution at the argmax of I. CC is the max of I, not its argmax. Likewise "I" in the legend for the colors should not be used to describe the insets, which are input distributions.
We have now changed this in the manuscript.
the data analysis of the IGF/FoxO pathway should be explained in the main text, not the SI. Otherwise it's impossible to understand how one arrives at, or how to intepret, figure 2E, which is central to the paper. For instance the fact that p(xu,theta) is assumed to be Gaussian, and how the variance and mean are estimated from the actual data is very important to understand the significance of the results.
While we have added more details in the manuscript in various places, for the sake of brevity and clarity, we have decided to keep the details of the calculations in the supplementary materials.
there's no Method's section. Most of the paper's theoretical work is hidden in the SI, while it should be described in the methods.
We thank the review of the comment. However, we believe that adding a methods section will break the narrative of the paper. The methods are described in detail in the supplementary materials with sufficient detail to reproduce our results. Additionally, we also provide a link to the github page that has all scripts related to the manuscript.
PS: please submit a PDF of the SI for review, so that people can read it on any platform (as opposed to a word document, especially with equations)
We have now done this.
Reviewer #3 (Recommendations For The Authors):
1. Subplots in Fig. 1, inset in Fig. 3 are not legible due to small font.
We have now increased the font.
1. Mean absolute error in Fig. S5 and relative error in related text should be clarified.
We have now clarified this in the manuscript.
1. Acronyms (MACO, MERIDIAN) should be defined.
We have now made these changes.
References
1. Gregor T, Tank DW, Wieschaus EF, Bialek W. Probing the limits to positional information. Cell. 2007;130(1):15364. doi: 10.1016/j.cell.2007.05.025. PubMed PMID: WOS:000248587000018.
2. CohenSaidon C, Cohen AA, Sigal A, Liron Y, Alon U. Dynamics and Variability of ERK2 Response to EGF in Individual Living Cells. Mol Cell. 2009;36(5):88593. doi: 10.1016/j.molcel.2009.11.025. PubMed PMID: WOS:000272965400020.
3. Gross SM, Dane MA, Bucher E, Heiser LM. Individual Cells Can Resolve Variations in Stimulus Intensity along the IGFPI3KAKT Signaling Axis. Cell Syst. 2019;9(6):5808 e4.
4. Loos C H, J. Mathematical modeling of variability in intracellular signaling. Current Opinion in Systems Biology. 2019;16:1724.
5. Dixit PD, Lyashenko E, Niepel M, Vitkup D. Maximum Entropy Framework for Predictive Inference of Cell Population Heterogeneity and Responses in Signaling Networks. Cell Syst. 2020;10(2):20412 e8.
6. Taniguchi Y, Choi PJ, Li GW, Chen H, Babu M, Hearn J, Emili A, Xie XS. Quantifying E. coli proteome and transcriptome with singlemolecule sensitivity in single cells. Science. 2010;329(5991):5338. doi: 10.1126/science.1188308. PubMed PMID: 20671182; PMCID: PMC2922915.
https://doi.org/10.7554/eLife.87747.3.sa3Article and author information
Author details
Funding
National Institute of General Medical Sciences (R35GM142547)
 Andrew Goetz
 Hoda Akl
 Purushottam Dixit
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
AG, HA, and PD were supported by NIGMS grant R35GM142547. The authors thank Andre Levchenko for useful discussions.
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure  PSL, France
Reviewing Editor
 AnneFlorence Bitbol, Ecole Polytechnique Federale de Lausanne (EPFL), Switzerland
Version history
 Preprint posted: March 13, 2023 (view preprint)
 Sent for peer review: March 30, 2023
 Preprint posted: June 22, 2023 (view preprint)
 Preprint posted: January 4, 2024 (view preprint)
 Version of Record published: January 31, 2024 (version 1)
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.87747. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, Goetz, Akl et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 282
 Page views

 24
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Microbiology and Infectious Disease
 Physics of Living Systems
Microsporidia are eukaryotic, obligate intracellular parasites that infect a wide range of hosts, leading to health and economic burdens worldwide. Microsporidia use an unusual invasion organelle called the polar tube (PT), which is ejected from a dormant spore at ultrafast speeds, to infect host cells. The mechanics of PT ejection are impressive. Anncaliia algerae microsporidia spores (3–4 μm in size) shoot out a 100nmwide PT at a speed of 300 μm/s, creating a shear rate of 3000 s^{1}. The infectious cargo, which contains two nuclei, is shot through this narrow tube for a distance of ∼60–140 μm (Jaroenlak et al, 2020) and into the host cell. Considering the large hydraulic resistance in an extremely thin tube and the lowReynoldsnumber nature of the process, it is not known how microsporidia can achieve this ultrafast event. In this study, we use Serial BlockFace Scanning Electron Microscopy to capture 3dimensional snapshots of A. algerae spores in different states of the PT ejection process. Grounded in these data, we propose a theoretical framework starting with a systematic exploration of possible topological connectivity amongst organelles, and assess the energy requirements of the resulting models. We perform PT firing experiments in media of varying viscosity, and use the results to rank our proposed hypotheses based on their predicted energy requirement. We also present a possible mechanism for cargo translocation, and quantitatively compare our predictions to experimental observations. Our study provides a comprehensive biophysical analysis of the energy dissipation of microsporidian infection process and demonstrates the extreme limits of cellular hydraulics.

 Computational and Systems Biology
 Physics of Living Systems
The adaptive dynamics of evolving microbial populations takes place on a complex fitness landscape generated by epistatic interactions. The population generically consists of multiple competing strains, a phenomenon known as clonal interference. Microscopic epistasis and clonal interference are central aspects of evolution in microbes, but their combined effects on the functional form of the population’s mean fitness are poorly understood. Here, we develop a computational method that resolves the full microscopic complexity of a simulated evolving population subject to a standard serial dilution protocol. Through extensive numerical experimentation, we find that stronger microscopic epistasis gives rise to fitness trajectories with slower growth independent of the number of competing strains, which we quantify with powerlaw fits and understand mechanistically via a random walk model that neglects dynamical correlations between genes. We show that increasing the level of clonal interference leads to fitness trajectories with faster growth (in functional form) without microscopic epistasis, but leaves the rate of growth invariant when epistasis is sufficiently strong, indicating that the role of clonal interference depends intimately on the underlying fitness landscape. The simulation package for this work may be found at https://github.com/nmboffi/spin_glass_evodyn.