Simple biochemical networks allow accurate sensing of multiple ligands with a single receptor

Cells use surface receptors to estimate concentrations of external ligands. Limits on the accuracy of such estimations have been well studied for pairs of ligand and receptor species. However, the environment typically contains many ligands, which can bind to the same receptors with different affinities, resulting in cross-talk. In traditional rate models, such cross-talk prevents accurate inference of concentrations of individual ligands. In contrast, here we show that knowing the precise timing sequence of stochastic binding and unbinding events allows one receptor to provide information about multiple ligands simultaneously and with a high accuracy. We show that such high-accuracy estimation of multiple concentrations can be realized with simple structural modifications of the familiar kinetic proofreading biochemical network diagram. We give two specific examples of such modifications. We argue that structural and functional features of real cellular biochemical sensory networks in immune cells, such as feedforward and feedback loops or ligand antagonism, sometimes can be understood as solutions to the accurate multi-ligand estimation problem.


Introduction
Cells obtain information about their environment by capturing ligand molecules with receptors on their surface and estimating the ligand concentration from the receptor activity. Limits on the accuracy of such estimation have been a subject of interest since the seminal work of Berg and Purcell [1], with several substantial extensions found recently [2][3][4][5][6][7][8]. Most of these assume one ligand species coupled to one receptor species, and the actual detection in most of these models is rather simple, involving counting the number or the duration of binding / unbinding events over a specific period of time. However, cells carry many types of receptors and have many species of ligands around them. The same ligand can bind to many receptors, albeit with different affinities, and vice versa. This is commonly referred to as cross-talk. At the same time, real cellular sensory systems are incredibly complex, involving many dozens of identified biochemical species downstream of a typical receptor [9]. Functionally many of such signaling motifs are probably related to solving the cross-talk problem [10,11], and are a topic of active research.
In traditional deterministic chemical kinetics, one cannot estimate concentrations of more ligands than there are receptor types. Further, even a weak cross-talk prevents determination of concentrations of individual chemical species since the activity of a receptor is a function of a weighted sum of concentrations of all ligands that can bind to it. In contrast, here we argue that, with cross-talk, concentration of more than one chemical species can be inferred from the activity of one receptor, provided that the stochastic temporal sequence of receptor binding and unbinding events is accessible instead of its mean occupancy. This is an important departure from the traditional view of cellular signaling that posits as many receptor types as there are ligand concentrations to be estimated. Indeed, previous works studying temporal sequences of receptor occupancy for ligand detection [11] and concentration estimation [5,13,12] have only considered the detection/estimation of a single ligand present in a mixture. We argue that the receptor occupancy sequence contains much more information about the mixture. In fact, based on the maximum likelihood techniques, which have been used previously to study receptor occupancy, we show that all components of the ligand mixture can be estimated by just one receptor, at least in principle. This surprising result can be understood by noting that a typical duration of time that a ligand remains bound to the receptors depends on its unbinding rate. Thus observing the statistics of the receptor's unbound time durations allows estimation of a weighted average of all chemical species that interact with it [5]. Then the statistics of the bound time durations tells how common each ligand is.
The result is very general and independent on the choice of a downstream biochemical kinetics scheme that actually performs the estimation. In this article, we derive the result for the simplest problem of this class, namely one receptor interacting with two ligand species. While the exact solution of the inference problem for finding both ligand concentrations is hard to implement using common biochemical machinery, we show that an accurate approximation is possible using simple extensions of the familiar kinetic proofreading mechanism [14,15]. We identify examples of such motifs implementing such estimation of multiple concentrations in signaling networks found downstream of many immune receptors [9], arguing that real biological systems may be implementing such multivariate concentration sensing. The kinetic schemes that we analyze detect rare ligands more accurately than a simple kinetic proofreading does, and we argue that the involved biochemical computation can explain properties like ligand antagonism, commonly observed in receptor signaling.
Overall, these different arguments support our main idea, that the temporal sequence of binding and unbinding on a single receptor can provide an accurate estimate of the concentration of multiple ligands that bind to the receptor, and that the involved calculations can be performed reliably by known biochemical networks.

