Introduction

Cells routinely encounter a wide variety of extracellular ligands and decode their identity with remarkable precision to generate context-specific responses. This selective processing of environmental cues is essential for regulating diverse biological processes, including development, immune surveillance, and tissue homeostasis (Cantley et al., 2014). Failures in ligand discrimination underlie many diseases, including diabetes and cancer (Madsen et al., 2025; Madsen and Vanhaesebroeck, 2020).

A key determinant of ligand specificity in biochemical networks is the thermodynamic stability of molecular complexes, such as ligand–receptor or substrate–enzyme pairs. At thermal equilibrium, the abundance of complexes is determined by their equilibrium binding constants. This imposes a fundamental limit on specificity: high-affinity ligands are inevitably favored over lower-affinity competitors, with complex abundances scaling in proportion to their association constants.

However, many signaling networks display paradoxical behaviors that cannot be explained by equilibrium affinity alone (Clark et al., 1999; Coombs et al., 2002; Freed et al., 2017; Madsen et al., 2025; Myers et al., 2023). For example, receptor tyrosine kinases (RTKs) can produce stronger downstream signaling outputs in response to intermediate-affinity ligands than to high-affinity ligands (Freed et al., 2017; Madsen et al., 2025; Myers et al., 2023). Additionally, partially inhibiting kinase activity can paradoxically increase receptor phosphorylation levels (Kiyatkin et al., 2020; Kleiman et al., 2011). These observations raise a fundamental question: how do cells overcome thermodynamic constraints to achieve nuanced, ligand-specific responses?

A classic solution to bypass equilibrium limits is kinetic proofreading (KPR), a mechanism first proposed by Hopfield (Hopfield, 1974) and Ninio (Ninio, 1975). KPR enhances specificity of high affinity ligands by introducing energy-consuming, irreversible steps—such as phosphorylation—that amplify differences between competing ligands. KPR has been invoked in diverse systems, including DNA replication (Hopfield, 1980), mRNA surveillance (Hilleren and Parker, 1999), protein folding (Gulukota and Wolynes, 1994), and immune receptor signaling (McKeithan, 1995; Huang et al., 2019; Lever et al., 2014). Yet, most KPR models (with a few exceptions (Lever et al., 2014; Murugan et al., 2014), see below) assume that correct recognition correlates with the highest affinity - a premise that fails in systems where low affinity ligands dominate signaling outputs compared to high affinity ligands.

In this study, we investigate how signaling networks achieve ligand specificity through non-equilibrium mechanisms that go beyond classical KPR. We develop a simple, biologically grounded model combining two ubiquitous features of cellular signaling: sequential multi-site phosphorylation and ligand-induced receptor degradation. These two motifs are found in many major receptor systems, including RTKs (Furdui et al., 2006; Sorkin and Goh, 2009), G protein-coupled receptors (GPCRs)(Koenig and Edwardson, 1997; Tobin, 2008), T cell receptors (McKeithan, 1995; Charpentier and King, 2021), and interleukin receptors (Kollewe et al., 2004; Cendrowski et al., 2016).

Our model shows that high-affinity ligand-receptor complexes are preferentially sorted toward degradation-prone states, while low-affinity complexes repeatedly dissociate the ligand, resulting in maximal signaling output only from intermediate-affinity ligands. Notably, this ligand specificity can be tuned by tuning easily controllable cellular parameters, e.g. enzyme abundances. This non-equilibrium kinetic sorting mechanism explains the paradoxical non-monotonic dependence of signaling activity on ligand affinity and phosphorylation rate observed in RTKs. More broadly, we propose that kinetic sorting provides a general strategy for achieving ligand discrimination across diverse natural and synthetic signaling networks.

Results

Classic kinetic proofreading always favors high-affinity ligands

Kinetic proofreading (KPR) is the standard model for non-equilibrium ligand discrimination. To set the stage, we first revisited the classic KPR model originally proposed by McKeithan to explain how T cell receptors avoid activation downstream of weak ligands (McKeithan, 1995) (Fig. 1a; see Supplementary Materials for equations).

Figure 1.

In this model, ligand-bound receptors undergo a series of phosphorylation steps, with the final state P_N representing the active, signaling-competent form. Importantly, ligand unbinding at any phosphorylation stage returns the receptor to the unbound state R. We parameterized the model using dimensionless quantities: the ligand dissociation rate δ = k_dτ, phosphorylation rate ω = k_pτ, and ligand concentration u = L/K_D, where K_D = k_d/k_on. Assuming saturating ligand (u → ∞), the steady-state abundance of the active state is:

As expected, increasing the phosphorylation cascade length N amplifies the preference for low-dissociation (high-affinity) ligands (Fig. 2a), reflecting the classical KPR outcome.

Figure 2.

Modified KPR schemes cannot explain paradoxical RTK behavior

Next, we examined two previously proposed extensions of KPR: the sustained signaling model and the limited signaling model (Lever et al., 2014) (Fig. 1b,c). Both models introduce an additional state to Mckeithan’s KPR scheme. The sustained signaling model adds an active but ligand-free state , while the limited signaling model introduces an inactivated state I downstream of P_N.

Of these, only the limited signaling model exhibits non-monotonic dependence on ligand dissociation rates at saturating ligand concentrations (Lever et al., 2014) (Fig. 2b), consistent with some paradoxical features observed in RTKs (Freed et al., 2017; Madsen et al., 2025; Myers et al., 2023). However, it fails to reproduce a second key observation: phosphorylation in this model increases monotonically with kinase activity, whereas RTK experiments show that partial kinase inhibition can paradoxically increase phosphorylation (Kiyatkin et al., 2020; Kleiman et al., 2011) (Fig. 2c). Thus, these models are insufficient to explain RTK signaling dynamics.

A kinetic sorting model integrates ligand-induced degradation

To address these gaps, we present a model incorporating two widespread signaling motifs: sequential multi-site phosphorylation and ligand-induced receptor degradation. In our model, receptors are delivered to the surface at a constant rate, internalized at a basal rate k_int, and degraded more rapidly when highly phosphorylated . Ligand-bound receptors undergo reversible phosphorylation. A key feature of our model is that all phosphorylated species are signaling competent. Indeed, in many signaling pathways all phosphorylation sites on the receptor (Schulze et al., 2005; Tobin, 2008; Kollewe et al., 2004; Lemmon and Schlessinger, 2010; Latorraca et al., 2020) have downstream effects. Therefore, we define the net activity A_n of phosphorylation site n as all receptor states where the site n is phosphorylated: A_n =∑_m≥n P_m.

Parameter ranges

To ensure that the phenomena captured by our model are relevant to real signaling networks, we selected ranges for the dimensionless parameters based on direct experimental measurements and model fits. Importantly, many of these kinetic processes have comparable rates across diverse receptor systems (Koenig and Edwardson, 1997; Subtil et al., 1994; Liu et al., 2000). Specifically, basal receptor internalization occurs at rates of k_int ≈ 10⁻⁴–10⁻³, s⁻¹(Wiley, 2003), while ligand-induced internalization is faster, at (Wiley, 2003; Lyashenko et al., 2020). Ligand dissociation rates typically fall in the range k_d ≈ 10⁻²–10⁻¹, s⁻¹(Chen et al., 2009; Lyashenko et al., 2020), and receptor phosphorylation and dephosphorylation occur at ∼ 10⁻¹–10⁰, s⁻¹(Kleiman et al., 2011; Chen et al., 2009; Lyashenko et al., 2020). For EGFR, equilibrium dissociation constants range from ∼ 0.1, nM for the high-affinity ligand Betacellulin to ∼ 25, nM for the low-affinity ligand AREG (Hu et al., 2022; Macdonald-Obermann and Pike, 2014). Based on these values, we set the following ranges for our dimensionless parameters: β ∈ [1, 100], ρ ∈ [0.01, 100], ω ∈ [1, 1000], and δ ∈ [1, 1000]. Finally, the number of phosphorylation sites with known functional roles typically ranges from 5 to 25 (Schulze et al., 2005). These broad ranges comfortably encompass experimentally measured estimates. Unless otherwise specified, our default parameter values are δ = 20, ω = 200, ρ = 0.01, β = 50, and N = 10.

Before examining how phosphorylation levels depend on model parameters, we illustrate the mechanism of kinetic sorting of receptor states, which tunes ligand specificity beyond pure thermodynamic preference, using a simple example. To that end, we consider a signaling network with N = 5 phosphorylation sites interacting with three ligands of distinct affinities—high, medium, and low. We assume the dissociation rates for these ligands are δ_H = 20, δ_M = 200, and δ_L = 1000, respectively. In order to compare our model with aforementioned paradoxical experimental observations which have been performed at saturating ligand concentration, we take the limit u → ∞.

