Phenotypic spandrel: absolute discrimination and ligand antagonism

Recent works in quantitative evolution combined with mathematical modelling have shown that evolution of biological networks is constrained by selected phenotypes in unexpected ways. Trade-offs between different functionalities increasingly appear as major forces shaping evolution of complex phenotypes moving on evolutionary Pareto fronts [1, 2]. Numerical experiments of in silico evolution of phenotypic models of gene networks [3] have further shown that, surprisingly, the selection of complex phenotypes leads to the appearance of complex traits that have not been explicitly selected for. Such phenomena are reminiscent of the architectural image of 'evolutionary spandrels' proposed by Gould and Lewontin [4]. They argued that many biological properties are necessary by-products of more fundamental adaptive traits, due to underlying constraints (e.g. tridimensional geometry in the case of architectural spandrels). Spandrels can themselves be tinkered by evolution into new functional structures, leading to the notion of 'exaptation' [5]. But we are still lacking a quantitative theory of such spandrels, which might explain that other scholars have questioned the notion of spandrels (see a summary in Gould's own rebuttal [6]).

A biological example of broad interest is the absolute discrimination between molecules of different 'qualities', one instance being early immune recognition of antigens by T cells [7–10]. T cells differentiate between self and not-self ligands based on a measure of some effective biochemical parameter μ that is characteristic of ligand quality. This quality has first been suggested to be defined by the typical binding time τ [11, 12] of ligands to T cell receptors (TCRs), defining the so-called 'life-time dogma' [8]. For higher association rates, the affinity K_D has also been proposed to be the parameter discriminating between self and not-self [13–15]. Discrimination has to be extremely specific to the quality μ (e.g. to prevent an auto-immune disease): insensitive to the very high number of self ligands with $\mu \lt {\mu }_{c}$, but also very sensitive to low ligand concentrations (e.g. to detect a nascent infection), thus the term 'absolute discrimination' [16]. Recent works in immunology have confirmed that discrimination by T cells is almost absolute, and furthermore shown how it can be modulated by external cytokines such as IL-2 [17].

In evolutionary simulations aiming at reverse-engineering absolute discrimination [18], evolved networks always present undesirable ligand antagonism. Antagonism is a 'dog in the manger' effect, due to the cross-talk between different ligands. Like the dog in Aesop's fable who cannot eat the grain but nevertheless prevents the horse from eating it as well, antagonist ligands that are unable to trigger a response actively prevent agonists from doing so. Such effects have been described experimentally in several immune recognition processes [19, 20]. It is intriguing that antagonism similar to nature spontaneously appears in evolutionary simulations without explicit selection. In both simulated and real networks, antagonism is a consequence of crucial negative interactions for specific and sensitive ligand detection [18, 20–22]. We could show that the main effect of kinetic proofreading in this system is to lower antagonism [18], but is it possible to construct an absolutely discriminating system without any antagonism? Below, we show from general arguments how ligand antagonism is a necessary consequence of absolute discrimination, thus qualifying as a 'phenotypic spandrel'. Our approach allows us to characterize generic solutions of absolute discrimination, to exhibit minimum networks and families of absolute discrimination models, and to further show how addition of proofreading steps leads to an 'inverted' hierarchy of antagonism that is not a generic feature of the antagonism spandrel. Finally, we generalize this theory by making connections with other biochemical examples, including endocrine signalling.

We consider a family of biochemical ligands with different quantitative properties encoded by a continuous parameter μ, subsequently called 'ligand quality'. Ideal absolute discrimination is performed when cells successfully discriminate between two categories of ligands, above and below critical quality ${\mu }_{c}$, irrespective of ligand quantity L. An idealized response diagram for absolute discrimination in the $(L,\mu )$ plane defines a vertical 'discrimination line' at ${\mu }_{c}$ (figure 1 (A)). We model absolute discrimination using biochemical networks with continuous, single-valued variables, at steady state, in a way similar to phenotypic models such as those reviewed in [9]. Ligands interact with receptors at the surface of the cell; we make a mean-field approximation so that all concentrations are averaged out over one cell, and all variables are assumed to be continuous and differentiable (due to biochemistry).

Zoom In Zoom Out Reset image size

We assume that a cellular network is responsible for discrimination between ligand qualities, so that the decision is eventually based on thresholding some internal variable T (which could for instance be the steady state value of a downstream protein concentration, or some total phosphorylation state). We however require T to be a differentiable function of parameters because of biochemistry. We will call $T(L,\mu )$ the value of T when L ligands of single quality μ are presented. While L is in principle a discrete variable, it can be treated as a continuous one for our calculations without loss of generality.

In vivo, cells are exposed to complex mixtures of ligands interacting with receptors at the surface of a single cell. We will use the notation $T(\{{L}_{n},{\mu }_{n}\})$ for the output value in the presence of a mixture, where L_n is the number of ligands presented of quality ${\mu }_{n}$ , and again insist on differentiability of this variable. The quality of ligands are thus indexed by n, e.g. n = 1, 2 if only two types of ligands are presented, L₁ ligands of quality ${\mu }_{1}$ and L₂ ligands of quality ${\mu }_{2}$. We will limit ourselves to a mixture of two ligands but our reasoning can be easily generalized.

We start by connecting the mixture case $T(\{{L}_{n},{\mu }_{n}\})$ to the pure case, assuming ${\mu }_{n}={\mu }_{c}+{\rm{d}}{\mu }_{n}$ with very small ${\rm{d}}{\mu }_{n}$. To build some intuition about the result, let us start with an ensemble of L identical ligands at ${\mu }_{c}$ and vary the quality of one single ligand presented to the cell by an infinitesimal quantity ${\rm{d}}{\mu }_{1}$. This should translate into an infinitesimal change of the output variable $T(L,{\mu }_{c})$ at linear order due to change of quality ${\rm{d}}{\mu }_{1}$ of a single ligand, defining the quantity ${ \mathcal A }$:

$$\begin{eqnarray}&&T(\{L-1,{\mu }_{c};1,{\mu }_{c}+{\rm{d}}{\mu }_{1}\})=T(L,{\mu }_{c})+{ \mathcal A }{\rm{d}}{\mu }_{1}.\end{eqnarray} \tag{ 1 }$$

${ \mathcal A }$ is a function of L and ${\mu }_{c}$, and is well defined in a mean-field approximation where all ligands are equivalent and T is differentiable as required.

Now let us consider another individual ligand. If we vary its quality by an infinitesimal quantity ${\rm{d}}{\mu }_{2}$, of same order of magnitude of ${\rm{d}}{\mu }_{1}$, at linear order this adds another contribution to T which is equal to ${ \mathcal A }(L,{\mu }_{c}){\rm{d}}{\mu }_{2}$. This is the same ${ \mathcal A }(L,{\mu }_{c})$ because all ligands are equivalent, and even though the first ligand changed of quality by ${\rm{d}}{\mu }_{1}$, this only infinitesimally changes all variables in the system by terms of order ${\rm{d}}{\mu }_{1}$ so that at lowest order there is no change for ${ \mathcal A }$.

We can then generalize this reasoning to any number of ligands: each infinitesimal variation of quality ${\rm{d}}{\mu }_{n}$ of one single ligand gives an equivalent infinitesimal contribution ${ \mathcal A }(L,{\mu }_{c}){\rm{d}}{\mu }_{n}$ at linear order, so that considering a mixture $\{{L}_{n},{\mu }_{n}\}$ we get directly (factorizing ${ \mathcal A }$)

$$\begin{eqnarray}&&T(\{{L}_{n},{\mu }_{n}={\mu }_{c}+{\rm{d}}{\mu }_{n}\})\\ &&\quad =T(L,{\mu }_{c})+{ \mathcal A }(L,{\mu }_{c})\displaystyle \sum _{n}{L}_{n}{\rm{d}}{\mu }_{n}\end{eqnarray} \tag{ 2 }$$

calling $L={\sum }_{n}{L}_{n}$. This is a completely generic result for mixtures that does not depend on the fact that the system is doing absolute discrimination.

As a simple example, let us assume a simple ligand–receptor system far from saturation, where ligands bind to receptors with binding time ${\mu }^{-1}$, with a local kinetic amplification mechanism (e.g. proofreading [23]), so that ${T}_{n}=\kappa {L}_{n}{\mu }_{n}^{2}$ is the average number of receptors bound to the ligands of quality ${\mu }_{n}$, assuming L_n ligands of this quality are presented, and some reaction rate κ (for simplicity we take $\kappa =1$ in the following). Let us take as an output the total number of receptors bound $T={\sum }_{n}{T}_{n}$. Then following our reasoning, starting with L ligands at ${\mu }_{c}$ so that $T=L{\mu }_{c}^{2}$, changing the quality of one ligand by ${\rm{d}}{\mu }_{1}$ gives a total number of bound receptors $T={(L-1){\mu }_{c}^{2}+1({\mu }_{c}+{\rm{d}}{\mu }_{1})}^{2}\simeq L{\mu }_{c}^{2}+2{\mu }_{c}{\rm{d}}{\mu }_{1}$ at linear order, i.e. with our notation, ${ \mathcal A }=2{\mu }_{c}$. Changing the quality of another ligand by ${\rm{d}}{\mu }_{2}$, we get $T={(L-2){\mu }_{c}^{2}+1({\mu }_{c}+{\rm{d}}{\mu }_{1})}^{2}+1{({\mu }_{c}+{\rm{d}}{\mu }_{2})}^{2}$ $\simeq L{\mu }_{c}^{2}+2{\mu }_{c}({\rm{d}}{\mu }_{1}+{\rm{d}}{\mu }_{2}).$ This reasoning can clearly be generalized to get equation (2). While this model looks extremely simple, the more complex examples with feedback presented in what follows work in a similar way.

Coming back to our general reasoning, we focus now on absolute discrimination, where the detection of ligand quality μ is specific and done independently from the total ligand concentration L. We assume some response is triggered if $T\geqslant {\rm{\Theta }}$. For L ligands at $\mu ={\mu }_{c}$ from figure 1(A), by continuity we necessarily have on the discrimination line:

$$\begin{eqnarray}&&\forall L\quad T(L,{\mu }_{c})={\rm{\Theta }}.\end{eqnarray} \tag{ 3 }$$

This means that T is independent from the concentration of ligand L when $\mu ={\mu }_{c}$, or in other words is biochemically adaptive with respect to variable L [24, 25]. Further assuming that the decision is made if $\mu \gt {\mu }_{c}$, and as a consequence of T becoming higher than the threshold, we necessarily have $T(L,{\mu }_{c}+{\rm{d}}\mu )\gt {\rm{\Theta }}=T(L,{\mu }_{c})$. Using equation (2) with a single type of ligand, we thus get ${ \mathcal A }(L,{\mu }_{c})\gt 0$.

Consider now the problem of mixtures of ligands of two kinds, sketched in figure 1(B). From equation (2) and positivity of ${ \mathcal A }$, if all ${\rm{d}}{\mu }_{n}$ are positive, then $T(\{{L}_{n},{\mu }_{c}+{\rm{d}}{\mu }_{n}\})\gt T(L,{\mu }_{c})={\rm{\Theta }}$, which in plain words means that a mixture of agonists (close to the threshold) always triggers a response. A more interesting case is to consider what happens with two different types of ligands, L₁ at ${\mu }_{c}$ and L₂ at ${\mu }_{c}-{\rm{d}}\mu$. From 2, we immediately get

$$\begin{eqnarray}&&T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})\\ &&=T({L}_{1}+{L}_{2},{\mu }_{c})-{L}_{2}{\rm{d}}\mu { \mathcal A }({L}_{1}+{L}_{2},{\mu }_{c}).\end{eqnarray} \tag{ 4 }$$