The model
Consider a single receptor interacting with a cognate and a non-cognate ligand (Fig 1) that have the concentrations c c and c nc , respectively. The binding rate of the ligands to the receptor are k c and k nc . The binding rates are diffusion limited and hence k c *k nc . It is the unbinding or off-rates, r c and r nc , that distinguish the two ligands: r nc > r c , and a cognate molecule typically stays bound for longer. The binding and unbinding rates (k α 's and r α 's) are fixed and can be assumed known for each receptor-ligand pair. Thus we are interested in the estimation of the ligand concentrations only, c c and c nc . Following Ref. [5], we estimate c c and c nc from the timeseries of binding, ft b i g, and unbinding, ft u i g, events of a total duration T using Maximum Likelihood techniques, paralleling a recent similar independent discussion, which focused on detection of a single ligand concentration [12]. The numbers of binding and unbinding events are different by, at most, one, which is insignificant since we consider T ! 1. Thus without loss of generality, we assume that the first event was a binding event at t b 1 , and the last one was the unbinding at t u n . We write the probability distribution of observing the sequence ft b 1 ; t u 1 ; . . . ; t b n ; t u n g, or alternatively the sequence of binding and unbinding intervals Here the first term under the product sign is the probability of the receptor staying unbound for t u i . The second term, which from now on we denote by Dðc c ; c nc ; t b i Þ, is proportional to the probability of staying bound for t b i . Dðc c ; c nc ; t b i Þ has contributions from binding events from both the cognate and the noncognate ligands, with odds of c c and c nc , respectively. Finally, Z is the normalization, where the sum is over all sequences of duration T and n binding-unbinding events. Note that here we define t u n ¼ t b 1 þ ðT À t u n Þ, so that the n'th unbound interval includes the "incomplete" unbound intervals before the first binding and after the last unbinding.
The log-likelihood of c c and c nc is the logarithm of P, Eq (1). Taking the derivatives of the log-likelihood w. r. t. c c and c nc and setting them to zero gives the Maximum Likelihood (ML) equations for the concentrations. Denoting by T u ¼ P n i¼1 t u i , the total time the receptor is . Two ligands, cognate and non-cognate having concentrations c c and c nc , bind to a receptor R with binding rates k c and k nc , respectively. The cognate unbinding rate is defined as lower than the non-cognate one (r c < r nc ). (b) Time series of receptor occupancy is used to determine both on-rates.
where c Ã c and c Ã nc denotes the ML solution. Multiplying Eqs (3) and (4) by c Ã c and c Ã nc , respectively, and adding them gives As in Ref. [5], the total on-rate (the weighted average of the external concentrations) is determined only by the average duration of the unbound interval, (n/T u ) −1 , because no binding is possible when the receptor is already bound. For the special case of k c % k nc % k (for ligands with binding rate determined by diffusion), Eq (5) determines the maximum likelihood estimate of the sum of the two concentrations, similar to the result in Ref. [5,12]: This shows that the estimates are negatively correlated. For general k i 's, a weighted sum of the concentrations is determined, but the negative correlation persists.
To get the individual concentrations, we need to solve the ML equations Eqs (3) and (4). In general, they can only be solved numerically. However, as all ML estimators, they are unbiased to the leading order in n (we verified this numerically). The standard errors of the ML estimates can be obtained by inverting the Hessian matrix, where greek indices stand for {c, nc}. Each term in the Hessian matrix is a sum of n numbers, each smaller than zero. The inverse of @ 2 log P @c a @c b , which scales as / 1/n, sets the minimum variance of any unbiased estimator according to the Cramer-Rao bound. It has straightforward analytical approximations in various regimes. For example, when the noncognate ligand is almost absent (c c /c nc ) 1), and its few molecules do not bind for long (r c /r nc ( 1), one gets , matching the accuracy of sensing one ligand with one receptor [5]. A regime relevant for detection of a rare, but highly specific ligand [11,12,16] can be investigated as well. For now, we focus on how the receptor estimates (rather than detects) concentrations of both ligands simultaneously, which requires us to explore the full range of on-and off-rates.
The estimates of the concentration c c and c nc are obtained by numerically solving ML equations, Eqs (3) and (4). We study the variability of these ML estimators in terms of their posterior variances. Notice that these posterior variances scale as 1/n, so we define the error of the ML estimators, E, as the squared coefficient of variation times the number of binding-unbinding events, n. Hence, we have, E c ¼ ns 2 ðc Ã c Þ=c 2 c and E nc ¼ ns 2 ðc Ã nc Þ=c 2 nc for cognate and noncognate ligands, respectively. These quantities have a finite limit at n ! 1. Specifically, E = 1 is the accuracy that a receptor that binds only a single ligand can obtain [5]. Thus E c and E nc compare the performance of our multi-ligand ML estimator to the limit achievable by a single ligand ML estimator. We show log 10 E c and log 10 E nc for different concentrations and off-rates in Fig 2. If the two ligands are readily distinguishable, r c ( r nc , then the ligand with the larger concentration has E * 1. When c c * c nc , E i * 4. . .5, and it grows to 10. . .30 for a ligand with a very small relative concentration. Emphasizing the importance of the time scale separation, E > 100 if the ligands are hard to distinguish, r c * r nc . Here the correlation coefficient ρ of the two estimates reaches −1 because the same binding event can be attributed to either ligand. Finally, the asymmetry of the plots w. r. t. the exchange of c c and c nc is because the cognate ligand can generate short binding events, while long events from the noncognate ligand are exponentially unlikely. In summary, it is possible to infer two ligand concentrations from one receptor, with the error of only 1. . .10 times larger than for ligand-receptor pairs with no cross talk, as long as the two off-rates are substantially different. This complements the findings of [12] that a single concentration can be inferred from a time series of "on" and "off" events in a background of noncognate bindings using Maximum Likelihood estimation. We have verified that the analytical expression for the estimation error derived in Ref. [12] for a single cognate ligand matches our numerical results (see Methods).

Approximate solution
It is not clear if there exist biochemical networks that can solve the ML equations, Eqs (3) and (4), exactly. Luckily, an approximate solution exists. Note that most of the long binding events come from the cognate ligand since the noncognate one dissociates faster. Defining long events as t b i ! T c and using Eq (5), we rewrite Eq (3) as Assuming that all long events are cognate, T c ) 1/r nc , gives  where n l is the number of long events, and the superscript "a" stands for the approximate solution. If further T is long enough so that there are many short events, and a single binding duration hardly affects k Ã c , then the sum in Eq (9) can be approximated by the expectation value: where Pðt b jc a c ; c a nc Þ is the probability of observing a binding event of duration τ b for the given binding rates, Plugging Eq (11) into Eq (10), we obtain Finally, since n l ( n, using Eq (5), we get (see Methods for a detailed derivation): In other words, the approximate cognate ligand concentration is proportional to the number of long events. We can estimate the bias and the variance of c a c and c a nc in a limiting case. If r c and r nc are not very different from each other, then one needs to focus on extremely long events in order to identify cognate bindings. This is only possible if T c is much larger than the inverse of both of the unbinding rates, T c ) fr À 1 nc ; r À 1 c g. Large T c ensures that the long binding events get no or minimal contribution from non-cognate ligands. However, since the time for which the receptor stays bound is exponentially distributed, under this condition, the number of "long" events (such that τ b > T c ) would be very small, n l ( n. Thus most of the variance of c a c and c a nc in Eqs (13) and (14) comes from the variability of n l , but not T u (since Further, the individual unbound periods are independent, so that hT u i = nhτ u i = n/(k c c c + k nc c nc ) (notice the use of c rather than c a here). Further, Thus for large T c , the bias of the approximate estimator, k nc c nc k c e À ðr nc À r c ÞT c , grows with the relative number of noncognate long bindings events. In turn, the latter is proportional to c nc , but decreases exponentially with T c .
Within the same approximation, the variance of the estimator is given by However, long binding events are rare, independent of each other, and hence obey the Poisson statistics. Thus σ 2 (n l ) = hn l i, so that The variance obviously grows with T c .
Knowing that the bias and the variance of the approximation change in opposite directions with T c , we can find the optimal cutoff (T c Ã ) by minimizing the overall error. We define such error L as the sum of the variance and the squared bias of the estimator, i. e., The optimal cutoff is obtained by minimizing L c or, in other words, solving the bias-variance tradeoff: Near the optimal cutoff, the bias is small, and we use c c instead of c a c for the variance of the estimator, Eq (16). Then solving Eq (19) gives: Plugging this into Eqs (15) and (16), we get the minimal error of the estimator, which we omit here for brevity. The optimal cutoff is proportional to 1/r nc if r nc ) r c , and it grows with r c , allowing for better disambiguation of cognate and noncognate events. Crucially, the off-rates are dictated by the ligand identities. In contrast, the concentrations, c c and c nc , are what the receptors measures. Therefore, it is encouraging that T c Ã depends only logarithmically on the concentrations (and also on the duration of the measurement, T u ). Thus even if T c is fixed as T c Ã at some fixed values of c c , c nc , it remains near-optimal for a broad range of external concentrations. To illustrate this, we use T c ¼ T c Ã ðk c c c ¼ k nc c nc ¼ 1=2Þ T 0 and analyze the quality of the  Simple biochemical networks allow accurate sensing of multiple ligands with a single receptor approximation in Fig 3, where we plot the ratio L c ðT 0 Þ=s 2 c c and L nc ðT 0 Þ=s 2 c nc . Notice that s 2 c c and s 2 c nc , the variances of the exact ML estimators, are proportional to E c and E nc , respectively. Since ML estimators are unbiased, the ratios L c ðT 0 Þ=s 2 c c and L nc ðT 0 Þ=s 2 c nc compare the errors of the approximate solution to the errors E c and E nc . Since these ratios approach 1 when r c /r nc ! 0 (specifically, for r c /r nc = 0.1, L c ðT 0 Þ=s 2 k c % 1:47, and L nc ðT 0 Þ=s 2 k nc % 1:21), we conclude that the approximation is accurate even at fixed T c = T 0 when its assumptions are satisfied. This happens even though T c Ã depends on c c and c nc , but apparently the approximate estimates are as good as the ML estimates even at fixed T c = T 0 and work well for a large range of concentration ratios. This is important, as the molecular mechanisms that sets the delays in the cell does not need to be modified for different ligand concentrations.
In contrast, when the ligands are nearly indistinguishable (r c /r nc * 1), both L c ðT 0 Þ=s 2 k c $ 100 and L nc ðT 0 Þ=s 2 k nc $ 100, but here one would not use one receptor to estimate two concentrations since even the ML solution is bad (cf. Fig 2). Note also that both L c and L nc are smaller for r c * r nc if c c ) c nc . This is because our main assumption (that almost all long events are cognate) holds better when cognate ligands dominate. Finally, the correlation coefficient between the approximate estimates, ρ a (right panel) reaches -1 earlier than in Fig 2. This is a direct consequence of Eqs (13) and (14).

Kinetic proofreading for approximate estimation
The approximate solution can be computed by cells using the well-known kinetic proofreading (KPR) mechanism [14,15,17,18]. In the simplest model of KPR [19], intermediate states between an inactive and an active state of a receptor delay the activation. Thus bound ligands can dissociate before the receptor activates, at which point it quickly reverts to the inactive state. Since r c < r nc , cognate ligands dominate among bindings that persist to activation. The resulting increase in specificity in various KPR schemes has led to their exploration in the context of detection of rare ligands [11,12,16,18]. Instead, here we analyze their ability to measure concentrations of both ligands simultaneously. We first consider the case where both the cognate and the non-cognate ligand concentration are comparable, c c * c nc and the dissociation rates are distinct, r c ( r nc . In the following sections, we explore another case, c c ( c nc and r c ≲ r nc , a situation common in immunology. Consider a biochemical network in Fig 4(a): the receptor, R, activates two messenger molecules, A and B. The former is activated with the rate k A only if the receptor stays bound for longer than a certain T c (with the delay achieved using the KPR intermediate states). The latter is activated with the rate k B whenever the receptor is bound. The molecules deactivate with the rates r A and r B , respectively, and all activations/deactivations are first-order reactions. Then the mean concentrations of the messenger molecules are (see Methods): Assuming again that most bindings longer than T c are cognate (T c ) 1/r nc ), Eq (21), can be written as: Further, it is easy to see that Eq (22) can be rewritten as: Now solving Eqs (23) and (24) for the on-rates, we get The corrections of the form " B=ðk B =r B À " BÞ appear because bindings only happen to unbound receptors, as emphasized in Ref. [5]. However, these nonlinear relations are still hard to implement with simple biochemical components. We solve this by further assuming ¼ " B=ðk B =r B Þ ( 1, which is true if the receptor is mostly unbound, which happens at low  concentrations. This gives These equations are analogous to Eqs (13) and (14). They are easy to realize biochemically (cf. Fig 4(a)): c c is related to the concentration of the proofread species A by a rescaling, and c nc comes from subtracting rescaled versions of B and A from each other. The subtraction can be done by the third species C, activated by B and suppressed by A. Since ( 1, then " A and " B are small, and many such activation-suppression schemes are linearized as the subtraction [8]. Note that such incoherent feedforward loops (the receptor activates A and B, which then affect C incoherently by suppressing and activating it, respectively) are ubiquitous in cellular networks downstream of receptors [9].
The bias of c a c and c a nc due to long, but noncognate binding events, Eq (15) [19,20]. This variability changes the rate of occurrence of long biding events, but they are still rare, nearly independent, and Poisson-distributed. Denoting by hÁi the averaging at a fixed T c , and by Á " the averaging over T c , we get where we have used the approximation " T c ) 1=r nc in the last step. Thus s 2 T c effectively renormalizes the cutoff to " T c À 1 2 r c s 2 T c . Replacing T c in Eqs (27) and (28) by its renormalized value, which is an easy change in the scaling factors, removes this additional bias due to the random T c in the KPR scheme.
Since long bindings are rare, the variance of the KPR estimator is dominated again generally by " A, but not " B. The intrinsic stochasticity in the production of molecules of A contributes to the variance. However, this contribution can be made arbitrarily small by increasing k A , and we neglect it here. A larger contribution comes from the random number of long bound intervals and a random duration of each of them. To calculate this, in the limit of rare long binding events, we use well-known results in the theory of noise propagation in chemical networks [21] This is a direct analog of Eq (16). In principle, one can measure more than two concentrations similarly, as long as all species have distinct off-rates. For example, to estimate three concentrations, one needs an additional branch downstream of the receptor that proofreads for an intermediate time. Then the branches with the strongest, intermediate, and no proofreading would measure approximately the highest affinity ligand, a combination of the two higher affinity ligands, and all three ligands, respectively. Appropriate activation and inhibition of downstream targets will then allow identifying individual concentrations from these combined readouts. The error (the variance of the ML estimator, and both the bias and the variance for the approximate and the KPR estimators) would grow with an increasing number of ligand species, largely because a larger range of off-rates would be required to disambiguate more ligands. However, this would still represent a dramatic increase in the information gained by a receptor that tracks its precise temporal dynamics, rather than just the average binding state.

Using precise timing to disambiguate two similar ligands
Here we depart slightly from our scenario and show how a KPR-based scheme relying on the entire temporal sequence of activation / deactivation events can estimate the concentration of a single cognate ligand even if the two ligands have very similar off-rates r c ≲ r nc , a situation common in immunology. In such a situation, the KPR branch gets activated not just by the cognate ligand, but also by the non-cognate ligand (though at a smaller rate). When the goal is the accurate estimation of the cognate ligand only, then the contribution to the KPR branch by the non-cognate ligand needs to be removed. To construct a signal transduction network able to do this, we abstract from the existing detailed model of FcRI immunological receptor [9], a well studied eukaryotic signal transduction system mediating many allergic reactions [22]. Here the main signaling branch gets activated through the Lyn-Syk kinase pathway following kinetic proofreading after a ligand binds to the receptor [9]. However, receptor binding excites an additional branch early on, after only one step in kinetic proofreading (a single phosphorylation on the β chain of the receptor). This branch activates Inpp5d (SHIP) phosphotase, which later dephosphorylates Phosphatidylinositol 3-phosphate (PIP3), a key downstream output of the main signaling branch, and sequesters the dephosphorylated product PtdIns(3, 4)P 2 [9]. The part of this signaling motif relevant for our analysis is summarized in a deliberately simplified signaling diagram in Fig 4(b), where A stands for PtdIns(3, 4, 5)P 3 (PIP3), I stands for PtdIns(3, 4)P 2 , and I is produced by SHIP. Further, R is the FcRI receptor bound to an antibody, and cognate and noncognate molecules are the antigens specific/nonspecific to the bound antibody.
In this network, we consider the main activator branch (A), activated after the usual KPR delay, and hence sensitive to long binding events only (which now have contributions both from k c and k nc ). The secondary inhibiting branch (I) is activated by many more binding events, though the shortest, nonspecific background binding events may be removed from both branches by additional proofreading steps (an early cross-phosphorylation event in the FcRI system). The messengers in both branches later form a complex AI, and only A not in the complex activates further downstream signaling. If the production rates of A and I are appropriately matched (which can be done if the off-rates are known a priori, which they should be for such a molecular signal detection system), this sequestration of A by I can effectively remove the contribution to the A branch coming from the non-cognate ligand. The kinetic diagram can be described with the following rate equations (where, for simplicity, we neglect the first proofreading common to both branches): where r A/I are the degradation rates of the messengers A and I, r AI is the sequestration rate, and β A/I are the messenger production rates, derived as above: Here k A/I are the rates of production of A and I, respectively, when the receptor has been bound for a sufficiently long time to produce either. We assume for simplicity r A = r I . Further, we choose r A = r I ( r AI A * r AI I, so that sequestration rather than degradation is primarily responsible for the disappearance of the messengers. Then the steady state solution of the rate equations (Eqs 31 and 32) is [23] (see Methods): " The numerators of both β A and β I are linear combinations of c c and c nc . If the parameters of the biochemical networks are such that the production rate of the proofread branch is ð1þk c c c =r c þk nc c nc =r nc Þ k I > 0, which has a c nc -independent numerator. Thus the contribution of non-cognate ligand to the activator branch is largely sequestered. Moreover, for large r AI , we have " A ss / ð1 þ k c c c =r c þ k nc c nc =r nc Þ À 1 , so that the activation of the A branch decreases as c nc increases. In contrast, if c c = 0 (no cognate ligands present), then , which grows with c nc . This behavior is reminiscent of the agonistantagonist picture in FcRI receptor activation [24]: a weak ligand by itself can activate the cellular response, but it inhibits (antagonizes) activation of the response by a stronger agonist if both are present.

Discussion
The realization of Refs. [5,13,12] and others that the detailed temporal sequence of binding and unbinding events carries more information about the ligand concentration than the mean receptor occupancy is a conceptual breakthrough. It parallels the realization in the computational neuroscience community that precise timing of spikes carries more information about the stimulus than the mean neural firing rate [25][26][27][28][29][30], and it has a potential to be equally impactful. This extra information when measuring one ligand concentration with one receptor [5,12] amounted to increasing the sensing accuracy by a constant prefactor, or, equivalently, getting only a finite number of additional bits from even a very long measurement [31]. In contrast, here we show that two concentrations can be measured with one receptor with the variance that decreases inversely proportionally to the number of observations, n, Eq (16), or to the integration time, 1/r B , Eq (30), so that the accuracy is only an (often small) prefactor lower than would be possible with one receptor per ligand species. Asymptotically, this doubles the information obtained by the receptor [31]. Crucially, such improvement would not be possible without the cross-talk, or binding among noncognate ligands and receptors. Normally, the cross-talk is considered a nuisance that must be suppressed [32,33]. Instead, we argue that cross-talk can be beneficial by recruiting more receptor types to measure the concentration of the same ligand. In particular, this allows having fewer receptor than ligand species, potentially illuminating how cells function reliably in chemically complex environments with few receptor types. Further, the cross-talk can increase the dynamic range of the entire system: a ligand may saturate its cognate receptor, preventing accurate measurement of its (high) concentration, but it may be in the sensitive range of non-cognate receptors at the same time. Finally, the increased bandwidth may lead to improvements in sensing a time-dependent ligand concentration [11,13]. In forthcoming publications, we plan to explore such many-to-many sensory schemes, extending ideas of Ref. [34] to tracking temporal sequences of activation of the receptor and to temporally varying environments.
While the exact maximum likelihood inference of multiple concentrations from a temporal binding-unbinding sequence is rather complex, we showed that when the cognate and the non-cognate off-rates are substantially different, there is a simpler, approximate, but accurate inference procedure for joint measurements of cognate and noncognate ligands. And even if the off-rates are close, one can still measure the cognate ligand concentration reliably. Crucially, this inference can be performed by biochemical motifs readily available to the cell. Namely, one needs two branches of activation downstream of the receptor, with at least one of them having a kinetic proofreading (KPR) time delay. Then the individual ligand concentrations can be obtained by mutual inhibition between the two branches, or by incoherent feedforward loops. We emphasize again that this allows estimation of multiple concentrations from activity of a single receptor, in contrast to a better estimation of just one concentration [12].
Our simple models only illustrate a wide class of models that can use the temporal structure of the receptor binding sequence to measure more that one ligand concentrations for various ligand combinations, including similar and dissimilar ligands. Additional branches from different points in the proofreading cascade provide additional information about the binding affinities of the mixture of ligands present in the environment, and then algebraic operations on these readouts can be performed by a large diversity of feedforward and feedback loops, competitions for the substrate and the enzyme, and so on. For example, in our simple model, the action of the antagonist is due to the competition for available receptors, while experiments suggest competition for a critical initiating kinase [24], which would require a straightforward modification of the model. Similarly, antagonists are usually "medium" affinity ligands, while very weak ligands do not antagonize receptors. As illustrated in Fig 4(b), this can be achieved by having an additional KPR time delay common to both A and I branches, which occurs in practice [9].
The kinetic diagram for the FcRI receptor is not unique, and similar (though not equivalent) structures exist for other immune cells and receptors as well [9]. Such common structural features result in a similar phenomenology of activation profiles, which are different for pure ligands and ligand mixtures, and depend nontrivially on the details of the binding affinities and concentrations of the ligands in the mixture [10,16,[35][36][37][38]. Interestingly, on longer time scales, a potentially related phenomenon in innate immune response is that of endotoxin tolerance (desensitization to commonly present ligands) [39], which also affects ligands of different affinity differently, and in this case also depends on the history of exposure to other ligands [40]. It is mediated by SHIP, a crucial player in our analysis of FcRI signaling [41], whose activity may be interpreted as setting the relative gain on the A and I branches of Fig 4(b), thus resulting in a more accurate signal estimation. In other words, one interpretation of the known results is that, as various feedback loops increase the activity of SHIP in response to frequent activation of signaling downstream of the receptor, the amount of I increases, thus sequestering more A, lowering its steady-state activity, and inducing tolerance. An important contribution of the understanding developed here is that one can try to interpret these various kinetic diagrams and their phenomenological consequences as implementing estimation of concentrations of potentially many ligands (rather detection of a single one [11,13,16]), and maybe even doing it in a (nearly) Maximum Likelihood optimal fashion, under various assumptions about the number of distinct ligands, their relative abundance, and the (dis) similarity of the off-rates. Exploring feasibility of such an interpretation is an additional interesting venue for future research.
In summary, monitoring precise temporal sequences of receptor activation/deactivation opens up new and exciting possibilities for environment sensing by cells.

Methods
Here we provide mathematical derivations of some of the steps ommitted in the Results.

Derivation of maximum likelihood equations
We start with: The log-likelihood of k c,nc is the logarithm of P: Taking the derivatives of the log-likelihood w. r. t. c c and c nc and setting them to zero gives the Maximum Likelihood (ML) equations for the concentrations. These are: Here, Dðc Ã c ; c Ã nc ; t b i Þ ¼ ðk c c Ã c r c e À t b i r c þ k nc c Ã nc r nc e À t b i r nc Þ, with Ã denoting the ML solution. Denoting by T u ¼ P n i¼1 t u i , the total time for which the receptor is unbound, these equations can be rewritten as Multiplying Eqs (41) and (42) by c Ã c and c Ã nc , respectively, and adding them gives Comparison of simulations with analytical results for single concentration estimation Here we compare the results obtained from the numerical simulations to the analytical expressions derived in Ref. [12] for detection of the concentration of the cognate ligand in a background of spurious ligands. The variance of the concentration estimation obtained from the simulations matches quite well with integral expression, Eq. (7) in Ref. [12], Fig 5a. Note that this expression is inverse of the (1, 1) term of the Hessian matrix, Eq 7. The analytical results obtained for the low concentration of the cognate ligand compared to a non cognate ligand (c c ( c nc ) also match the simulations , Fig 5b. Finally, using Eq 43 we get c a nc ¼ for longer than a certain T c (with the delay achieved using the KPR intermediate states). The latter is activated with the rate k B whenever the receptor is bound. The molecules deactivate with the rates r A and r B , respectively, and all activations/deactivations are first-order reactions. The rate equation for the two molecules can be written as: The Θ functions represent the fact that A is produced only when the receptor has been bound for longer than the cutoff time T c , and B is produced only when the receptor is bound. The steady state value of " A can be obtained by equating the average deactivation rate r A " A to k A times the fraction of time the receptor occupancy was larger than the cutoff, T c , i.e., Similarly, " B can be obtained as: Therefore, the mean concentrations of the messenger molecules are: " B ¼ k c c c =r c þ k nc c nc =r nc 1 þ k c c c =r c þ k nc c nc =r nc k B r B : ð59Þ Using precise timing to disambiguate two close ligands: Derivation of Eqs (35) and (36) The rate equations are: Equating the r. h. s. to zero gives the steady state conditions: b A À r A " A ss À r AI " A ss " I ss ¼ 0; ð62Þ b I À r I " I ss À r AI " A ss " I ss ¼ 0: The latter of these can be rewritten as: Plugging this in Eq (62), we get which can be simplified to: This quadratic equation has the solution: Now sssuming r A = r I and r A = r I ( r AI A * r AI I, we get: One can similarly can get the equation for " I ss as well.