Fig. 3 shows that low affinity ligands (δ_L = 1000) predominantly sort receptors towards the inactive state B and early phosphorylation states P_n, n ∼ 1 as frequent ligand unbinding prevents progression to later phosphorylation states. This behavior resembles the traditional KPR mechanism described by McKeithan (McKeithan, 1995). In contrast, receptors bound to high affinity ligands are sorted toward later phosphorylation states, which mark them for enhanced degradation.

Figure 3.

Here, similar to traditional KPR, the fraction of receptors reaching the final phosphorylation state is highest. Yet, the overall receptor pool is reduced due to ligand-induced degradation, lowering net phosphorylation activity. Strikingly, receptors bound to intermediate affinity ligands (δ_M = 200) are sorted towards intermediate phosphorylation states, resulting in maximal phosphorylation output. Below, we kinetic parameters govern the ability of the network to overcome thermodynamic preference and acquire ligand specificity.

Early phosphorylation sites show ligand-specificity

Figure 4a illustrates how total phosphorylation activity at each site, A_n, n ∈ [1, N] varies with ligand dissociation rate δ. Notably, early phosphorylation sites (n ∼ 1) exhibit maximal activity at intermediate values of δ while both high- and low-affinity ligands suppress net receptor phosphorylation. Our model predicts that this ligandspecificity diminishes for later sites, where outputs increasingly resemble traditional KPR which favors high-affinity ligands.

Figure 4.

To examine how model parameters shape ligand specificity, we focused on the activity at the first phosphorylation site, A₁, which exhibits the strongest discriminatory behavior (Fig. 4a). As shown in Fig. 4b, achieving ligand specificity at high dissociation rates δ requires sufficiently high phosphorylation rates ω. Notably, our model captures a puzzling observation from EGFR signaling: the high-affinity ligand EGF produces lower/comparable steady-state phosphorylation compared to lower-affinity ligands such as Epigen and Epiregulin (Freed et al., 2017; Myers et al., 2023; Madsen et al., 2025). Experimental estimates place the basal EGFR internalization rate at k_int ≈ 1.3 × 10⁻³, s⁻¹ (Chen et al., 2009), the EGF dissociation rate at k_d ≈ 3 × 10⁻², s⁻¹ (Chen et al., 2009), and the phosphorylation rate at k_p ≈ 10⁻¹ −10⁰, s⁻¹, yielding δ_EGF ≈ 10–20 and ω_EGFR ≈ 100–1000. Low-affinity ligands such as Epigen (EPGN) and Epiregulin (EREG) have equilibrium dissociation constants about 10-fold higher than EGF (Hu et al., 2022), corresponding to δ_EPGN ≈ δ_EREG ≈ 100–200. The effective degradation rate of fully activated receptors is estimated to be 10–50 times higher than that of inactive receptors (Lyashenko et al., 2020), implying β = 50. Under these conditions, our model predicts a switch in phosphorylation levels: as δ increases from δ_EGF to δ_EPGN, receptor phosphorylation increases—reversing the expectation based purely on thermodynamic affinity. This effect arises because EGF-bound receptors are efficiently sorted towards degradation-prone states compared to those bound to lower-affinity ligands.

Our model also explains another paradox in EGFR signaling. Experimental studies have shown that EGF-stimulated receptors exhibit higher steady-state phosphorylation when kinase activity is partially inhibited (Kiyatkin et al., 2020; Kleiman et al., 2011). As shown in Fig. 4c, at low δ values (e.g., δ = 16), decreasing the phosphorylation rate ω from levels typical of EGFR (ω_EGFR ≈ 100–1000) paradoxically increases overall receptor phosphorylation. A similar effect is observed when receptor dephosphorylation is enhanced (Fig. 4d). Importantly, our model makes a testable prediction: the reversal of thermodynamic preference observed between EGF and EPGN/EREG will disappear when kinase activity is mildly suppressed (see, e.g., the curves for ω = 256 and ω = 16 over δ ∈ [10, 100]), such as by treatment with low doses of the kinase inhibitor gefitinib (Herbst et al., 2004).

Multi-site phosphorylation and ligand-induced degradation are both essential for ligandspecificity