Now from equation (3), we have $T({L}_{1}+{L}_{2},{\mu }_{c})={\rm{\Theta }}=T({L}_{1},{\mu }_{c})$, and since ${ \mathcal A }\gt 0$ we thus have $T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})-T({L}_{1},{\mu }_{c})$ $=-{L}_{2}{\rm{d}}\mu { \mathcal A }({L}_{1}+{L}_{2},{\mu }_{c})\lt 0$, so that we get our main result:

$$\begin{eqnarray}&&T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})\lt T({L}_{1},{\mu }_{c}).\end{eqnarray} \tag{ 5 }$$

This expression establishes ligand antagonism: a mixture of sub-threshold ligands L₂ with critical agonist ligands L₁ yields lower signalling variable $T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})$ than the case where the same quantity L₁ of critical ligands are presented (response $T({L}_{1},{\mu }_{c})$). Thus if a decision is based on the thresholding of output T, the response disappears. Antagonism is therefore established as a general property of systems performing absolute discrimination based on one parameter μ, completely irrespective of internal biochemistry.

An intuitive explanation of the above reasoning can be made by comparison with general models for signalling that do not perform absolute discrimination. Let us consider the behaviour of an output variable T of a given signalling pathway, before any thresholding-based decision. We study how its position varies with respect to a reference level, where only one type of ligand is presented but where the T level is such that the cell responds to the external signal. For many models of signalling pathways with independent receptors, we expect that the output function T is monotonic in both L and μ. For instance, in kinetic proofreading proposed in [23], each receptor contributes additively to signalling once it is bound, so that more ligands, or ligands with longer binding times, necessarily give stronger signals. In such models, starting from a critical ligand concentration L₁ triggering a response, any addition of ligands L₂ thus gives a higher output T value, and thus if the response is based on the thresholding of T, the response is maintained. Furthermore, addition of L₂ ligands with critical quality ${\mu }_{c}$ gives a higher output than addition of the same quantity L₂ of ligands with lower quality $\mu \lt {\mu }_{c}$ (figure 1(C)). Considering now absolute discrimination and again addition of L₂ ligands on figure 1(C), the constraint encoded by adaptation from equation (3) means that addition of extra L₂ ligands of quality $\mu ={\mu }_{c}$ does not change the value of output T. But if discrimination is based on μ, we expect that output T still is a monotonic function of μ. Thus, from this point with ${L}_{1}+{L}_{2}$ ligands at ${\mu }_{c}$, if ligand quality of the extra L₂ ligands is lowered ($\mu \lt {\mu }_{c}$), we still expect a lower contribution of those ligands to the output T, just like the monotonic example (figure 1(C)). Thus if we started right at the threshold for T, the response is now below the threshold of activation, corresponding to antagonism.

Antagonism thus appears as a direct consequence of biochemical adaptation at $\mu ={\mu }_{c}$. As long as such adaptive behaviour is observed and variable T is differentiable, antagonism ensues, irrespective of the details of biochemistry (such as receptor complexations, non-linearities in networks, etc.). Counterexamples can be built if the differentiation hypothesis does not hold: for instance if the cell could measure the maximum binding time of individual ligands, which is clearly not a differentiable function, then there would be no antagonism. Equation (2) further tells us that antagonism occurs in ligand mixtures at linear order as soon as $\sum {L}_{n}{\rm{d}}{\mu }_{n}\lt 0$. Of course it will also be observed for some range of ${\rm{d}}\mu$ even at non-linear order. Furthermore, adaptation does not necessarily have to be perfect: non-monotonic response curves [26], with T varying around Θ, essentially approximate well the necessary adaptation for absolute discrimination and display antagonism [18, 22]. Mathematically, the position with respect to the threshold depends on the competition between the 'flatness' of $T(L,{\mu }_{c})$ as a function of L and the mixture term ${ \mathcal A }(L,{\mu }_{c}){\sum }_{n}{L}_{n}{\rm{d}}{\mu }_{n}$. Examples of this effect are given in the next section.

We now exhibit and study two interesting classes of models performing absolute discrimination but with different antagonistic behaviour.

A 'homeostatic model' (figure 2) is inspired by a ligand–receptor adaptive network evolved in [24], and implements both absolute discrimination and linear antagonism at all orders as described above. Receptors are produced with fixed rate (rescaled to 1 without loss of generality). Receptor–ligand complexes (D_i) can form irreversibly with association rate κ but are degraded inside the cell with time-scale ${\mu }_{i}$, defining quality of ligands (see figure 2(A)). The output is the total sum of ligand–receptor complexes. The steady state equations of such a homeostatic model for the mixture $\{({L}_{n},{\mu }_{n})\}$ are

$$\begin{eqnarray}&&R=\displaystyle \frac{1}{\kappa {\displaystyle \sum }_{n}{L}_{n}}\quad T=\displaystyle \sum _{n}{D}_{n}=\kappa R\displaystyle \sum _{n}{\mu }_{n}{L}_{n}\end{eqnarray} \tag{ 6 }$$