To assess the importance of sequential multi-site phosphorylation on ligand specificity, we analyzed , the phosphorylation of the first site for signaling networks with N phosphorylation sites. Fig. 5a shows that multi-site phosphorylation is essential to endow signaling networks with ligand specificity and ligand-induced receptor degradation alone is not sufficient. This is because the non-monotonic preference for intermediate affinity ligands arises only when the receptors can be sorted among multiple phosphorylation sites: earlier ones for low affinity ligands and later ones for high affinity ligands.

Figure 5.

To assess how receptor degradation shapes ligand specificity for a multi-site phosphorylation network, we examined how altering receptor turnover influences model behavior. As shown in Fig. 5b, the optimal dissociation rate δ_opt —which maximizes receptor phosphorylation levels—increases with ligand-induced degradation rate β. When kinase activity ω, low affinity ligands (high δ) over high affinity ligands. Crucially, this optimal δ_opt emerges only when receptor degradation is strong (β ≫ 1). These predictions can be tested by blocking receptor degradation, e.g., via mutation of ubiquitination sites (Gerritsen et al., 2023).

To quantify ligand specificity, we computed receptor phosphorylation in response to ligands differing by at least one k_BT in binding free energy from the optimal ligand. Figure 5c shows that as β increases, phosphorylation downstream of suboptimal ligands (red line in inset) declines relative to the optimal ligand. This enhanced specificity is further amplified by increasing kinase activity ω.

These results show that both multi-site phosphorylation and ligand induced degradation are essential for ligand specificity.

Discussion

Cells face the formidable task of decoding multiple extracellular signals to generate appropriate, context-specific responses. This challenge is especially acute for receptors like RTKs, GPCRs, and interleukin receptors, which bind multiple cognate ligands and yet elicit distinct downstream outcomes. While equilibrium affinity provides a baseline for ligand specificity, it cannot fully explain the rich and often counterintuitive behaviors observed in many signaling systems.

Here, we show that non-equilibrium mechanisms—specifically, kinetic sorting through multi-site phosphorylation and ligand-induced degradation—can explain how signaling networks achieve lig- and specificity beyond equilibrium limits. In the model, high-affinity ligand–receptor complexes are sorted toward degradation-prone states, low-affinity complexes are sorted towards inactivated states, and intermediate-affinity ligands strike the optimal balance between progression and degradation to maximize signaling. This framework explains paradoxical features observed in RTK systems, including the non-monotonic dependence of phosphorylation on ligand affinity and kinase activity.

Importantly, our model predicts that early phosphorylation sites show the strongest ligand discrimination, consistent with recent experimental observations. It also makes the testable prediction that impairing receptor degradation should reduce specificity by eliminating the kinetic sorting effect.

More broadly, our findings suggest that ligand-specific sorting of multidimensional receptor states—across phosphorylation, degradation, localization, and engagement fates—may be a general strategy for encoding ligand identity. While our model focused on the activity–degradation axis, real-world receptors operate in even richer state spaces, offering exciting directions for future work.

Our findings complement prior studies on mechanisms of ligand specificity that operate at thermal equilibrium, such as those described in the Bone Morphogenetic Protein (BMP) pathway (Antebi et al., 2017; Su et al., 2022; Parres-Gold et al., 2025). BMP signaling relies on promiscuous ligand–receptor interactions, with specificity emerging from differences in receptor abundance, binding affinity, and complex activity. By contrast, our work shows that non-equilibrium mechanisms—such as phosphorylation cycles and ligand-induced receptor degradation—can achieve ligand discrimination even for a single receptor type. Given that ligand–receptor promiscuity, multisite phosphorylation, and receptor turnover are common features across signaling systems (e.g., in the EGFR/ErbB family (Linggi and Carpenter, 2006)), it is likely that biological networks integrate both equilibrium and non-equilibrium strategies to achieve robust and tunable ligand specificity.

In recent years, there has been growing interest in engineering synthetic physical and chemical circuits capable of carrying out complex computational tasks, including input discrimination, classification, prediction, and the generation of multiple stable cell states (Shakiba et al., 2021; Ma et al., 2022; Benzinger et al., 2022; Zhu et al., 2022; Floyd et al., 2024; Parres-Gold et al., 2025; Aoki et al., 2019). Some of these synthetic strategies rely on equilibrium thermodynamics (Parres-Gold et al., 2025), while others exploit non-equilibrium steady states (Floyd et al., 2024). We propose that non-equilibrium kinetic sorting, which harnesses receptor synthesis and degradation, could provide synthetic biologists with a powerful framework for achieving precise control over molecular abundances and dynamic system behavior.