where we defined the output variable T as the sum of all possible ligand–receptor complexes ${\sum }_{n}{D}_{n}$. When only one type of ligand $(L,{\mu }_{c})$ is present, $T={\sum }_{n}{D}_{n}={\mu }_{c}$ is adaptive (i.e. independent from L) as expected. At steady state, R is inversely proportional to the total ligand concentration presented irrespective of their μ. Importantly, this is not due to a titration effect of a fixed pool: R dynamically buffers any ligand addition, so that steady state output T is a weighted linear combination of ${\mu }_{n}$ in equation (6). Such a combination is always inferior to the maximum of ${\mu }_{n}$, so if the maximum indeed is ${\mu }_{c}$, antagonism ensues. If we consider the mixture of L₁ ligands (critical time ${\mu }_{c}$) with L₂ ligands (${\mu }_{c}-d\mu$) we get

$$\begin{eqnarray}&&T={\mu }_{c}-{\rm{d}}\mu \displaystyle \frac{{L}_{2}}{{L}_{1}+{L}_{2}}\quad {\rm{i}}.{\rm{e}}.\ { \mathcal A }(L,{\mu }_{c})={L}^{-1}.\end{eqnarray} \tag{ 7 }$$

This contribution is exact and linear in ${\rm{d}}\mu$ at all orders. Note that the R variable biochemically encodes the antagonistic strength ${ \mathcal A }$ (modulo the κ prefactor that sets up time-scale but compensates in steady state expression).

Zoom In Zoom Out Reset image size

Figures 2(B) and (C) show values of T for pure ligands and a mixture, with vertical discrimination lines at ${\mu }_{c}$. This simple model thus represents a perfect implementation of the principles sketched in figure 1. Surprisingly, absolute discrimination for this model works at steady state even in the limit $L\to 0$, however it should be noted from the equations in figure 2 that the speed of convergence slows down proportionally to L in this limit, so that for small L one essentially never reaches steady state.

Antagonism strength is proportional to the deviation of antagonists' quality ${\rm{d}}\mu$ from ${\mu }_{c}$ (equation (7)) . In particular, antagonism is very weak for ligand quality close to ${\mu }_{c}$ and gets stronger as quality of ligands gets further below the threshold. Interestingly the hierarchy of antagonism with increasing antagonism for lower quality μ is actually completely opposite to antagonism observed in the immune examples. It is known there that antagonism is maximized for ligands close to critical quality, and vanishes for very small binding time, corresponding to self [20]. It is a biological necessity that self ligands, which are the most frequent ones, do not antagonize immune response (otherwise the immune system would not be able to detect any foreign agent), so it is not clear a priori how to reconcile this with the constraint that absolute discrimination implies antagonism, and a special study in such a context is thus required.

Many immune models are based on kinetic proofreading, first proposed in this context by McKeithan [23]. This model was based on the observation that one of the main differences between self and not-self ligands is their binding time τ to the T cell receptors. Kinetic proofreading in this context assumes that a ligand receptor complex can undergo subsequent steps of phosphorylations, but that upon release of the ligand, the phosphorylation state of the receptor is reinitialized. It has been well known since the original proposal of kinetic proofreading by Hopfield and Ninio [27, 28] that each of these phosphorylation steps can contribute geometrically up to a factor τ, so that the number of the last complex in a proofreading cascade with N steps, in the absence of saturation, is roughly proportional to $L{\tau }^{N}$ at steady state, with L the number of ligands presented. Following the observation of sequential phosphorylations on internal tails of immune receptors [11], many current models for immune recognition [20, 22, 29, 30] have elaborated on these first ideas developed by McKeithan, see also [9] for a review/comparison between models.

The 'life-time dogma' [8] posits that immune recognition is almost an absolute discrimination process where the quality of ligand μ is encoded by the binding time of ligands τ [17]. This notion first came from a series of experimental data suggesting that the typical ligand concentration to trigger a response falls by at least 4 to 5 orders of magnitude with a moderate change of binding time τ [12, 31]. Various theoretical and experimental studies of this behaviour have led authors to propose that kinetic proofreading should be complemented with combinations of internal positive/negative feedbacks to explain such observed high specificity in τ [8, 16, 20].

We start with a general derivation for proofreading-based models realizing absolute discrimination, to connect them to biochemical adaptation. Receptors exist in different phosphorylation states (denoted C_j, figure 3(A), where j indicates the step in the phosphorylation cascade). The transition rates between states are assumed to be functions of internal variables (denoted ${\bf{M}}$), accounting for all signalling inside cells (kinases, phosphatases). These effectors are assumed to diffuse freely and rapidly, so that each receptor sees the same value for their concentrations. ${\bf{M}}$ values depend on the total occupancy of some C_j. The downstream decision is usually assumed to be effectively controlled by thresholding on one state, index t in the proofreading cascade (i.e. $T={C}_{t}$). Using standard assumptions, it is easy to show (appendix) that for such proofreading based models, in the limit of unsaturated receptors, the number of receptor states C_j bound to ligands $(L,\tau )$ takes the functional form:

$$\begin{eqnarray}&&{C}_{j}(L,\tau )={c}_{j}({\bf{M}},\tau )L\end{eqnarray} \tag{ 8 }$$

i.e., we can deconvolve influence of ligands into a linear contribution (L), and an indirect one (c_j) purely due to internal variables.

Zoom In Zoom Out Reset image size

In the context of absolute discrimination, with quality $\mu =\tau$, equation (3) imposes that $T={C}_{t}$ is tuned to threshold Θ irrespective of L at critical quality ${\tau }_{c}$. Equation (8) thus implies that the indirect contribution c_t to output variable T verifies that

$$\begin{eqnarray}&&{c}_{t}({\bf{M}},{\tau }_{c})={\rm{\Theta }}/L.\end{eqnarray} \tag{ 9 }$$

c_t is a pure function of internal variables that decreases as a function of L. But from equation (8), c_t contributes multiplicatively to the signal. This means that internal variables are used to compensate for the 'direct' linearity of L (from equation (8)), thus necessarily implementing an incoherent feedforward or feedback loop [32] via c_t (which here plays a role similar to R in the homeostatic model above).

This is exactly what happens in the so-called adaptive sorting model (figure 3(B), see full equations in the appendix). In this model, a kinase K responsible for one of the phosphorylation steps of the cascade is negatively regulated by an earlier step in the cascade. Because of this, effective phosphorylation rate mediated by K takes the functional form of equation (9) for high enough L, and consequently the concentration of the last complex of the cascade is a pure function of τ [18]. A decision can then be made by thresholding on C_N (vertical asymptotic discrimination line in figure 3(B)).

Another model with similar behaviour, inspired by the T cell immune recognition network [22], is displayed in figure 3(C) (see full equations in the appendix). While adaptive sorting relies on repression of an internal kinase K, this model instead uses activation of a phosphatase (variable S) to provide the incoherent feedback term, allowing for adaptation similar to equation (9). One caveat though is that because of the non-specificity of the phosphatase that acts on all steps in the cascade, adaptation is not perfect, but still antagonism is observed (figure 3(C) , see [22] for an explicit study).

An important difference from the homeostatic model is that the transition rates between variables (and thus c_j) depend explicitly on ligand binding time τ in proofreading models. It is then informative to consider a slightly more general ansatz with only one internal variable M and separable influence of ligand quality μ, so that the total output variable from a mixture of ligands $\{{L}_{n},{\mu }_{n}\}$ is

$$\begin{eqnarray}&&T(\{{L}_{n},{\mu }_{n}\})=\displaystyle \sum _{n}c({\mu }_{n},M){L}_{n}={M}^{\beta }\displaystyle \sum _{n}\,f({\mu }_{n}){L}_{n}\end{eqnarray} \tag{ 10 }$$

with $M={\sum }_{n}g({\mu }_{n}){L}_{n}$ at steady state for mixture $\{{L}_{n},{\mu }_{n}\}$. $\beta =-1$ gives perfect adaptation, but different values of β can give realistic biological effects such as loss of response at high ligand concentration, as observed in [22]. We get by direct differentiation:

$$\begin{eqnarray}&&{ \mathcal A }(L,{\mu }_{c})={L}^{\beta }{\left.\displaystyle \frac{{\rm{d}}}{{\rm{d}}\mu }({{fg}}^{\beta })\right|}_{\mu ={\mu }_{c}}.\end{eqnarray} \tag{ 11 }$$

This term is very similar to equation (7), with an extra μ dependency coming from $f,g$. If we impose that a mixture of agonists triggers a response, we necessary have ${ \mathcal A }\gt 0$. If $f\propto {\mu }^{t+1}$ and $g\propto {\mu }^{m+1}$, which is approximately the case for adaptive sorting (see the appendix), we thus have by substituting in equation (11), $t+1+\beta (m+1)\gt 0$, and so for $\beta =-1$ we get $t\gt m$. In many immune recognition models [18, 20, 22], this constraint is naturally realized because the internal variable (kinase or phosphatase) implementing approximate biochemical adaptation is regulated by a much earlier step (m) in the same proofreading cascade than the output (t). This shows how the structure of proofreading cascades with additional feedback constrains both discrimination and antagonism. Furthermore, as soon as $\beta \lt 0$, in a similar way to c_t above, ${ \mathcal A }$ encodes a negative feedforward effect.

To understand what happens for $\mu \ll {\mu }_{c}$, it is more illuminating to compute a ratio of responses (calling r the fraction of total ligands L with quality $\mu \lt {\mu }_{c}$)

$$\begin{eqnarray}\rho & = & \displaystyle \frac{T(\{(1-r)L,{\mu }_{c};{rL},\mu \})}{T(L,{\mu }_{c})}\\ & = & \left(1-r+\displaystyle \frac{f(\mu )}{f({\mu }_{c})}r\right){\left(1-r+\displaystyle \frac{g(\mu )}{g({\mu }_{c})}r\right)}^{\beta }.\end{eqnarray} \tag{ 12 }$$

In models of immune detection based on proofreading, such as adaptive sorting, f and g are powers of quality $\mu =\tau$ (binding time), so that $f(\mu )/f({\mu }_{c}),g(\mu )/g({\mu }_{c})\ll 1$ when $\mu \lt {\mu }_{c}$. So if $\beta =-1$, we see that $\rho \to 1$ when $\mu /{\mu }_{c}\to 0$, explaining why antagonism disappears in this parameter regime (corresponding to self in immune detection). Proofreading models interpolate between no antagonism for $\mu =\tau \to 0$ and linear antagonism close to threshold. Thus antagonism strength increases as τ increases from 0, as observed experimentally [20], before quickly collapsing again very close to the threshold. Quality of ligands τ for maximum antagonism is closer to the threshold ${\tau }_{c}$ with increasing m and t − m (figure 4).

Zoom In Zoom Out Reset image size

Interestingly, far from threshold, the behaviour of antagonism with respect to quality μ is a power of $m+1$, i.e., the negative internal contribution g dominates, while close to threshold, the power t − m dominates, which quantifies the competition between the direct positive (f) and indirect negative (g) contributions to the signal. So practically, the measurement of antagonism for different μ provides a way to quantify direct and indirect contributions to the response, which suggests an experimental strategy to quantify feedbacks inside such systems via measurements of antagonism.