Finally, we address a major concern in non-equilibrium signaling circuits: the energetic cost of operation. Previous theoretical work has shown that free energy dissipation places fundamental constraints on the performance of signaling networks (Bryant and Machta, 2023; Lan et al., 2012; Mehta and Schwab, 2012; Qian and Reluga, 2005; Cao et al., 2015; Azeloglu and Iyengar, 2015; Floyd et al., 2024; Mahdavi et al., 2024). These studies typically focus on futile cycles of reversible modifications such as phosphorylation or methylation. In contrast, ligand-induced receptor degradation—a central feature of many signaling networks—is a far more energy-intensive process. For example, MCF10A cells maintain approximately 10⁵ EGFR molecules on the surface (each 1,210 amino acids in length)(Shi et al., 2016), with a synthesis rate of about 15 receptors per second (Lyashenko et al., 2020), corresponding to an energetic cost of roughly ∼ 8 × 10⁴ ATP/sec (assuming 4.5 ATP per peptide bond (Milo et al., 2010)). By comparison, EGFR dephosphorylation occurs over ∼ 15 seconds (Kleiman et al., 2011), and only 5−10% of receptors are phosphorylated at steady state (Shi et al., 2016; Feng et al., 2023), resulting in a much lower energetic cost of ∼ 6 × 10² ATP/sec for dephosphorylation. Thus, the energetic burden of receptor turnover can exceed that of reversible modification cycles by up to two orders of magnitude. These estimates suggest that, at least in eukaryotic systems, the energetic demands of non-equilibrium modification cycles are unlikely to pose a fundamental limitation on the functionality of signaling networks.

Supplementary materials

Equations for proofreading models

The equations describing species abundances in the traditional KPR model similar to that of McKeithan (McKeithan, 1995) are as follows:

For the limited signaling model, the dynamics of B, and P_i, i ∈ [1, N − 1] are identical to the traditional KPR model. The dynamics of R and P_N are modified as follows:

Equations for the model with receptor degradation

Signaling receptors participate in a variety of complex regulatory processes, including non-linear ligand binding dynamics (Limbird et al., 1975; Macdonald and Pike, 2008), receptor oligomerization (Mudumbi et al., 2024; Huang et al., 2016), context-specific interactions with adapter proteins (Madsen and Vanhaesebroeck, 2020; Feng et al., 2023), and trafficking between cellular compartments leading to degradation (Sorkin and Goh, 2009; Wiley, 2003; Irannejad and Von Zastrow, 2014).

While computational models that incorporate these mechanistic details are powerful tools for hypothesis generation (Chen et al., 2010; Qiao et al., 2025), they often require large-scale datasets for accurate parameterization (Feng et al., 2023). As an alternative, simplified models that intentionally omit certain mechanistic details can still yield deep qualitative insights, even if they cannot quantitatively reproduce experimental data.

In this study, we present such a simplified model aimed at explaining two paradoxical features of receptor tyrosine kinase (RTK) signaling: (1) the non-monotonic relationship between ligand-receptor affinity and steady-state receptor phosphorylation (Freed et al., 2017; Madsen et al., 2025; Myers et al., 2023), and (2) the counterintuitive increase in receptor phosphorylation following mild kinase inhibition (Kleiman et al., 2011; Kiyatkin et al., 2020).

To keep the model simple and tractable, we neglect receptor recycling and oligomerization. Previously, we showed that the combined effects of endocytosis, recycling, and degradation can be captured by a single effective dimensionless parameter, β in this study, which reflects the degradation bias of fully phosphorylated receptors compared to partially phosphorylated receptors(Lyashenko et al., 2020). Similarly, receptor dimerization and negative cooperativity can be abstracted into a Hill coefficient η < 1 (Lyashenko et al., 2020). For the phenomena explored here, including oligomerization would modify the shape of the response curves but not their qualitative behavior.

Under these assumptions, the governing equations for the model are given by:

All equations are solved at steady state and in the limit u → ∞. All codes required to generate the figures in the manuscript can be found at https://github.com/BarriosJer0/KineticSorting.

Additional information

Funding

National Institutes of Health (R35GM142547)