Finally, if β is not exactly equal to −1, the system is not perfectly adaptive, but the positivity of ${ \mathcal A }$ still implies antagonism of similar magnitude to models with perfect adaptation. Assuming a power-law dependency so that $T\propto {L}^{\epsilon }$ close to $\mu ={\mu }_{c}$, then one can show (see the appendix) that upon addition of sub threshold ligands with quality ${\mu }_{c}-{\rm{d}}\mu$, antagonism appears as soon as ${\rm{d}}\mu$ is of order $\epsilon$ (in proper units). This can be seen in the model of [22], illustrated in figure 4(A) 'immune model'. For ligand quality close to the rescaled quality 1, there is a synergistic effect (represented here by negative antagonism) due to the small increase of response upon the addition of ligands of the same quality ($\epsilon$ power-law). However, as soon as the quality of the L₂ ligands is sufficiently low, we see a decrease in response linear in ${\rm{d}}\mu$, until one reaches the antagonistic regime when the (negative) effect of lowering ligand quality (${ \mathcal A }{\rm{d}}\mu$) dominates the (positive) effect of ligand concentration increase $\left(\epsilon \mathrm{log}\left(1+\tfrac{{L}_{2}}{{L}_{1}}\right)\right)$.

We have established how antagonism is a necessary by-product of strict specific and sensitive detection, thus qualifying as a 'phenotypic spandrel'. While there have been recent studies of spandrels at the level of protein structure [33], to our knowledge, this is the first study of an explicit spandrel structuring a complex signalling phenotype.

With the exception of [34], our modelling hypotheses correspond to most phenotypic models of immune ligand discrimination we are aware of [9], as well as more complex models such as the one in [20]. Noise can be included as in [18, 22] via time-integration of some output variable, replacing all concentrations by expectation values, e.g. $C\to \langle C\rangle$. To possibly disentangle antagonism from underlying biochemical adaptation, one would need to examine limiting cases or complexify hypotheses. One could consider for instance the time-course of a variable far from steady-state [34] or coupled to multi stability [35] to 'break' the differentiability hypothesis.

Antagonism appears as a direct consequence of biochemical adaptation. Further relaxing assumptions to look for alternative mechanisms, it is useful to write the Taylor expansion of output $T(L,\mu )$ for pure ligands close to $\mu ={\mu }_{c}$:

$$\begin{eqnarray}&&{\rm{d}}{T}=\displaystyle \frac{\partial T}{\partial L}{\rm{d}}{L}+\displaystyle \frac{\partial T}{\partial \mu }{\rm{d}}\mu .\end{eqnarray} \tag{ 13 }$$

As said before, strict specificity imposes that for $\mu ={\mu }_{c}$, $T={\rm{\Theta }}$ independently of L, which mathematically implies from equation (13) that for ${\rm{d}}\mu =0$, dT = 0 so that ${\left.\tfrac{\partial T}{\partial L}\right|}_{\mu ={\mu }_{c}}=0$, which is biochemical adaptation as assumed throughout this manuscript. Imperfect adaptation means non zero $\tfrac{\partial T}{\partial L}$, but if it is small enough, we saw previously that we still expect antagonism. Notice in this case from equation (13) that dT a priori stays small as both dL and ${\rm{d}}\mu$ vary. When we relax the strict specificity hypothesis so that dT = 0 for small but non zero ${\rm{d}}\mu$, another possible limit appears, most visible by computing the slope of the discrimination line (for dT = 0):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{L}}{{\rm{d}}\mu }=-\displaystyle \frac{\tfrac{\partial T}{\partial \mu }}{\tfrac{\partial T}{\partial L}}.\end{eqnarray} \tag{ 14 }$$

To get a vertical discrimination line in the $(L,\tau )$ plane, the slope needs to go to $\infty ,$ so one could either take $\tfrac{\partial T}{\partial L}=0$ in the denominator (adaptation), or directly take $\tfrac{\partial T}{\partial \mu }\to \infty$, corresponding to infinite kinetic amplification. This could be approximated by kinetic proofreading with a high number N of proofreading steps, where an amplification up to ${\tau }^{N}$ in magnitude can be obtained. For instance, in a model explicitly designed to provide better sensitivity [30], a rather high number of steps (N = 25) is taken to account for specificity compatible with biological ranges, and indeed such models do not yield antagonism. Notice however that if $\tfrac{\partial T}{\partial \mu }\to \infty$, dT in equation (13) varies very rapidly when ${\rm{d}}\mu$ changes. So in this limiting case, the output itself is (by construction) infinitely sensitive to changes of μ. This is why there is no antagonism: subthreshold ligands yield infinitely small contributions compared to the ones at threshold, irrespective of their concentration, and in particular cannot antagonize them. As pointed out in [20], depending on the system considered, it is actually not clear that many proofreading steps are available to the cell. This observation led to the proposal that specificity is rather due to the presence of internal feedbacks, effectively reducing instead $\tfrac{\partial T}{\partial L}$ to perform specific detection, which is the framework of this article.

At least two examples of absolute discrimination/antagonism are offered in the immune context. For early immune recognition mediated by TCRs (the discrimination parameter being the binding time τ), agonists simultaneously presented with ligands just below threshold fail to trigger a response, while agonists alone do [20, 21]. There is some debate in immunology about the parameter defining ligand quality: some argue that the binding time τ defines quality, while other authors [13–15] define an effective parameter (e.g. dwelling time t_a in [14]) to encode ligand strength. However, if all ligands can be hierarchically ordered on one 'quality' axis effectively defining μ, whether $\mu =\tau$ or $\mu ={t}_{a}$, our reasoning applies and antagonism should ensue. For FcRIs mediated response in mast cells [19], where the discrimination parameter is also the binding time τ, a very similar 'dog in the manger' effect is observed where low affinity antagonists titrate the Lyn kinase responsible for proofreading steps (exactly like adaptive sorting evolved in [18]).

Possible instances outside of immunology include antagonism via negative feedback in Hh signalling [36] or hormonal pathways, which present properties very reminiscent of immune recognition [22]. Non-monotonic response activity with varying ligand concentration [26, 37] could correspond in our framework to approximate adaptation where $\tfrac{\partial T}{\partial L}$ is kept small. It has been established that in such contexts, antagonists differ from agonists purely based on slower binding kinetics [38], suggesting a hierarchical ordering of ligand quality as hypothesized here. A recent meta-study shows that most of these pathways indeed use internal negative feedbacks to change the monotonicity of the response [39], thus similar to the negative feedforward contributions predicted from equations (9) and (11).

Is antagonism an evolutionary spandrel related to immune detection? Networks evolved in [18] systematically show local biochemical adaptation and antagonism, as expected from the present derivation. Therefore evolutionary spandrels can be observed and studied in evolutionary simulations, and their emergence studied theoretically. In nature, it is difficult to definitely know if a trait is a spandrel without detailed historical data and comparisons [6]. Other possibilities could be that absolute discrimination has been selected in conjunction with other properties (such as information transmission [40, 41] or statistical decision [42]). But as detailed mathematically in this work, the presence of a phenotypic spandrel here is due to the non-trivial computation performed by the cell to disentangle the kinetics of the binding and ligand concentration. Experiments probing internal feedbacks (such as [19]) connect spandrels to cellular computations, similar to symmetries connected to physical laws.

We thank Jean-Benoît Lalanne and Grégoire Altan-Bonnet for useful discussions. This project has been funded by the Natural Science and Engineering Research Council of Canada (NSERC), Discovery Grant Program, and a Simons Investigator Award for Mathematical Modelling of Biological Systems to PF. LS has been supported by an Undergraduate Student Research Award from NSERC.

Appendix. Equations for adaptive sorting (AS)

$$\begin{eqnarray}&&K={K}_{T}/({C}_{m}+{C}^{* })\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\dot{C}}_{0}=-({\phi }_{0}+{\tau }^{-1}){C}_{0}+\kappa \left(L-\displaystyle \sum _{i}{C}_{i}\right)\left(R-\displaystyle \sum _{i}{C}_{i}\right)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\dot{C}}_{i}=-({\phi }_{i}+{\tau }^{-1}){C}_{i}+{\phi }_{i-1}{C}_{i-1},\hspace{0.15cm}1\leqslant i\end{eqnarray} \tag{ 17 }$$

where ${\phi }_{m}={\phi }_{K}K$, ${\phi }_{N}=0$ and ${\phi }_{i}=\phi$ for other values of i.

The parameters used for simulation in the main text are $\kappa ={10}^{-5}$, $R=3\times {10}^{4}$, $\phi ={\phi }_{K}{K}_{T}=0.09$, ${C}^{* }=3,m=2,N=3$. Output is C_N, i.e. with conventions defined in the main text t = N. We used ${\tau }_{c}=4\ s$ and defined threshold ${\rm{\Theta }}=0.31$ to plot a discrimination line.

We get immediately at steady state

$$\begin{eqnarray}&&{C}_{j}={\tau }^{j}{\lambda }_{j}{C}_{0}\end{eqnarray} \tag{ 18 }$$

setting ${\gamma }_{j}=\tfrac{{\phi }_{j-1}}{{\phi }_{j}\tau +1}$ and ${\lambda }_{0}=1,{\lambda }_{j}={{\rm{\Pi }}}_{1\leqslant k\leqslant j}{\gamma }_{k}$. Neglecting ${\sum }_{i}{C}_{i}$ in front of R we have

$$\begin{eqnarray}&&{C}_{0}=\displaystyle \frac{\kappa {RL}\tau }{\kappa R\tau {\displaystyle \sum }_{j}{\lambda }_{j}{\tau }^{j}+{\phi }_{0}\tau +1}.\end{eqnarray} \tag{ 19 }$$

Assuming ${\phi }_{j}\tau ,\kappa R\tau \ll 1$, we then get the following scaling laws for $i\leqslant m$

$$\begin{eqnarray}&&{C}_{i}=\kappa R{\tau }^{i+1}{\phi }^{i}L\end{eqnarray} \tag{ 20 }$$

and for $i\gt m$

$$\begin{eqnarray}&&{C}_{i}=\kappa R{\tau }^{i+1}{\phi }^{i-1}{\phi }_{K}{KL}.\end{eqnarray} \tag{ 21 }$$

Since $K\propto {C}_{m}^{-1}$, we recover the scaling law of the Ansatz from the main text.

Appendix. Equations for the immune model

$$\begin{eqnarray}&&S={S}_{T}\displaystyle \frac{{C}_{1}}{{C}_{1}+{C}_{* }}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\dot{{C}_{0}} & = & \kappa \left(L-\displaystyle \sum _{i=0}^{N}{C}_{i}\right)\left(R-\displaystyle \sum _{i=0}^{N}{C}_{i}\right)\\ & & +(b+\gamma S){C}_{1}-(\phi +{\tau }^{-1}){C}_{0}\\ {\dot{C}}_{j} & = & \phi {C}_{j-1}+(b+\gamma S){C}_{j+1}\\ & & -(\phi +b+\gamma S+{\tau }^{-1}){C}_{j}\\ {\dot{C}}_{N} & = & \phi {C}_{N-1}-(b+\gamma S+{\tau }^{-1}){C}_{N}.\end{eqnarray} \tag{ 23 }$$

Parameters used for simulation in the main text are $\kappa ={10}^{-4}$, $R=3\times {10}^{4}$, $\phi =0.09$, $b=0.04,\gamma {S}_{T}=0.72$, ${C}_{* }=300$. Output is C_N. We used ${\tau }_{c}=4\ s$ and defined threshold ${\rm{\Theta }}=0.09$ to plot a discrimination line.

Appendix. Origin of linear separation (equation (8))

Considering a family of models similar to the ones defined above (e.g. the ones in [9]), it is clear that if we assume unsaturated receptors, i.e. ${\sum }_{i=0}^{N}{C}_{i}\ll R$, calling ${\bf{C}}$ the vector of occupancies and ${\bf{M}}$ the internal variables , we have:

$$\begin{eqnarray}&&\dot{{\bf{C}}}=\kappa ({\bf{M}})R{\bf{L}}+{ \mathcal T }(R,{\bf{M}}){\bf{C}}\end{eqnarray} \tag{ 24 }$$

${\bf{L}}=(L,0,\ldots ,0)$ is the vector corresponding to ligand input, $\kappa ({\bf{M}})$ the association rate of ligands to receptors—by definition here ligands and receptors bind into the first state C₀. ${ \mathcal T }(R,{\bf{M}})$ is a matrix defining linear rates between occupancy states, depending on internal variables and parameter τ. Dynamics and steady state value of ${\bf{M}}$ are given by extra equations that are model-specific (e.g. equations (15) and (22) above).

For such systems, we have at steady state

$$\begin{eqnarray}&&{\bf{C}}=-\kappa {({\bf{M}})R{ \mathcal T }(R,{\bf{M}})}^{-1}{\bf{L}}\end{eqnarray} \tag{ 25 }$$

from which we can directly compute the c_js as a function of ${\bf{M}}$ to get a functional form similar to equation (8). Then we can use equations defining steady-state values of ${\bf{M}}$ as functions of ${\bf{C}}$ to close the system.

When several ligand types are present, the independence of ligand binding means that a system similar to equation (24) holds for every single vector of occupancy ${{\bf{C}}}^{\tau }=({C}_{j}^{\tau })$ of receptor states bound to ligands with quality τ. Coupling between different types of ligands only occurs via internal variables ${\bf{M}}$. Note that we can also generalize this formalism so that transition rates depend on occupancies (i.e. effectively giving non linear transition rates between states) by assuming that internal variables ${\bf{M}}$ are occupancies themselves. The strong underlying assumption here is that the coupling is global, via total occupancies.

Appendix. Antagonism when adaptation is not perfect

In this section, we briefly illustrate what happens to antagonism when adaptation is not perfect, using the Ansatz presented in the main text as a case-study. As said in the main text, antagonism then relies on a competition between the 'flatness' of T as a function of L and variation of ${\rm{d}}\mu$. To see this, let us consider the Ansatz of the main text with $\beta =-1+\epsilon$, with $\epsilon \gt 0$ so that $T(L,{\mu }_{c})\propto {L}^{\epsilon }$ slowly increases as a function of L for $\mu ={\mu }_{c}$. Writing the equivalent of equation (4) to get:

$$\begin{eqnarray}&&T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})={({L}_{1}+{L}_{2})}^{\epsilon }{ \mathcal B }({\mu }_{c})\\ &&\quad -{L}_{2}{({L}_{1}+{L}_{2})}^{-1+\epsilon }{\rm{d}}\mu { \mathcal B }^{\prime} ({\mu }_{c})\end{eqnarray} \tag{ 26 }$$

calling ${ \mathcal B }={{fg}}^{\beta }$. Then for the difference of output ${\rm{\Delta }}T=T(\{{L}_{1},{\mu }_{c};{L}_{2},{\mu }_{c}-{\rm{d}}\mu \})-T({L}_{1},{\mu }_{c})$ we get

$$\begin{eqnarray}{\rm{\Delta }}T & = & ({({L}_{1}+{L}_{2})}^{\epsilon }-{L}_{1}^{\epsilon }){ \mathcal B }({\mu }_{c})\\ & & -{L}_{2}{({L}_{1}+{L}_{2})}^{-1+\epsilon }{\rm{d}}\mu { \mathcal B }^{\prime} ({\mu }_{c})\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\simeq \epsilon \mathrm{log}\left(1+\displaystyle \frac{{L}_{2}}{{L}_{1}}\right){ \mathcal B }({\mu }_{c})\\ &&-{L}_{2}{({L}_{1}+{L}_{2})}^{-1+\epsilon }{\rm{d}}\mu { \mathcal B }^{\prime} ({\mu }_{c}).\end{eqnarray} \tag{ 28 }$$

One can see that the first term in equation (28) is of order $\epsilon$, while the second term is the usual negative antagonistic term of order ${\rm{d}}\mu$, and in the limit of small L₂ compared to L₁, both terms are linear in L₂. So for ${\rm{d}}\mu$ big enough compared to $\epsilon$, the second term should dominate the first one, and we expect antagonism to occur even without perfect adaptation. This can be observed on an even less general model with the non-monotonic dose response curve from [22], illustrated in figure 3(C), where the output can vary over one decade while the input varies over four decades [22], but still the system displays antagonism as soon as τ is sufficiently below ${\tau }_{c}$.