Critical lack of equilibrium in stochastic kinetic proofreading

TL;DR

This paper investigates how the absence of detailed balance in stochastic kinetic proofreading models influences their ability to discriminate between ligands, revealing a critical threshold necessary for high specificity.

Contribution

It establishes a critical level of detailed balance violation required for strong ligand discrimination and shows that exceeding this threshold does not guarantee improved specificity.

Findings

Existence of a critical lack of detailed balance for high specificity

Models with similar lack of detailed balance can differ in discrimination ability

Strong discrimination is not solely determined by the degree of detailed balance violation

Abstract

In this paper we study the relation between the property of detailed balance and the ability of discriminating between different ligands for a class of stochastic models of kinetic proofreading. We prove the existence of a critical amount of lack of detailed balance that the kinetic proofreading models must have in order to have strong specificity for a value of the binding energy $\sigma$. We also prove that the fact that a kinetic proofreading model has a lack of detailed balance that is larger than the critical one does not necessarily yield strong discrimination properties. Indeed, there exist different sets of chemical rates, leading to the same amount of lack of detailed balance, that have strong discrimination property in some cases and not in others.

Figures

Equations

$$\prod _{R\in \omega }\frac{K_R}{K_{-R}}=1,$$

$$S\stackrel{R_1^{(k)}}{⟶}{C}_k\stackrel{R_2^{(k)}}{⟶}{C}_{k+1}\stackrel{R_3^{(k)}}{⟶}S.$$

$$e^{{\Delta }_k}:=\frac{K_{R_1^{(k)}}{K}_{R_2^{(k)}}{K}_{R_3^{(k)}}}{K_{-R_1^{(k)}}{K}_{-R_2^{(k)}}{K}_{-R_3^{(k)}}}.$$

$$T_N\sim \exp \left({\lambda }_{\Delta }(\sigma ,E)N\right)\text{ as }N\to \infty ,$$

$$T_{life}\sim \exp \left(bN\right)\text{ as }N\to \infty ,$$

$$p_{res}(\sigma )\sim \frac{e^{-{\lambda }_{\Delta }(\sigma ,E)N}}{e^{-{\lambda }_{\Delta }(\sigma ,E)N}+e^{-bN}}\text{ as }N\to \infty$$

$$|\sigma -{\sigma }_c|\lesssim \frac{1}{N}.$$

$$R_1^{(k)}:S\to C_k$$

$$-R_1^{(k)}:C_k\to S$$

$$R_2^{(k)}:C_k\to C_{k+1}$$

$$-R_2^{(k)}:C_k\to C_{k-1}$$

$$\mu \sim e^{-bN}\text{ as }N\to \infty$$

$$S\stackrel{R_1^{(k)}}{⟶}{C}_k\stackrel{R_2^{(k)}}{⟶}{C}_{k+1}\stackrel{-R_1^{(k+1)}}{⟶}S$$

$$e^{\Delta }=\frac{K_{R_1^{(k)}}{K}_{R_2^{(k)}}{K}_{-R_1^{(k+1)}}}{K_{-R_1^{(k)}}{K}_{-R_2^{(k)}}{K}_{R_1^{(+1k)}}}.$$

$$\frac{d{n}_S}{dt}$$

$$\frac{d{n}_0}{dt}$$

$$\frac{d{n}_k}{dt}$$

$$\frac{d{n}_N}{dt}$$

$$S_N(E):=\sum _{k=1}^N{e}^{-kE}=\frac{1-e^{-NE}}{e^E-1}.$$

$$R(t):=\alpha e^{\Delta }{n}_N(t).$$

$$p_{deg}(t):=-\mu M(t)$$

$$M(t):=n_S(t)+\sum _{k=0}^N{n}_k(t).$$

$$\frac{dM}{dt}=-\alpha e^{\Delta }{n}_N-\mu M.$$

$$p_{deg}(t)+R(t)+M(t)=1,\forall t\ge 0.$$

$$p_{res}(\sigma ):={\int }_0^{\infty }R(t)dt=\alpha e^{\Delta }{\int }_0^{\infty }{n}_N(t)dt$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\log \left(\frac{1}{p_{res}(\sigma )}-1\right)={\lambda }_{\Delta }(\sigma ,E)-b$$

$${\lambda }_{\Delta }({\sigma }_c,E)=b.$$

$${\lambda }_{\Delta }(\sigma ,E)=\{\begin{array}{ll}E& \text{ if }\Delta <{\Delta }_c\\ {\psi }_{\Delta }(\sigma ,E)& \text{ if }\Delta >{\Delta }_c\end{array}$$

$${\psi }_{\Delta }(\sigma ,E):=E-\log \left[\frac{1}{2\alpha }\left(e^{\sigma }+\alpha (e^E+e^{\Delta })-\sqrt{(e^{\sigma }+\alpha (e^E+e^{\Delta }){)}^2-4{\alpha }^2{e}^{\Delta +E}}\right)\right]$$

$${\Delta }_c={\Delta }_c(\sigma ):=\log \left(1+\frac{e^{\sigma }}{\alpha (e^E-1)}\right).$$

$$p_{res}(\sigma )\sim {\left(1+C{e}^{N(\lambda (\sigma ,E)-b)}\right)}^{-1}\text{ as }N\to \infty$$

$$|\lambda (\sigma ,E)-\lambda ({\sigma }_c,E)|=|\lambda (\sigma ,E)-b|\lesssim \frac{1}{N}.$$

Full text

Keywords: detailed balance, criticality, kinetic proofreading, matched asymptotic expansions, Laplace transforms.

1. 1 Introduction

2. 1.1 Notation

3. 2 Main result of the paper

4. 2.1 The model

5. 2.2 Strong discrimination properties for $\Delta>\Delta_{c}$

6. 2.3 Necessary conditions on $b$ for strong discrimination when $\Delta>\Delta_{c}$

7. 3 Kinetic proofreading models out equilibrium

8. 3.1 Kinetic proofreading model involving molecules of ATP, ADP and phosphate groups

9. 3.2 Lack of detailed balance for constant non equilibrium concentration of ATP

10. 4 Matched asymptotic expansions

11. 4.1 Outer region: $N-k\gg 1$ as $N\to\infty$.

12. 4.1.1 Asymptotic behaviour of $n_{k}^{(out)}$ when $\Delta<\Delta_{c}$ and $k\gg 1$, $N-k\gg 1$.

13. 4.1.2 Asymptotic behaviour of $n_{k}^{(out)}$ when $\Delta>\Delta_{c}$ and $k\gg 1$, $N-k\gg 1$.

14. 4.1.3 Asymptotic behaviour in the critical case $\Delta=\Delta_{c}$ when $k\gg 1$ and $N-k\gg 1$

15. 4.2 Inner region: $N-k\approx 1$ as $N\to\infty$

16. 4.2.1 Subcritical case $\Delta<\Delta_{c}$.

17. 4.2.2 Supercritical case $\Delta>\Delta_{c}$.

18. 4.2.3 Critical case $\Delta=\Delta_{c}$.

19. 4.3 Discrimination property

20. 5 Approximated model on $\mathbb{N}$

21. 6 Analysis of the model with a finite number of kinetic proofreading steps

22. 7 Some PDE limits

23. 7.1 PDE limit when $E\approx 1$, $\Delta\approx 1$, $\Delta>E$, $e^{\sigma}\approx\frac{1}{N}$ and $\mu\approx\frac{1}{N}$

24. 7.2 PDE limit when $E\approx\frac{1}{N}$, $\Delta\approx 1$, $e^{\sigma}\approx\frac{1}{N}$ and $\mu\approx\frac{1}{N}$

25. 8 Some models without detailed balance and lack of strong discrimination properties

26. 8.1 Detachment rate depending on $\Delta$

27. 8.2 Attachment rate depending on $\Delta$

28. 8.3 Dephosphorylation rate depending on $\Delta$

29. 8.4 $\Delta\to\infty$ and $\alpha e^{\Delta}\approx 1$

30. 9 Conclusions and open problems

1 Introduction

Kinetic proofreading models are chemical reaction networks that can distinguish between ligands leading to a correct output and ligands leading to a wrong output. A family of kinetic proofreading models were introduced by Hopfield and Ninio in the 70’s (see [14, 19]). Generalizations of this model, including inhibition effects, feedbacks and rebindings, have been studied by several authors (see for instance [5, 7, 11, 12, 13, 24, 25]).

It was noticed in [14] that the kinetic proofreading networks must work in out of equilibrium conditions in order to perform their function. The kinetic proofreading model analysed in [14] is a deterministic linear chemical reaction network modelled by a system of ODEs. Computing the flux solutions of this system, it is possible to observe that the specificity to detect ligands can be much higher in systems that are out of equilibrium compared to the optimal specificity that can be achieved in equilibrium conditions. The relation between the dissipation of free energy and the specificity properties of ODEs models of kinetic proofreading systems has been analysed in several papers in the biophysics literature (for instance in [3, 22, 27, 29]).

In this paper we study a probabilistic version of the classical deterministic linear model of kinetic proofreading proposed by Hopfield. The ligands are assumed to be characterized by their binding energy $\sigma$ with the receptor. We assume that when ligands bind to a receptor they form a complex that can be at different phosphorylation states. In particular we assume that $N>1$ is the number of states. Once the complex is formed the ligand can detach from the receptor, with a rate that depends on its binding energy $\sigma$ with the receptor, or it can either gain a phosphate group (i.e. a phosphorylation event takes place), or it loses an inorganic phosphate group (i.e. a dephosphorylation event takes place). If the complex ligand receptor reaches the phosphorylation state $N$, then it produces an output. The chemical reaction network that we consider can be summarized as follows

{\ \ \ \ \ S\ \ \ \ \ }$${C_{0}}$${C_{1}}$${\ldots\ldots}$${C_{N-1}}$${C_{N}}$${output}

here $S$ represents the ligands and $\{C_{k}\}_{k=0}^{N}$ are the complexes ligand receptor at different phosphorylation states.

In this paper we are interested in applications of the kinetic proofreading mechanism to recognition processes of the immune system. There is experimental evidence that low densities of foreign antigens are able to trigger an immune response (e.g. [1]). Therefore, in this paper we consider a stochastic version of the classical deterministic kinetic proofreading model that allows to study the probability that one ligand induces a response of the immune system. Other stochastic models of kinetic proofreading have been studied in [6, 16, 17, 18, 28]. However, the question that we address in this paper, i.e. whether strong specificity can take place in equilibrium systems (or in systems that are close to equilibrium in a suitable sense that will be clarified later) is not studied in these papers.

In the context of immune recognition processes, a ligand is a complex, pMHC, made of a peptide attached on a major histocompatibility complex (MHC). The receptors are the receptors of the $T$-cells and the output corresponds to the response of the immune system (we refer to [15] for details on the role of kinetic proofreading mechanisms in the immune system). The goal of the kinetic proofreading mechanism is to discriminate between self-ligands, that should not induce an immune response, and foreign-ligands, corresponding for instance to pathogens, that are expected to produce an immune response.

The chemical reaction network that we consider in this paper is a generalization of the model that we analysed in [10]. In that paper we assume that most of the chemical reactions are one-directional. In particular we assume that dephosphorylation does not occur and that the attachment of a ligand with a receptor always forms a dephosphorylated complex. In the model that we study in this paper, instead, all the reactions, except for the reaction which produces the output, are assumed to be bidirectional. This generalization allows us to study the role of the lack of detailed balance for the linear kinetic proofreading networks proposed by Hopfield. In Subsection 8.4 we will also show that the model studied in [10] can be obtained as a limit case of the models studied in this paper.

The property of detailed balance in chemical reaction networks states that at steady state every chemical reaction is balanced by its reverse reaction. In particular this implies that, at steady state, there are no fluxes of chemicals in the system. We stress that the property of detailed balance must hold in every chemical reaction network at constant temperature that does not exchange substances with the environment. However, many biological systems, exchange energy with the environment and can be modelled by effective systems for which the property of detailed balance fails (see [9]).

The model of kinetic proofreading that we study in this paper does not satisfy the detailed balance property as it is an effective model that describes a chemical network that is in contact with reservoirs of ATP, ADP and phosphate groups that we denote with $P$, which are out of equilibrium. In Section 3 we introduce a model of kinetic proofreading in which we explicitly take into account of the role of ATP, ADP and phosphate groups in the chemical reactions. This model satisfies the property of detailed balance. We then show that if we assume that the concentration of ATP, ADP and phosphate groups is kept at constant non equilibrium values by means of an active mechanism, then we recover the kinetic proofreading model without detailed balance that we study in this paper. Notice that here we do not try to model the active mechanisms that keep the concentration of ATP, ADP and phosphate groups out of equilibrium.

In the previous paper [10] we studied a class of kinetic proofreading models, simpler than the ones considered in this paper, for which we prove that they have strong discrimination properties (strong specificity). In this paper we say that a kinetic proofreading network has strong discrimination properties if it is able to distinguish with low error rate between ligands with differences of binding energies that are of order $1/N$, where we recall that $N$ is the number of kinetic proofreading steps (i.e. the number of phosphorylation states). One of the goals of this paper is to prove that a minimal amount of lack of detailed balance is necessary in order to obtain this strong discrimination property. We consider a class of linear kinetic proofreading networks of the same type as the network introduced in [14]. We prove that there exists a critical lack of detailed balance that the system must have in order to have strong specificity.

We now explain how we measure in a quantitative manner the lack of detailed balance. It is well known that the fact that the detailed balance property holds or not depends on the chemical rates of the reactions that belong to the cycles of the network. Indeed the Wegscheider criterium (see [27]) guarantees that the property of detailed balance holds in a chemical reaction network if and only if

[TABLE]

for every cycle $\omega$ in the chemical network. Here a cycle is a set of chemical reactions that when applied to a vector of concentrations has a null net effect. In (1.1) we are denoting with $K_{R}$ the rate of the reaction $R$ and with $K_{-R}$ the rate of $-R$, the reverse of the reaction $R$.

Notice that the set of the cycles of a chemical reaction network has the structure of Abelian group. For the specific kinetic proofreading network that we study in this paper it is possible to define a basis of cycles that contains three reactions. Indeed we consider the basis of cycles $\{\omega_{k}\}_{k=1}^{N}$, where for every $k$ we have that the cycle $\omega_{k}=[R_{1}^{(k)},R_{2}^{(k)},R_{3}^{(k)}]$ has the form

[TABLE]

A convenient way to measure the lack of detailed balance of the kinetic proofreading network is to introduce the parameter $\Delta_{k}$, associated to each cycle $\omega_{k}$, defined as

[TABLE]

Notice that if $\Delta_{k}=0$ for every $k=1,\dots,N$, then the detailed balance property holds. We assume that every cycle has the same lack of detailed balance, i.e. for every $k$ it holds that $\Delta_{k}=\Delta$. We will use the parameter $\Delta$ as a measure of the lack of detailed balance of the kinetic proofreading model.

We prove that the class of kinetic proofreading models studied in this paper can exhibit strong discrimination properties only if $\Delta>\Delta_{c}$ where $\Delta_{c}$ depends on the parameters of the model, i.e. on the phosphorylation rates, the energies of the substances in the network and the binding energy $\sigma$ of the ligand with the receptor. Hence, the fact that the chemical rates satisfy the condition $\Delta>\Delta_{c}$ is a necessary condition in order to have strong discrimination properties. Since the value of $\Delta$ measures the lack of detailed balance (lack of equilibrium) one could naively think that this value could characterize the behaviour of the system. In other words one could think that the condition $\Delta>\Delta_{c}$ is sufficient to have strong discrimination properties. It turns out that this is not the case. In Section 8 we show that some choices of chemical rates $\{K_{R_{j}}\}$ yielding lack of detailed balance, and specifically satisfying $\Delta>\Delta_{c}$, produce strong discrimination properties while other choices of chemical rates, that give the same value of $\Delta$, do not. Therefore, the fact that the network does not satisfy the detailed balance property, is not a sufficient condition for strong specificity, even if $\Delta$ is very large. Indeed, it turns out that the details of the chemical reactions rates are also crucial.

We already anticipated that we say that a chemical reaction network has strong discrimination properties if it can distinguish with low error rates ligands that are characterized by binding energies whose difference is of order $1/N$. We clarify now what we mean with that in this paper. We prove that the characteristic time at which a complex ligand-receptor reaches the state $C_{N}$ and therefore produces a response scales like

[TABLE]

where $\lambda_{\Delta}(\sigma,E)$ is a function of $\sigma,E$ and $\Delta$. Suppose that the ligands have a characteristic life time that satisfies

[TABLE]

for a positive $b$. Then, the probability $p_{res}(\sigma)$ that one ligand produces a response satisfies

[TABLE]

and will be close to $1$ when $b-\lambda_{\Delta}(\sigma,E)\gtrsim 1/N$ and close to [math] when $\lambda_{\Delta}(\sigma,E)-b\gtrsim 1/N$.

In this paper we prove that for any binding energy $\sigma_{c}\in\mathbb{R}$ there exists a critical amount of detailed balance $\Delta_{c}(\sigma_{c})$ that has the following property. If $\Delta>\Delta_{c}(\sigma_{c})$ then we can find a parameter $b\in(0,E)$ such that $\lambda_{\Delta}(\sigma_{c},E)=b<E.$ Around the parameter $\sigma_{c}$ we have an abrupt transition from response to non-response. Indeed, if $\sigma-\sigma_{c}\gtrsim\frac{1}{N}$ we have that $p_{res}(\sigma)\approx 0$ while if $\sigma_{c}-\sigma\gtrsim\frac{1}{N}$, then $p_{res}(\sigma)\approx 1$. The region of transition from response to non-response, where we have that $p_{res}(\sigma)\in(0,1)$, is of order $1/N$ and is the region where we have

[TABLE]

Hence the accuracy of the discrimination is of order $1/N$. If instead we have that $\Delta<\Delta_{c}(\sigma_{c})$ it is not possible to find a $b$ such that $\lambda_{\Delta}(\sigma_{c},E)=b$. Hence, in this case, we do not have a sharp transition from non response to response around a critical binding energy separating the energies leading to a response from the energies that do not lead to a response. In this case the network does not have strong discrimination properties.

Notice that in this paper we assume that $T_{life}$ satisfies (1.3). However, also if $T_{life}\lesssim N$ it would be possible to study the discrimination properties of the kinetic proofreading model. Also in this case we would have some discrimination properties, meaning that higher $\sigma$ implies a smaller probability of response. However in that setting we would not obtain the existence of a sharp transition from $1$ to [math] for the probability of response around a critical value of $\sigma$.

An interpretation of the critical behaviour found in this paper is the following. In the subcritical regime the response is achieved by means of a direct jump from the signal $S$ to the final state $N$. Due to this fact, there is no dependence on $\sigma$ in the probability of response as $N\to\infty$. Instead in the supercritical regime the response is achieved through a transition of the system along the phosphorylation chain. In the second case the process takes place faster due to the fact that $\lambda_{\Delta}(\sigma,E)$ is smaller than $E$. Notice that the direct transition from $S$ to $N$ is very unlikely and is expected to take place in time scales of order $e^{EN}$. However, in the supercritical regime the process is faster and takes place at the time scale $e^{\lambda_{\Delta}(\sigma,E)N}$, that is much shorter than $e^{EN}$ for $N$ large.

1.1 Notation

In this paper we use the notation $\mathbb{R}_{*}=(0,\infty)$ and $\mathbb{N}_{0}=\mathbb{N}\setminus\{0\}$. Moreover we use the notation $f\sim g$ as $N\to\infty$ to indicate that the two functions $f$ and $g$ are asymptotically equivalent, i.e. $\lim_{N\to\infty}\frac{f(N)}{g(N)}=1$. We use the notation $f\approx g$ to indicate that $f$ and $g$ are roughly of the same order. We say that $f\gg g$ as $N\to\infty$ if it holds that $\lim_{N\to\infty}\frac{g(N)}{f(N)}=0$. Finally we use $f\gtrsim g$ to say that there exists a positive constant $c>0$ such that $f\geq cg$.

2 Main result of the paper

In this section we present the model that we study in this paper and the main results that we prove.

2.1 The model

We now explain the assumptions on the chemical rates that we make in this paper. The particular choice of parameters that we consider allows us to define a chemical reaction network in which $\Delta\geq 0$ measures the lack of detailed balance. To define the chemical reaction rates we start associating to each substance in the network an energy. Hence we define a vector of energies $\overline{E}=(E_{S},E_{0},E_{1},\dots,E_{N-1},E_{N})\in\mathbb{R}^{N+2}$. We assume that the attachment of a ligand to the receptor, i.e. the reactions of the form

[TABLE]

take place at rate $K_{R_{1}^{(k)}}=e^{-E_{k}}e^{\sigma}$ for every $k\in\{0,\dots,N\}$. The detachment rate depends on the binding energy between the ligand and the receptor, i.e. the reactions of the form

[TABLE]

take place at rate $K_{-R_{1}^{(k)}}=e^{\sigma}$ for every $k\in\{0,\dots,N\}$ where $\sigma\in\mathbb{R}$. The reactions

[TABLE]

take place at rate $K_{R_{2}^{(k)}}=\alpha e^{\Delta}$ for every $k\in\{0,\dots,N\}$. On the other hand we assume that the reaction

[TABLE]

takes place at rate $K_{-R^{(k)}_{2}}=\alpha e^{E_{k}-E_{k-1}}$ for every $k\in\{0,\dots,N\}$.

In order to simplify the model we make the following assumptions

$E_{S}=0$ and $E_{0}=\sigma$ and $E_{k}=\sigma+kE$ for every $k\in\{1,\dots,N\}$ for $E>0$. As a consequence the rate of the reactions $(k)\rightarrow(k-1)$ is $\alpha e^{E}$. 2. 2.

The ligand is degraded at rate $\mu$ independently on its state. Moreover we assume that the degradation rate $\mu$ is such that

[TABLE]

where $b>0$.

These assumptions guarantee that for every cycle $\omega_{k}=[R_{1}^{(k)},R_{2}^{(k)},-R_{1}^{(k+1)}]$ of the form

[TABLE]

it holds that

[TABLE]

Hence when $\Delta\neq 0$, the detailed balance property does not hold. We stress that we are assuming here that the only reactions that are affected by the parameter $\Delta$ are the phosphorylation reactions $(k)\rightarrow(k+1)$. The fact that detailed balance fails along these reactions implies that the rate $\alpha e^{\Delta}$ of the phosphorylation reactions is larger than the rate that we have when detailed balance holds, i.e. $\alpha$. In other words, the lack of detailed balance accelerates the phosphorylation reactions.

Let $n_{k}(t)$ be the probability that a ligand, characterized by the binding energy $\sigma$, reaches the state $k\in\{0,\dots,N\}$ in the time interval $(0,t]$ and let $n_{S}(t)$ be the probability that it reaches the state $S$ in the time interval $(0,t]$. The dynamics of the vector of the probabilities $n(t)=(n_{S}(t),n_{0}(t),n_{1}(t),\dots,n_{N}(t))^{T}\in\mathbb{R}_{*}^{N+1}$ is described by the following system of ODEs

[TABLE]

with initial datum $n(0)=(1,0,\dots,0)$ and where

[TABLE]

The fact that we are considering as initial datum $n(0)=(1,0,\dots,0)$ means that at time $t=0$ the only substances present in the network are the free ligands. The probability that a ligand, characterised by the binding energy $\sigma$, produces a response in the time interval $(0,t]$, denoted with $R(t)$, is given by

[TABLE]

The probability $p_{deg}(t)$ that a ligand is degraded in the time interval $(0,t]$ is

[TABLE]

where $M(t)$ can be interpreted as the probability that the ligand is the network (either in a complex or as a free ligand) in the time interval $(0,t]$

[TABLE]

Later we refer to $M$ as the total mass of ligands in the network. According to the system of equations introduced above we have that $M$ is decreasing in time. More precisely we have that

[TABLE]

Finally notice that, as expected, we have that

[TABLE]

In this paper we analyse the probability $p_{res}(\sigma)$ that one ligand characterized by its binding energy $\sigma$ with the receptor produce a responses, i.e.

[TABLE]

as the number of kinetic proofreading steps $N\to\infty$.

2.2 Strong discrimination properties for $\Delta>\Delta_{c}$

The goal of the paper is to analyse the relation between the lack of detailed balance and strong discrimination properties of kinetic proofreading mechanisms. We state here precisely our definition of strong discrimination.

Definition 2.1** (Strong discrimination).**

Let $\Delta,\sigma,b,E,\alpha$ be positive constants. The kinetic proofreading network described by the system of ODEs (2.1) has strong discrimination properties with respect to the energy $\sigma$ if

[TABLE]

where $\lambda_{\Delta}(\sigma,E)$ is a parameter that depends non trivially on $\sigma$ and there exists a solution $\sigma_{c}$ to

[TABLE]

In this paper we prove that, for the kinetic proofreading networks introduced in Section 2.1, we have that

[TABLE]

where $\psi$ is a function which depends non trivially on $\sigma$, $E$ and $\Delta$

[TABLE]

and where

[TABLE]

In the latter, for notational convenience, we will remove the subscript $\Delta$ from $\lambda_{\Delta}$ and from $\psi_{\Delta}$.

The fact that $\lambda(\sigma,E)$ does not depend on $\sigma$ for $\Delta$ and $\sigma$ such that $\Delta<\Delta_{c}(\sigma)$ and depends on $\sigma$ for $\Delta$ and $\sigma$ such that $\Delta>\Delta_{c}(\sigma)$ is the most important result of this paper. It implies that the kinetic proofreading networks described by (2.1) can have strong discrimination properties, in the sense of Definition 2.1, only when $\Delta$ and $\sigma$ are such that $\Delta>\Delta_{c}(\sigma)$.

2.3 Necessary conditions on $b$ for strong discrimination when $\Delta>\Delta_{c}$

The goal of this section is to explain how to interpret the consequences of the definition of strong discrimination given in Section 4.3 and in particular to introduce a condition on the parameter $b$ that is necessary in order to have that equation (2.8) admits a solution, hence that the strong discrimination property holds.

The fact that the parameter $\lambda(\sigma,E)$ depends on $\sigma$ is crucial in order to have a sharp transition of $p_{res}(\sigma)$ from $1$ to [math] around a critical binding energy $\sigma_{c}$. Indeed (2.7) implies that the probability $p_{res}(\sigma)$ that one ligand characterized by the binding energy $\sigma$ satisfies

[TABLE]

where $C$ is a positive constant that depends on the parameters. When $\Delta<\Delta_{c}(\sigma)$ we have that $\lambda(\sigma,E)=E$. As a consequence, if $\lambda(\sigma,E)=E<b$, then $\lim_{N\to\infty}p_{res}(\sigma)=1$, while when $\lambda(\sigma,E)=E>b$, then $\lim_{N\to\infty}p_{res}(\sigma)=0$. As will be explained in detail in Section 8.4, when $\Delta\to\infty$ and $\alpha,E,\sigma$ are of order one the function $\psi$ does not depend on $\sigma$ and we have a similar behaviour. The model in this case does not have strong discrimination properties in the sense of Definition 2.1.

Consider a binding energy $\sigma_{c}\in\mathbb{R}$. Assume that $\Delta>\Delta_{c}=\Delta_{c}(\sigma_{c})$. By the continuity of the function $\Delta_{c}(\sigma)$ defined by (2.11) there exists a neighborhood $U$ of $\sigma_{c}$ such that $\lambda(\sigma,E)$ is a strictly increasing function of $\sigma$ in $U$. Then we can find a parameter $b>0$ such that (2.8) holds. Precise conditions that the parameter $b$ must satisfy are given in (4.28). In particular we must have that $b<E$. The reason is that $\lambda(\sigma_{c},E)=\psi(\sigma_{c},E)$ is a decreasing function of $\Delta$. Therefore, in order to have that $(\sigma_{c},\Delta)$ are in the region of strong discrimination in Figure 1 we must have $\psi(\sigma_{c},E)=b<E$.

The condition $b<E$, needed to have strong discrimination has an interpretation: the characteristic life time of the ligands $e^{bN}$ must be strictly smaller than the characteristic time to reach the state $N$ directly from the signal $S$, i.e. $e^{EN}$ as $N\to\infty$.

The fact that there exists a $\sigma_{c}$ such that (2.8) holds implies that there exists a region of transition from response ($p_{res}(\sigma)=1$) to non-response ($p_{res}(\sigma)=0$) when $\sigma$ is such that

[TABLE]

Moreover we have that if $\sigma>\sigma_{c}$ is such that $|\lambda(\sigma,E)-\lambda(\sigma_{c},E)|\gtrsim\frac{1}{N}$, then $\lim_{N\to\infty}p_{res}(\sigma)=0$ while when $\sigma<\sigma_{c}$ is such that $|\lambda(\sigma_{c},E)-\lambda(\sigma,E)|\gtrsim\frac{1}{N}$, then $\lim_{N\to\infty}p_{res}(\sigma)=1$.

In Definition 2.1 we say that the kinetic proofreading model has strong discrimination properties if $\lambda(\sigma,E)$ depends on $\sigma$ and when there exists a solution to (2.8), because in this case there exists a critical binding energy $\sigma_{c}$ which separates the binding energies that produce response from the binding energies that do not produce a response. Moreover the size of the region of transition from response to non-response, in the space of binding energies $\sigma$, where errors in the detection of ligands can occur, is of order $1/N$, hence tends to zero as $N\to\infty$. In Figure 1 we compare the probability of response $p_{res}(\sigma)$ that we obtain selecting the parameters in such a way that $\Delta>\Delta_{c}$ with the probability of response when $\Delta<\Delta_{c}$. Notice that in the latter case the probability is almost constantly equal to [math] while in the first case we have a sharp transition from response to non response around a critical binding energy $\sigma_{c}$.

Finally let us consider the limiting case in which $\Delta\to\infty$ and $\alpha,\sigma$ and $E$ are of order one. Under these assumptions we have that

[TABLE]

Therefore in this regime $\psi(\sigma,E)$ does not depend on $\sigma$, hence we do not have strong discrimination properties. More details on the limiting cases with $\Delta\to\infty$ will be presented in Section 8.

3 Kinetic proofreading models out equilibrium

In this paper we prove that the class of linear kinetic proofreading models described by (2.1) must have a sufficiently large amount of lack of detailed balance in order to have strong discrimination properties. However, at the fundamental level, every chemical reaction network at constant temperature and that does not exchange chemicals with the environment must satisfy the property of detailed balance (see [9, 23]). The kinetic proofreading models that we consider in this paper can be seen as effective systems that approximate the dynamics of larger chemical networks that exchange chemicals with the environment and where the detailed balance property holds. Hence the lack of detailed balance is justified by the fact that the system is assumed to be in contact with reservoirs of chemicals, in particular of molecules of $ATP$, $ADP$ and phosphate groups $P$.

3.1 Kinetic proofreading model involving molecules of ATP, ADP and phosphate groups

In this section we derive the kinetic proofreading model introduced in Section 2.1 starting from a chemical reaction network in which we include explicitly the molecules of $ATP$, $ADP$ and the phosphate groups. This enlarged chemical reaction network satisfies the detailed balance property and reduces to the model studied in this paper if we assume that the concentrations of $ATP$, $ADP$ and $P$ are kept at constant levels by an active process, for instance $ATP$ production by the mitochondria (that we do not model in detail in this paper). We will also explain that the parameter $\Delta$, measuring the lack of detailed balance, depends on the values of the frozen concentrations of $ATP$, $ADP$ and $P$.

Since we are interested in the property of detailed balance it is convenient to write the linear kinetic proofreading model described in Subsection 2.1 as a chain of cycles $\{\omega_{k}\}_{k=0}^{N}$. These cycles have the following form

[TABLE]

Notice that if $\Delta=0$, then we have that the Wegscheider condition (1.1) holds for every cycle, hence the property of detailed balance holds for the chemical reaction network.

In order to obtain the cycles $\{\omega_{k}\}_{k=0}^{N}$ in which the Wegscheider criterium fails, starting from cycles for which the Wegscheider criterium holds, we consider the chemical reaction network with substances $\Omega_{c}=\Omega\cup\{ATP,ADP,P,S\}$. We assume that these substances interact via the chemical reactions $\{R_{1}^{(k)},-R_{1}^{(k)}\}_{k=0}^{N}$, $R_{2},-R_{2}$, $\{R_{3}^{(k)},-R_{3}^{(k)}\}_{k=0}^{N}$ defined as

[TABLE]

for $k\in\{0,\dots,N\}.$ Here we are using the notation

[TABLE]

We assume that the chemical rates of the reactions $\{R_{1}^{(k)},-R_{1}^{(k)}\}_{k=0}^{N}$ are independent on $k=0,\dots N$. More precisely, the chemical rates of the reactions $R_{1}^{(k)},R_{2}$ and of the reverse reactions $-R_{1}^{(k)}$ and $-R_{2}$ are given by

[TABLE]

for some constants $E_{T},E_{D},E_{P}\in\mathbb{R}$. The chemical rates of the reactions $R_{3}^{(k)}$ and of the reverse reactions $R_{-3}^{(k)}$ are given by

[TABLE]

Cycles of the enlarged chemical reaction network.

Notice that if we apply the reactions $R^{(k)}_{1},R_{-2}$, $R_{-3}^{(k)},R_{3}^{(k+1)}$ to any vector of concentrations, then we will have a null effect, indeed

[TABLE]

Therefore the set of the reactions $\overline{\omega}_{k}=\{R_{1}^{(k)},R_{-2},R_{-3}^{(k)},R_{3}^{(k+1)}\}$ for $k=0,1,\dots,N$ are cycles of the enlarged chemical reaction network. Notice that by construction we have that the circuit condition (1.1) holds along every cycle $\overline{\omega}_{k}$, indeed the rates were selected in such a way that

[TABLE]

As a consequence, this enlarged chemical reaction network satisfies the property of detailed balance.

System of ODEs associated with the enlarged chemical reaction network and its conserved quantities.

The system of ODEs associated with the chemical network that includes ATP, ADP and P is the following,

[TABLE]

Notice that the vector $\textbf{E}\in\mathbb{R}^{N+4}$ defined as

[TABLE]

is such that $e^{-\textbf{E}}$ a steady state for the system of equations (3.1).

Moreover notice that we have some conserved quantities in our system. In particular we have that if $n(t)=(n_{0},n_{1},\dots,n_{N-1},n_{N},n_{T},n_{D},n_{P},n_{S})\in\mathbb{R}^{N+4}$ is the vector of the concentrations, then we have that

[TABLE]

where $m_{1}=(\overbrace{1,1,\dots,1,1}^{N\text{ times}},0,0,0,1)^{T}$ and $m_{2}=(0,1,2,3,\dots,N-1,N,2,1,1,0)^{T}$. In other words we have that the number of phosphate groups is constant in time, hence

[TABLE]

and the number of ligands is constant in time, i.e.

[TABLE]

Notice that all the conserved quantities, i.e. all the vectors $m\in\mathbb{R}^{N+4}$ such that $m^{T}n(t)=m^{T}n(0)$, can be written as linear combinations of $m_{1}$ and of $m_{2}$.

Since the chemical reaction network (3.1) satisfies the property of detailed balance all the steady states of (3.1) are of the form

[TABLE]

for suitable constants $\mu_{1},\mu_{2}$, see Lemma 3.10 in [9].

Therefore the concentrations of $ATP$, $ADP$ and $P$ at any steady state $N\in\mathbb{R}_{+}^{N+4}$ of (3.1), denoted respectively by $N_{T},N_{D},N_{p}$, satisfy the following equality

[TABLE]

3.2 Lack of detailed balance for constant non equilibrium concentration of ATP

We now assume that the concentrations of $ATP$, $ADP$ and $P$ are constant in time. This can be justified when the chemical network (3.1) is in contact with a reservoirs of $ATP$, $ADP$ and $P$. Let $\overline{n}_{T}$, $\overline{n}_{D}$ and $\overline{n}_{P}$ be the constant concentrations of $ATP$, $ADP$ and $P$ respectively. Let us define $\Delta$ as

[TABLE]

Since the concentrations of $ATP$, $ADP$ and $P$ are constant in time, the chemical reactions in (3.2) can be reduced to

[TABLE]

where we assume that the chemical rates $\overline{K}_{1},\overline{K}_{-1}$, $\overline{K}_{3},\overline{K}_{-3}^{(k)}$ of the reactions depend on $\overline{n}_{T},\overline{n}_{D},\overline{n}_{P}$, more precisely

[TABLE]

Notice that if we choose that the concentrations of $ATP$, $ADP$ and $P$ are given by

[TABLE]

then the chemical reaction network (3.4) satisfies the property of detailed balance if and only if $\Delta=0$. Moreover notice that, under these assumptions on the frozen concentrations, we have that the chemical network (3.4) and the chemical reaction network (3.1) coincide.

Notice that the concentration of molecules of ATP, ADP and phosphate groups are assumed to be constant in time because the system is in contact with reservoirs of these substances. Therefore we must have fluxes of ATP, ADP and P between the chemical network and the environment. We can rewrite the equations for the change in time of ATP, ADP and P in equation (3.1) taking into account of these fluxes. The equations that we obtain are the following

[TABLE]

In the equations above we have used the fact that $\overline{n}_{P}=e^{-E_{P}}$, $\overline{n}_{D}=e^{-E_{D}}$ and $\overline{n}_{T}=e^{-E_{T}+\Delta}$. Evaluating the fluxes at the steady states $N$ with $N_{T}=e^{-E_{T}}$, $N_{D}=e^{-E_{D}}$, $N_{P}=e^{-E_{P}}$ that the fluxes must be given by

[TABLE]

The fluxes $J^{ext}_{T}(N)$, $J^{ext}_{D}(N)$, $J^{ext}_{P}(N)$ can be used in order to compute the energy needed in order to main the system at steady state.

4 Matched asymptotic expansions

In this section we explain, using formal arguments, why the kinetic proofreading model has strong discrimination properties if and only if $\Delta>\Delta_{c}$.

In order to do this we study the probability of response $p_{res}(\sigma)$ defined as in (2.6) for large values of $N$. Therefore we analyse the behaviour of $n_{N}$ as $N\to\infty$. To this end we use the method of matched asymptotic expansion. More precisely we study the behaviour of $n_{k}$ as $N\to\infty$ in two main regions: the outer region, where $k$ is such that $N-k=m$ and $m\gg 1$ as $N\to\infty$ and the inner region, where $k$ is such that $N-k\approx 1$ as $N\to\infty$. To obtain the asymptotic behaviour of $n_{N}$ as $N\to\infty$ we match the behaviours that we obtain in the inner and in the outer region.

Notice that in this paper we study the discrimination properties of the linear kinetic proofreading model when the number of kinetic proofreading steps $N$ tends to infinity. The number of kinetic proofreading steps has been measured experimentally by many authors and has been estimated to be between 2 and 20 (see for instance [1, 4, 21, 26]). Therefore the discrimination property in concrete biological situations will not be as strong as in the model analysed here. However, considering the limit as $N\to\infty$ allows to understand in detail how the kinetic proofreading mechanism performs its functions. In particular allows us to prove that a critical lack of detailed balance is needed in order to have strong discrimination properties.

Before analysing the behaviour of $n_{k}$ in the inner and in the outer region we study how to approximate the concentration of ligands, $n_{S}$ as $N\to\infty$. To this end we consider the limit as $N\to\infty$ in the equation for $n_{S}$ in (2.1) and using the fact that $\mu$ satisfies (2.1) we obtain

[TABLE]

where we recall that $M(t)$ is defined as (2.4).

As we will show later, the concentrations $n_{k}$ stabilize to exponential solutions that behave as $e^{-\lambda(\sigma,E)k}$ for $k\approx N$ and $N$ large, at time scales that are of order $N$. Here $\lambda(\sigma,E)$ is positive and depends also on $\Delta$. We now argue heuristically to infer that the probabilities $n_{k}$ can be approximated by quasi stationary solutions $\tilde{n}_{k}$ solving

[TABLE]

and such that their total mass decreases in time scales $t$ larger than $N$. Indeed, assume that $k=Nx$ and that $g(t,x)=\frac{n_{k}(t)}{\tilde{n}_{k}}$. Then, in order to obtain the time scale for the stabilization of $n_{k}$ we approximate the discrete equation in (2.1) with the following PDE

[TABLE]

Here we are using the rough approximation $\frac{\tilde{n}_{k}}{n_{k}}\approx 1$ in order to approximate the incremental quotient by derivatives. Notice that the time scales at which the characteristic curves of this equation stabilize for $t\approx N$. Moreover it stabilizes to $\frac{n_{S}(t)}{\tilde{n}_{Nx}}e^{-NxE}$. In Section 7 we explain in detail under which conditions the discrete system of ODEs (2.1) can be approximated by means of PDEs. A more detailed analysis of the time dependent discrete equation (2.1) gives the same order of magnitude for the stabilization time. In Section 6 we obtain rigorously that the stabilization of $n_{k}e^{-k\lambda(\sigma,E)}$ takes place for times of order $N$.

In particular the equation for $M$, i.e. (2.5), implies that the variation of mass occurs at time scales of order $e^{cN}$ for a suitable positive constant $c>0$ that depends on $E$, $\Delta$ and on $\sigma$. Notice also that $n_{S}$ stabilizes in times $t$ of order $1$. Therefore in order to obtain the variation of $(n_{S},n_{k})_{k=0}^{N}$ in times scales of order $N$ we can use a quasi steady state approximation to replace $(n_{S},n_{k})_{k=0}^{N}$ by a steady state having mass $M(t)$ and to use the equation for $n_{N}$ in (2.1) in order to compute the change of the mass in exponentially long time scales. Later we will show that these assumptions are consistent with the results that we obtain via this heuristic argument, see (4.25).

In particular, the quasi steady state approximation for $n_{S}$ implies that $n_{S}\approx\tilde{n}_{S}$ where

[TABLE]

where

[TABLE]

4.1 Outer region: $N-k\gg 1$ as $N\to\infty$.

In this section, we study the behaviour of $n_{k}$ when $N\to\infty$ for $k\approx 1$ such that $k=N-m$ with $m\to\infty$ as $N\to\infty$. Since in this region we have that $k\ll N$, we assume that the dynamics is given by the system of equations obtained taking the limit as $N\to\infty$ in (2.1), i.e.

[TABLE]

The initial datum for equation (4.2) is the same initial datum that we have for the original system of ODEs (2.1), i.e. we have $n^{(out)}(0)=(1,0,\dots,0)$. Finally notice also that

[TABLE]

As before, we assume that $n_{S}$ and $\{n^{(out)}_{k}\}_{k=0}^{\infty}$ can be approximated by quasi stationary solutions $\tilde{n}_{S}$ and $\{\tilde{n}^{(out)}_{k}\}_{k=0}^{\infty}$. Later we will see that this assumption is consistent with our results. Therefore we look for solutions to the following system of equations

[TABLE]

In order to study the solutions to this equation we consider separately the case of $\Delta=\Delta_{c}$ and $\Delta\neq\Delta_{c}$.

As a first step we compute $\tilde{n}_{k}^{(out)}$ for $\Delta\neq\Delta_{c}$. Then we will see that, as $N\to\infty$, we have different behaviours for $\tilde{n}_{k}^{(out)}$ depending on whether $\Delta>\Delta_{c}$ or $\Delta<\Delta_{c}$. For completeness we also study the case $\Delta=\Delta_{c}$.

Explicit solution to (4.3) in the subcritical and supercritical case $\Delta\neq\Delta_{c}$.

Notice that the solution to equation (4.3) can be written as a linear combination of a particular solution and of the homogeneous solutions. Here with homogeneous solutions we mean solutions to equation (4.3)-(4.4) where $\tilde{n}_{S}=0$, i.e.

[TABLE]

It is easy to see that a particular solution to (4.4) is

[TABLE]

where the constant $A(0)$ is given by

[TABLE]

Notice that since we have that $\Delta\neq\Delta_{c}$, then we have that $e^{\sigma}-\alpha(e^{\Delta}-1)(e^{E}-1)\neq 0.$

We now look for solutions to (4.5) of the form $n_{k}^{hom}=\theta(0)^{k}$. Substituting $n_{k}^{hom}$ in equation (4.5) we deduce that $\theta(0)$ must satisfy the following equation

[TABLE]

where

[TABLE]

As a consequence we have that the solutions $\theta_{1}(0)$ and $\theta_{2}(0)$ to the equation (4.6) are

[TABLE]

Now notice that

[TABLE]

Indeed we have that the function $f(\theta)=\theta^{2}-\frac{\Omega(0)}{\alpha e^{E}}\theta+e^{\Delta-E}$ is such that $f(1)=-e^{\sigma}<0$. Since we have that $M^{(out)}(t)$ is finite, then the only admissible solution is $\theta_{1}(0)$, indeed $\theta_{2}^{k}(0)\rightarrow\infty$ as $k\to\infty$.

We deduce that

[TABLE]

where $F$ is a suitable function that varies slowly in time, as a consequence of the fact that, as explained at the beginning of Section 4, the time scale at which the total mass of chemicals in the system starts to decrease is exponential in $N$. In order to obtain the function $F$ we notice that the formula above implies that

[TABLE]

On the other hand, equation (4.3) together with the fact $\tilde{n}_{1}^{(out)}(t)=e^{-E}A(0)\tilde{n}_{S}(t)+F(t)\theta_{1}(0)$ implies that

[TABLE]

We deduce that

[TABLE]

where we have that

[TABLE]

Summarizing, we obtain the following approximation for $n_{k}^{(out)}$ valid when $t\approx N$

[TABLE]

for

[TABLE]

Notice that the quasi steady state approximation $n_{k}^{(out)}(t)\approx\tilde{n}_{k}^{(out)}(t)$ and $n_{S}(t)\approx\tilde{n}_{S}(t)$ is valid for $t\approx N$.

4.1.1 Asymptotic behaviour of $n_{k}^{(out)}$ when $\Delta<\Delta_{c}$ and $k\gg 1$, $N-k\gg 1$.

We now want to study the asymptotic behaviour of $n_{k}^{(out)}$ as $k\to\infty$. To this end we first of all assume that $\Delta<\Delta_{c}$. A direct computation shows that in this case we have that $\varphi(\sigma)<1$. As a consequence we deduce that, in the subcritical case (i.e. when $\Delta<\Delta_{c}$) we have that

[TABLE]

Using the quasi steady state approximation we deduce that $n_{k}^{(out)}(t)\sim A(0)\overline{n}_{S}M(t)e^{-kE}$ as $k\gg 1$, $N-k\gg 1$ and $t\approx N$.

4.1.2 Asymptotic behaviour of $n_{k}^{(out)}$ when $\Delta>\Delta_{c}$ and $k\gg 1$, $N-k\gg 1$.

On the contrary, in the supercritical case in which $\Delta>\Delta_{c}$, we have that $\varphi(\sigma)>1$, hence

[TABLE]

where we recall that $B$ is defined by (4.9) and $\psi$ is given by (2.10) and it depends non trivially on $\sigma$. The quasi steady state approximation implies that $n_{k}^{(out)}(t)\sim A(0)\overline{n}_{S}M(t)Be^{-k\psi(\sigma,E)}$ as $k\gg 1$, $N-k\gg 1$ and $t\approx N$.

Notice that both in the subcritical and in the supercritical case we have that the asymptotic behaviour of $n_{k}^{(out)}$ as $k\to\infty$ is given by $F(t)e^{-\lambda(\sigma,E)k}$ for a suitable function $F(t)$ encoding the change of mass in the system and a suitable exponent $\lambda(\sigma,E)$. The main difference between the supercritical and the subcritical case is that in the first case, when $\Delta>\Delta_{c}$, we have that $\lambda(\sigma,E)=E-\log(\varphi(\sigma))$, hence $\lambda(\sigma,E)$ depends on the binding energy $\sigma$. Instead in the case $\Delta<\Delta_{c}$ we have that $\lambda(\sigma,E)$ does not depend on the binding energy $\sigma$. Indeed in that case we have that $\lambda(\sigma,E)=E$. As we will see later in Subsection 4.3 the fact that $\lambda(\sigma,E)$ depends on $\sigma$ when $\Delta>\Delta_{c}$ is crucial in order to have strong discrimination properties.

4.1.3 Asymptotic behaviour in the critical case $\Delta=\Delta_{c}$ when $k\gg 1$ and $N-k\gg 1$

By completeness we also study the behaviour of $n_{k}$ as $k\to\infty$ in the critical case $\Delta=\Delta_{c}$. Under this assumption, we have that $\alpha(e^{\Delta}-1)(e^{E}-1)=e^{\sigma}$. As in the case $\Delta\neq\Delta_{c}$ we have that the behaviour in the outer region is described by equation (4.3). We will now compute an explicit solution.

Explicit solution to (4.3).

First of all we notice that when $\Delta=\Delta_{c}$ we have that

[TABLE]

is a particular solution to (4.4). We now aim at finding a solution to the homogeneous equation corresponding to (4.4), i.e. to the following equation

[TABLE]

In this case we have homogeneous solutions of the form $n_{k}^{(hom)}=\theta(0)^{k}$ where $\theta(0)$ is a solution to

[TABLE]

Notice that the two solutions $\theta_{1}(0)$ and $\theta_{2}(0)$ to the second order equation are given by

[TABLE]

As before, the only admissible solution is $\theta_{1}(0)<1$. Therefore, the solution to equation (4.4) is

[TABLE]

where $F(t)=\tilde{n}_{0}^{(out)}(t)$ and where $\tilde{n}_{0}^{(out)}(t)$ can be computed by solving equation (4.3). Indeed we have that

[TABLE]

which implies that

[TABLE]

Asymptotic behaviour of $n_{k}^{(out)}$ as $k\to\infty$ when $\Delta=\Delta_{c}$.

Equality (4.12) implies that

[TABLE]

The quasi steady state approximation implies that $n_{k}^{(out)}(t)\approx\frac{ke^{-kE}\tilde{n}_{S}(t)}{\alpha(e^{\Delta+E}-1)}$ as $k\gg 1$ and $N-k\gg 1$ and for $t\approx N$. As in the subcritical case $\Delta<\Delta_{c}$, in this case we have that $\lambda(\sigma,E)=E$, hence $\lambda(\sigma,E)$ does not depend on the binding energy $\sigma$.

4.2 Inner region: $N-k\approx 1$ as $N\to\infty$

We now study the behaviour of $n_{k}$ in the region in which $k=N-m$ and $m\approx 1$. In this region $k$ is large and the dynamics is described by the equation for $n_{k}$ in (2.1), i.e.

[TABLE]

Notice that we are rewriting the equation for $n_{N}$ in (2.1) introducing artificially the state $N+1$ in the model and imposing that $n_{N+1}=0$. This allows to have the same equation to describe the evolution of $n_{k}$ for $k<N$ and for $k=N$.

Recall that $M$ changes in exponentially long times, i.e. for $t\approx e^{cN}$ for a positive constant $c$. This allows us to approximate $n_{k}$ with a quasi steady state solution $\tilde{n}_{k}$ Since we assume that $\mu\sim e^{-bN}$ as $N\to\infty$, this quasi steady state approximation for the solution to equation (4.14) satisfies the following system of equations

[TABLE]

Notice that the condition $\tilde{n}_{N+1}=0$ plays the role of an absorbing boundary conditions, indeed when a complex reaches the state $N+1$, it is lost. This boundary condition, together with the degradation term $\mu$, is the reason why the total mass in our system slowly decreases in time when $t\approx e^{cN}$. In order to analyse the solution to (4.15) we consider separately the subcritical case from the supercritical case.

4.2.1 Subcritical case $\Delta<\Delta_{c}$.

In order to study the solution to (4.15) it is convenient to consider the following change of variables

[TABLE]

where we recall that $N-k=m$. Therefore we now study the following equation for $h_{m}$

[TABLE]

We stress that the condition $h_{-1}(t)=0$ corresponds to $\tilde{n}_{N+1}(t)=0$. Instead the condition on the limit of $h_{m}$ as $m\to\infty$ is a consequence of the matching with the behaviour of $\tilde{n}_{k}$ in the outer region, i.e. (4.10).

Notice that $h^{(p)}_{m}=A(0)n_{S}(t)$ is a particular solution to equation (4.17). We now study the homogeneous solutions to equation (4.17), i.e. solutions to (4.17) without the source term. In particular we look for homogeneous solutions of the form

[TABLE]

Substituting this ansatz in equation (4.17) we obtain that $\Theta$ must satisfy the following equation

[TABLE]

This equation has two solutions $\Theta_{1}$ and $\Theta_{2}$ given by

[TABLE]

Notice that since we expect the limit as $m\to\infty$ of $h_{m}$ to be finite for every fixed time $t>0$ then $\Theta^{m}$ is an admissible homogeneous solution only when $\Theta\leq 1$. Moreover, since $\Delta<\Delta_{c}$, then we have that

[TABLE]

Hence there exists a unique admissible homogeneous solution, which is

[TABLE]

As a consequence, the solution to (4.17) is given by

[TABLE]

Moreover, the condition $h_{-1}=0$, implies that

[TABLE]

Summarizing we have that

[TABLE]

Hence we obtain that

[TABLE]

Using the quasi steady state approximation $\tilde{n}_{N}(t)\approx n_{N}(t)$ for $t\approx N$ we deduce that

[TABLE]

Here we stress again that the asymptotic behaviour of $n_{N}$ as $N\to\infty$ is of the form $F(t)e^{-N\lambda(\sigma,E)}$ where the parameter $\lambda(\sigma,E)$, in this subcritical case $\Delta<\Delta_{c}$, does not depend on $\sigma$ and where $F(t)$ is a suitable function of time.

4.2.2 Supercritical case $\Delta>\Delta_{c}$.

As we did in the subcritical case (see (4.16)) we remove the leading exponential behaviour by using the change of variables

[TABLE]

where $\lambda$ is defined as (2.9). Performing the change of variable in equation (4.15) we deduce that $\ell_{k}$ must satisfy the following equation

[TABLE]

Since we are analysing the behaviour in the inner region where $k\gg 1$ we have that $\tilde{n}_{S}\ll e^{k\log(\varphi(\sigma))}\ell_{k}$, $\tilde{n}_{S}\ll e^{\Delta+E+(k-1)\log(\varphi(\sigma))}\ell_{k-1}$, and $\tilde{n}_{S}\ll e^{(k+1)\log(\varphi(\sigma))}\ell_{k+1}$. Therefore we can approximate the equation above as follows

[TABLE]

Moreover, we define $h_{m}=\ell_{N-m}$. Then $h_{m}$ satisfies the following equation

[TABLE]

where the last equality comes from the matching condition, indeed we have that for $k\gg 1$ and $N-k\gg 1$ we have that the behaviour of $\tilde{n}_{k}$ is given by (4.11).

We make the ansatz $h_{m}(t)=\Theta^{m}$ for every $m\in\mathbb{N}$ and obtain that

[TABLE]

We deduce that

[TABLE]

In particular we obtain that

[TABLE]

and

[TABLE]

Both these solutions are admissible as they are both smaller or equal than one. We deduce that

[TABLE]

Matching the behaviour in the inner with the behaviour in the outer region implies that $F_{1}(t)=BA(0)n_{S}(t)$. On the other hand, since $h_{-1}(t)=0$, we also have that

[TABLE]

where we recall that $B$ is given by (4.9). Summarizing we have that

[TABLE]

therefore we have that when $t\approx N$ it holds that

[TABLE]

Notice that the asymptotic behaviour of $n_{N}$ as $N\to\infty$ is of the form $F(t)e^{-N\lambda(\sigma,E)}$ where $\lambda(\sigma,E)$ depends on $\sigma$ where $F(t)$ is a function of time that takes into account the change of mass, notice that this function changes very slowly in time, more precisely in time scales much larger than $N$.

4.2.3 Critical case $\Delta=\Delta_{c}$.

For completeness, we now analyse the behaviour of $\tilde{n}_{k}$ in the inner region under the assumption that $\Delta=\Delta_{c}$. To this end we make the change of variables $\ell_{k}:=\frac{e^{kE}}{k}\tilde{n}_{k}$ in equation (4.15) and deduce that $\ell_{k}$ satisfies

[TABLE]

Since in the inner region we have that $k\gg 1$, then $\tilde{n}_{S}$ is of lower order with respect to the other terms in the equation. Therefore, we approximate equation (4.20) with the following homogeneous equation

[TABLE]

In order to match the behaviour in the inner region with the behaviour in the outer region it is convenient to consider the change of variable $k=N-m$ and to study $h_{m}:=(N-m)\ell_{N-m}$. Notice that, due to its definition, $h_{m}$ satisfies the following equation

[TABLE]

Notice that the last condition comes from the matching with the outer region, i.e. by the fact that the behaviour in the outer region is given by (4.13).

Since the equation for $h_{m}$ is homogeneous it has solutions of the form $h_{m}=\theta^{m}$ where $\theta$ satisfies

[TABLE]

We deduce that both $\theta_{1}^{m}=1$ and $\theta_{2}^{m}=e^{-m(\Delta+E)}$ are admissible solutions to (4.22) . Therefore we have that the sequence $\{h_{m}\}_{m}$ given by

[TABLE]

with $F_{2}=-F_{1}e^{-(\Delta+E)}$ and with

[TABLE]

is the solution to (4.22). In particular we conclude that $h_{0}=\frac{n_{S}(t)}{\alpha e^{\Delta+E}}$ and therefore

[TABLE]

4.3 Discrimination property

We now study the discrimination property of the model under different assumptions on $\Delta$. To this end we analyse the probability of response $p_{res}$ defined as in (2.6) as $N\to\infty$ in the different regimes, subcritical, critical and supercritical using the asymptotics (4.18), (4.19), (4.23) for $n_{N}$.

Subcritical case $\Delta<\Delta_{c}$.

Using the asymptotics for $n_{N}$ in (4.18) we approximate the probability of response as follows

[TABLE]

In the approximation above, we are using that the region where $t\approx N$ gives the main contribution to the integral $\int_{0}^{\infty}n_{N}(t)dt$. Moreover, notice that we have that

[TABLE]

Notice that equation (4.24) implies that the mass changes in exponentially large time scales. More precisely we have that

[TABLE]

Using this approximation for $M$ we deduce that

[TABLE]

Using the fact that $\mu\sim e^{-bN}$ as $N\to\infty$ we obtain that

[TABLE]

where $C_{1}:=\left[\alpha e^{\Delta}A(0)\overline{n}_{S}\left(1-\varphi(\sigma)e^{-E-\Delta}\right)\right]^{-1}$. In particular we conclude that $p_{N}\approx 0$ when $E>b$, $p_{N}\approx 1$ as $N\to\infty$ if $E<b$. Notice that in this case if we plot the probability $p_{res}(\sigma)$, defined as (2.6), as a function of $\sigma$ we do not have a sharp transition from [math] to $1$ around a critical binding energy $\sigma_{c}$ (see Figure 2). The dependence on $\sigma$ is encoded only in the constant $C_{1}$.

Critical case $\Delta=\Delta_{c}$.

Arguing as in the subcritical case $\Delta<\Delta_{c}$ we deduce that (4.23) implies

[TABLE]

We conclude that $p_{res}(\sigma)\approx 0$ as $N\to\infty$ when $E>b$, while $p_{res}(\sigma)\approx 1$ as $N\to\infty$ if $E\leq b$. The situation is similar to the one that we have for the subcritical case $\Delta<\Delta_{c}$.

Supercritical case $\Delta>\Delta_{c}$

Using the approximation for $n_{N}$ as $N\to\infty$ given in (4.19) we deduce that

[TABLE]

where $C_{2}=\left[\alpha e^{\Delta}BA(0)\overline{n}_{S}\left(1-(\varphi(\sigma))^{2}e^{-(\Delta+E)}\right)\right]^{-1}$. Recall that $\varphi(\sigma)$ is a function of $\sigma$, indeed

[TABLE]

Moreover notice that

[TABLE]

hence the function $\varphi$ is strictly decreasing. In particular this implies that the function $\lambda(\sigma,E)=E-\log(\varphi(\sigma))$ is increasing in $\sigma$. As a consequence if we assume that

[TABLE]

then we have that there exists $\sigma_{c}\in\mathbb{R}$ such that $\lambda(\sigma_{c},E)=b$. Notice that if we consider $p_{res}$ as a function of $\sigma$ we see a sharp transition from $1$ to [math], i.e. from response to non-response, around the critical binding energy $\sigma_{c}$. The region in which the transition takes place is the region in which we have that

[TABLE]

This region tends to [math] as $N\to\infty$. If instead we have that $\sigma$ is such that $\lambda(\sigma,E)-\lambda(\sigma_{c},E)>\frac{1}{N}$ then we have $p_{res}(\sigma)\approx 1$ as $N\to\infty$ and if $\lambda(\sigma_{c},E)-\lambda(\sigma,E)>\frac{1}{N}$ then we have $p_{res}(\sigma)\approx 0$ as $N\to\infty$.

Notice that the discrimination property that we obtain for $\Delta>\Delta_{c}$ relies on the fact that we assume that ligands are degraded at a rate $\mu$ that satisfy (2.1). If we assume that ligands are not degraded, hence $\mu=0$, the discrimination property is lost.

5 Approximated model on $\mathbb{N}$

In this Section we study in detail the model that we obtain taking the limit as $N\to\infty$ in (2.1). Analysing this model rigorously using Laplace transform methods is easier than analysing the model with a finite number of states $N$. Moreover the behaviour of the model on the half line gives insights on the behaviour of the model with a finite number of states $N$ when $N$ is large. Indeed we prove that if $\Delta<\Delta_{c}$, then the solution to the problem on the half line is such that as $t\to\infty$

[TABLE]

We refer to Theorem 5.8 for the rigorous statement. If instead $\Delta>\Delta_{c}$, then

[TABLE]

while

[TABLE]

We refer to Theorem 5.9 for the rigorous statement. Notice that the results that we obtain for this model are consistent with the formal arguments in Section 4. Moreover notice that, when $k\approx t$ we have that $n_{k}(t)$ is at steady state. This is consistent with the assumption that $n_{N}(t)$ stabilizes in time scales of order $N$ that we make in Section 4.

The system of ODEs that we study in this section is the one that we obtain taking the limit as $N\to\infty$ in (2.1), i.e.

[TABLE]

with initial datum $n(0)=(1,0,\dots,0)$. Notice that since $(n_{S},n_{0},\dots,n_{N})$ is the vector of the probabilities of reaching the different states of the system in the time interval $(0,t]$, we have that $M(t)=n_{S}(t)+\sum_{k=0}^{\infty}n_{k}(t)=1.$

Taking Laplace transforms in all the terms in (5) we obtain the following systems of equations for the vector of Laplace transforms $(\hat{n}_{S}(z),\hat{n}_{0}(z),\hat{n}_{1}(z),\dots,\hat{n}_{N-1}(z))$ of $(n_{S}(t),n_{0}(t),n_{1}(t),\dots,n_{N-1}(t))$

[TABLE]

where $z\in\mathbb{C}$ and where we use the notation

[TABLE]

We will compute now explicit solutions to (5). As a first step we compute $\hat{n}_{S}$.

Lemma 5.1**.**

Assume that the vector $(\hat{n}_{S}(z),\hat{n}_{0}(z),\hat{n}_{1}(z),\dots,\hat{n}_{N-1}(z))$ is a solution to (5). Then the function $\hat{n}_{S}$ is analytic for every $z\in\mathbb{C}\setminus\{z_{S},0\}$ where

[TABLE]

More precisely, we have that

[TABLE]

Proof.

Since in this model we have that $\frac{dM(t)}{dt}=0$ we deduce that

[TABLE]

Hence the equation for $\hat{n}_{S}$ can be rewritten as

[TABLE]

This implies that (5.5) holds. ∎

We find now a particular solution to (5). To this end it is convenient to define the function $A:\mathbb{C}\rightarrow\mathbb{C}$, which is given by

[TABLE]

where

[TABLE]

Lemma 5.2**.**

Assume that the function $\hat{n}_{0}:\mathbb{C}\rightarrow\mathbb{C}$ is given and let $\hat{n}_{S}$ be given by (5.5). We define the sequence $\{\hat{n}_{k}^{(p)}\}_{k=1}^{\infty}$ of functions $\hat{n}_{k}^{(p)}:\mathbb{C}\rightarrow\mathbb{C}$ given by

[TABLE]

The sequence $\{\hat{n}_{k}^{(p)}\}_{k=1}^{\infty}$ is a particular solution to

[TABLE]

Proof.

We substitute (5.8) in the equation for $\hat{n}_{k}$ and obtain that

[TABLE]

Using (5.6) we deduce that

[TABLE]

∎

Notice that for every $k\in\mathbb{N}_{0}$ the function $\hat{n}_{k}^{(p)}$ is analytic for every $z\in\mathbb{C}\setminus\{z_{A},z_{S}\}$.

We now study the homogeneous solutions of equation (5.9), i.e. we find the solutions to the following equation

[TABLE]

Lemma 5.3**.**

Assume that $\hat{n}_{0}:\mathbb{C}\rightarrow\mathbb{C}$ is given and let $\hat{n}_{S}$ be given by (5.5). Then equation (5.10) admits two solutions $\theta_{1}^{k}$ and $\theta_{2}^{k}$ where

[TABLE]

Proof.

Notice that the homogeneous equation (5.10) is a linear finite difference equation of order $2.$ Therefore it has two linearly independent solutions. Moreover the solutions are of the form $\theta^{k}$. Substituting $n_{k}^{(h)}=\theta^{k}$ in (5.10) we deduce that

[TABLE]

Hence the desired result follows. ∎

Lemma 5.4**.**

Let $\theta_{1}:\mathbb{C}\rightarrow\mathbb{C}$ and $\theta_{2}:\mathbb{C}\rightarrow\mathbb{C}$ be given by (5.11). Then we have that

[TABLE]

Proof.

We start by proving that $|\theta_{2}(z)|\geq 1$. To this end it is convenient to prove that if $\rho\geq R_{0}$ such that $R_{0}\leq 1$ and $\psi\in(0,2\pi)$ then it holds that

[TABLE]

To prove this we first of all notice that $\xi:=\sqrt{1+\rho e^{i\psi}}$ is such that $\mathop{}\!\mathrm{Re}\,(\xi)\geq\sqrt{1-R_{0}}$. Indeed, let $\xi=\xi_{1}+i\xi_{2}$ for $\xi_{1},\xi_{2}\in\mathbb{R}$. This implies that $1+\rho\cos(\psi)+i\rho\sin(\psi)=\xi_{1}^{2}-\xi_{2}^{2}+2i\xi_{1}\xi_{2}$. Hence we deduce that $\rho\sin(\psi)=2\xi_{1}\xi_{2}$ and $1+\rho\cos(\psi)=\xi_{1}^{2}-\xi_{2}^{2}$. In particular this implies that

[TABLE]

and the desired conclusion follows. Using the fact that $\mathop{}\!\mathrm{Re}\,(\xi)\geq\sqrt{1-R_{0}}$ we obtain that

[TABLE]

In order to apply (5.13) to prove that $|\theta_{2}(z)|>1$ it is convenient to use the notation $\beta:=e^{E-\Delta}$ and $H(z):=\frac{\Omega(z)}{\alpha(1+\beta)}$. We deduce that

[TABLE]

Let us assume now that $\Delta\geq E$, hence $\beta\leq 1$. In this case,

[TABLE]

Hence we choose $R_{0}=\frac{4}{1+\frac{2}{\beta}+\frac{1}{\beta^{2}}}$ and using (5.13) we deduce that $|\theta_{2}(z)|\geq\frac{1+\beta}{2\beta}|H|\geq 1$.

Assume now that $\beta\geq 1$ because $E\geq\Delta$. Define $\gamma=\frac{1}{\beta}\leq 1$. Notice that

[TABLE]

Therefore we can repeat the argument above used to prove that $|\theta_{2}(z)|\geq 1$ when $\beta\leq 1$.

Let us prove now that $|\theta_{1}(z)|=\left|\frac{1+\beta}{2\beta}H\left(1-\sqrt{1-\frac{4\beta}{(1+\beta)^{2}H^{2}}}\right)\right|\leq 1$. Notice that $|\theta_{1}(z)\theta_{2}(z)|=\frac{2}{1+\beta}$. Assume now that $\beta\geq 1$. In this case we have that $|\theta_{1}\theta_{2}|\leq 1$ and $|\theta_{2}|\geq 1$. As a consequence we must have $|\theta_{1}|\leq 1$. Assume now that $\beta\leq 1$. Let us consider as before $\gamma=\frac{1}{\beta}\geq 1$. We can now write

[TABLE]

and use the argument above to deduce that $|\theta_{1}(z)|\leq 1$. ∎

Remark 5.5**.**

In this paper we use the following definition of square root function in the domain $\mathbb{C}\setminus(-\infty,0)$

[TABLE]

where $z=|z|e^{i\psi}$ with $\psi\in(-\pi,\pi]$.

The solution to (5) can be written as the sum of the particular solution and the homogeneous solution. Hence we have the following statement.

Proposition 5.6**.**

Let $\hat{n}_{S}$ be given by (5.5). Let $A$ be given by (5.6) and let $\Omega$ be defined as (5.3). Then the sequence of functions $\{n_{k}\}_{k=0}^{\infty}$ defined as

[TABLE]

with

[TABLE]

is the solution to (4.2).

Proof.

The solution to (4.2) is given by the linear combination of the particular solution and of the homogeneous solutions. In particular, the solution must have the following form

[TABLE]

Recall that the solution to (5) are such that $\hat{M}(z)=\hat{n}_{S}(z)+\sum_{k=0}^{\infty}\hat{n}_{k}(z)<\infty$. Since $|\theta_{2}(z)|\geq 1$ for every $z\in\mathcal{C}$ we must take $D(z)=0$ for every $z\in\mathbb{C}$. We now compute $C(z)$. To this end we notice that the equation for $\hat{n}_{0}$ combined with (5.15) implies that

[TABLE]

As a consequence we deduce that

[TABLE]

We then obtain (5.14). ∎

Lemma 5.7**.**

The function $\hat{n}_{k}$ is analytic in $\mathbb{C}\setminus\left(\{z_{S},0\}\cup[z_{1},z_{2}]\right)$ where

[TABLE]

and where we recall that $z_{S}$ is defined in (5.4).

Proof.

Recall that $z_{S}$ is a pole of $\hat{n}_{S}$. Since we have that $A(z_{S})\left(1+B(z_{S})\varphi(z_{S})^{k}\right)\neq 0$, then we deduce that $z_{S}$ is a pole of $\hat{n}_{k}$.

Notice now that the function $A$ has one pole in $z_{A}$ given by (5.7). However notice that $1+B(z_{A})\varphi(z_{A})^{k}=0$. Indeed we have that $\Omega(z_{A})=\alpha(e^{E+\Delta}+1)$ and therefore $\varphi(z_{A})=1$. Moreover we also have that $B(z_{A})=-1$ hence $B(z_{A})\varphi(z_{A})^{k}=-1$. As a consequence we have that $z_{A}$ is not a pole of $\hat{n}_{k}$.

Moreover we can prove that if the function $B$ has a pole, then the pole is in $z_{B}=-e^{\sigma}$. However notice that $\hat{n}_{S}(z_{B})=0$. As a consequence $z_{B}$ is not a pole of $\hat{n}_{k}$. Finally notice that the function $\varphi$ has a branch cut when $z\in\mathbb{R}$ is such that $\Omega(z)^{2}-4\alpha^{2}e^{E+\Delta}\leq 0$. Since we have that $\Omega(z)^{2}-4\alpha^{2}e^{E+\Delta}\leq 0$ if and only if $[z_{1},z_{2}]$ the desired conclusion follows. ∎

We now aim at studying the behaviour of $n_{k}(t)$ for large values of $t$ and large values of $k$. To this end we study the behaviour of $n_{k}(t)$ when $k\approx t$ as $t\rightarrow\infty$. More precisely we will study the behaviour of $n_{\lfloor 2\alpha e^{\frac{E+\Delta}{2}}\theta t\rfloor}(t)$ where $\theta\geq 0$.

Theorem 5.8**.**

*Let us define $\tau:=2\alpha e^{\frac{E+\Delta}{2}}t$ and let $k:=\lfloor\theta\tau\rfloor$ for $\theta\geq 0$.

Assume that $\Delta<\Delta_{c}$. Then*

[TABLE]

Proof.

By the inverse Laplace transform formula we have that

[TABLE]

where $\delta>0$. Moreover recall that since we assume that $\Delta<\Delta_{c}$ we have that $0>z_{A}$. Notice that $n_{k}=\ell_{k}+h_{k}$ where

[TABLE]

and

[TABLE]

Let us compute $\ell_{k}$. To this end notice that

[TABLE]

where $C$ is a closed contour surrounding $z_{A},z_{S}$ and [math]. Since the function $z\mapsto A(z)\hat{n}_{S}(z)$ is a meromorphic function we use Cauchy’s residue theorem and deduce that

[TABLE]

Notice that since $z_{A},z_{S}<0$ the dominant term in $\ell_{k}(t)$ as $t\rightarrow\infty$ is $e^{-kE}A(0)\overline{n}_{S}$. Hence

[TABLE]

Let us now compute $h_{k}$. To this end it is convenient to perform the change of variables

[TABLE]

in the integral (5.19). This change of variables maps [math] in $w_{0}$, $z_{S}$ in $w_{S}$, $z_{A}$ to $w_{A}$, $z_{1}$ to $w_{1}$ and $z_{2}$ to $w_{2}$ where

[TABLE]

and

[TABLE]

Using the fact that

[TABLE]

as well as the fact that $\tau=2\alpha e^{(E+\Delta)/2}t$ we deduce that

[TABLE]

where $\tilde{\delta}=\frac{1}{2\alpha}e^{-\frac{E+\Delta}{2}}\left(\delta+e^{\sigma}+\alpha(e^{E}+e^{\Delta})\right)$, $C(\alpha,E,\Delta)=-4\alpha e^{\Delta/2+E}\sinh(\Delta/2)$, and

[TABLE]

Notice that, as expected, the function $f$ defined as

[TABLE]

is holomorphic in $\mathbb{C}\setminus\left(\{w_{S},w_{0},w_{A}\}\cup[0,1]\right)$. We now compute

[TABLE]

Notice that

[TABLE]

where $\gamma=\cup_{j=1}^{3}\gamma_{j}\cup(\tilde{\delta}-i\infty,\tilde{\delta}+i\infty)$ and where $\gamma_{1}=(-r+i\infty,\tilde{\delta}+i\infty)$, $\gamma_{3}=(-r-i\infty,\tilde{\delta}-i\infty)$ and where $\gamma_{2}$ is a contour that connects $-r-i\infty$ with $-r+i\infty$ passing by $w_{\theta}:=\sqrt{1+\theta^{2}}$ in such a way that the segment $[0,1]$ remains outside $\gamma$. We can compute the first integral. Indeed

[TABLE]

where $\Gamma$ is the set of the poles of $f$ that are in the interior of the closed curve $\gamma$. Moreover, for $j=1,3$ we have that

[TABLE]

Therefore in order to obtain $g_{k}$ we only need to compute the integral of $f$ on $\gamma_{2}$. Notice that by the definition of $f$ we have that, if $k=\theta\tau$, then

[TABLE]

with

[TABLE]

Notice now that the function $\Phi$ has a minimum in $w_{\theta}=\sqrt{1+\theta^{2}}$. Indeed, $\Phi^{\prime}(w_{\theta})=1-\frac{\theta}{\sqrt{w_{\theta}^{2}-1}}=0$. Moreover notice that

[TABLE]

hence $\Phi^{\prime\prime}(w_{\theta})=\frac{\sqrt{1+\theta^{2}}}{\theta^{2}}>0$.

Notice that $\Phi(w_{\theta})=\sqrt{1+\theta^{2}}-\theta\log(w_{\theta}+\theta)\in\mathbb{R}$. We deform the contour $\gamma_{2}$ in such a way that $\gamma_{2}=\gamma_{2}^{\prime}\cup\gamma_{2}^{\prime\prime}\cup\gamma_{2}^{\prime\prime\prime}$ is such that every $w\in\gamma_{2}^{\prime}$ is of the form $w=w_{\theta}+iv$ for some $v\in\mathbb{R}$ and $\gamma_{2}^{\prime\prime}$ and $\gamma_{2}^{\prime\prime\prime}$ connect respectively $\gamma_{2}^{\prime}$ with $-r-i\infty$ and with $-r+i\infty$. Then, using Watson’s Lemma (see [2, Section 6.4]) we deduce that

[TABLE]

By Taylor expansion we have that for $w=\sqrt{1+\theta^{2}}+iv$ it holds that

[TABLE]

As a consequence we deduce that

[TABLE]

We conclude that

[TABLE]

Since we have that $h_{\tau\theta}(\tau)=C(\alpha,E,\Delta)e^{\frac{\theta\tau(\Delta-E)}{2}}e^{-\tau w_{0}}g_{\tau\theta}(\tau)$ we deduce that

[TABLE]

In order to study the behaviour of $n_{\lfloor\tau\theta\rfloor}=\ell_{\lfloor\tau\theta\rfloor}+h_{\lfloor\tau\theta\rfloor}$ we use the fact that $h_{\lfloor\tau\theta\rfloor}\sim h_{\tau\theta}$ as $\tau\to\infty$ and we compare the leading term of $h_{\tau\theta}$ with the leading term of $\ell_{\tau\theta}$ which is

[TABLE]

Now notice that, since $\Delta<\Delta_{c}$, it holds that

[TABLE]

Indeed, we have that

[TABLE]

which implies that

[TABLE]

Therefore $w_{0}>\Phi(w_{0})+\theta\frac{E+\Delta}{2}$ for every $\theta$. Moreover let us define the function $\Psi_{M}$ as

[TABLE]

Notice that this function has a maximum in $\theta_{M}:=\sinh\left(\frac{\Delta+E}{2}\right)$. Moreover we have that $\Psi_{M}(\theta_{M})=\cosh\left(\frac{\Delta+E}{2}\right)$. Now notice that $w_{0}-\cosh\left(\frac{\Delta+E}{2}\right)=e^{\sigma}-\alpha(e^{E}-1)(e^{\Delta}-1)>0$ because $\Delta<\Delta_{c}$. Hence $w_{0}>\Phi(w_{\theta})+\theta\frac{E+\Delta}{2}$ for every $\theta>0$. To prove (5.25) we only need to prove that $w_{0}>\Phi(w_{A})+\theta\frac{E+\Delta}{2}$ for every $\theta>0$. This follows by the fact that

[TABLE]

Hence $\Phi(w_{A})+\theta\frac{E+\Delta}{2}=w_{A}<w_{0}$.

Recall that $\Gamma\subset\{w_{A},w_{S},w_{0}\}$ is the set of the poles of the function $f$ defined as in (5). The fact that equality (5.25) holds implies that the dominant term in $n_{\tau\theta}=\ell_{\tau\theta}+h_{\tau\theta}$ is the leading term of $\ell_{\tau\theta}$ which is (5.24). This implies that (5.18) holds. ∎

Theorem 5.9**.**

Let us define $\tau:=2\alpha e^{\frac{E+\Delta}{2}}t$ and let $k:=\lfloor\theta\tau\rfloor$ for $\theta\geq 0$. Assume that $\Delta>\Delta_{c}$.

If $\theta$ is such that $w_{\theta}<w_{S}$, then we have that

[TABLE] 2. 2.

If instead $\theta$ is such that $w_{\theta}>w_{S}$, then

[TABLE]

Proof.

As in the proof of Theorem 5.8 we apply the inverse Laplace transform formula to deduce that

[TABLE]

where $z_{A}>\delta>0$. As in the proof of Theorem 5.8 $n_{k}=\ell_{k}+h_{k}$ where

[TABLE]

and

[TABLE]

Notice that since $z_{S}<0$ the dominant term in $\ell_{k}(t)$ as $t\rightarrow\infty$ is $e^{-kE}A(0)\overline{n}_{S}$. Hence

[TABLE]

Moreover

[TABLE]

where $g_{k}$ and $C(\alpha,E,\Delta)$ are defined as in the proof of Theorem 5.8. With the same arguments as the ones used in the proof of Theorem 5.8 we deduce that

[TABLE]

where $\gamma$ is as in the proof of Theorem 5.8. In order to compute $\int_{\gamma}f(w)dw$ we need to distinguish between the following three cases.

$w_{0}>w_{\theta}>w_{S}$, 2. 2.

$w_{0}>w_{S}>w_{\theta}$, 3. 3.

$w_{\theta}>w_{0}>w_{S}$.

Case 1: $w_{0}>w_{\theta}>w_{S}$.

Assume that 1. holds, hence that $\theta$ is such that $w_{S}<w_{\theta}<w_{0}$. Then we have that

[TABLE]

Now notice that

[TABLE]

where we are using the notation

[TABLE]

and where the function $a$ is given by (5.21) and the function $\Phi$ by (5.23). Therefore

[TABLE]

Hence

[TABLE]

In order to study the behaviour of $n_{\theta\tau}(t)$ as $\tau\to\infty$ we need to identify the leading order term in $\ell_{\theta\tau}+h_{\theta\tau}$. Since $\Delta>\Delta_{c}$ and $w_{0}>w_{\theta}$, we have that

[TABLE]

Indeed, the fact that $\Delta>\Delta_{c}$ implies that

[TABLE]

Moreover notice that by the definition of $w_{\theta}$ we have that $\Phi(w_{0})>\Phi(w_{\theta})$. Hence the desired conclusion follows. We conclude that in this case

[TABLE]

Therefore using that

[TABLE]

we deduce that (5.26) holds.

Case 2: $w_{0}>w_{S}>w_{\theta}$.

The main difference between case $1$ and case $2$ is the value of the integral $\int_{\gamma}f(w(dw)$. Indeed, when $\theta$ is such that $w_{S}>w_{\theta}$ we have that

[TABLE]

As a consequence, computations similar to the ones made for case 1 imply that

[TABLE]

In order to find the leading order term in $\ell_{\theta\tau}+h_{\theta\tau}$ we prove that

[TABLE]

Indeed, the fact that $\Delta>\Delta_{c}$ implies that

[TABLE]

Therefore the asymptotics (5.26) follow.

Case 3: $w_{\theta}>w_{0}>w_{S}$.

When $\theta$ is such that $w_{\theta}>w_{0}>w_{S}$ we have that

[TABLE]

As a consequence,

[TABLE]

Notice that the fact that $w_{\theta}>w_{0}>w_{S}$ and $\Delta>\Delta_{c}$ imply that

[TABLE]

Hence

[TABLE]

Therefore (5.27) follows.

∎

Theorem 5.10**.**

*Let us define $\tau:=2\alpha e^{\frac{E+\Delta}{2}}t$ and let $k:=\lfloor\theta\tau\rfloor$ for $\theta\geq 0$.

Assume that $\Delta=\Delta_{c}$. Then*

if $\theta$ is such that $w_{\theta}<w_{S}$, then we have that

[TABLE] 2. 2.

If instead $\theta$ is such that $w_{\theta}>w_{S}$, then

[TABLE]

Proof.

Notice that if $\Delta=\Delta_{c}$ we have that $z_{A}=z_{0}$ and as a consequence we have that

[TABLE]

As in the proof of Theorem 5.8 and Theorem 5.9 we deduce that the function $\ell_{k}$ defined as in (5.28) is given by

[TABLE]

In order to obtain the asymptotics of $n_{k}$ we have to compute the asymptotics of the function $h_{\theta\tau}(\tau)=\frac{C(\alpha,E,\Delta)}{2\pi i}e^{-w_{0}\tau}g_{\tau\theta}(\tau)$ where

[TABLE]

where $\gamma$ and $\gamma_{2}$ are defined as in Theorem 5.10 and $f$ is given by (5).

When $w_{0}>w_{\theta}>w_{S}$ then

[TABLE]

where we are using the fact that $\Phi(w_{0})=w_{0}-\theta E.$ This together with detailed asymptotics of $\int_{\gamma_{2}}f(w)dw$ as $\tau\rightarrow\infty$ implies that as $\tau\to\infty$

[TABLE]

Moreover notice that

[TABLE]

and we recall that the maximum of the function $\Psi_{M}=\Phi(w_{\theta})+\theta E$ is attained at $\theta_{M}=\sinh(E)$, hence when $w_{\theta}=\sqrt{1+\theta_{M}^{2}}=\cosh(E)=w_{0}$. Therefore when $w_{0}>w_{\theta}$ we have that

[TABLE]

As a consequence we deduce that

[TABLE]

where we used the fact that $C(\alpha,E,E)=2\alpha(e^{2E}-e^{E})$. When $w_{0}=w_{\theta}>w_{S}$ we obtain the same result. When $w_{0}>w_{S}>w_{\theta}$, then we have that

[TABLE]

and we can prove that (5.34) holds also in this case.

Finally when $w_{\theta}>w_{0}>w_{S}$, then

[TABLE]

This implies (5.33). ∎

6 Analysis of the model with a finite number of kinetic proofreading steps

In this section we study in detail, using Laplace transforms, the model described by (2.1) when the number of kinetic proofreading steps $N$ is large, but finite. We recover the same results obtained using matched asymptotics expansions in Section 4.

Performing the Laplace transform to all the terms in equation (2.1) we obtain the following equation for the vector of the Laplace transforms $(\hat{n}_{S},\hat{n}_{0},\dots,\hat{n}_{N})$ of the concentrations $(n_{S},n_{0},\dots,n_{N})$

[TABLE]

where $\Omega(z)$ is defined as (5.3) and where

[TABLE]

Since the probability of response is given by

[TABLE]

we aim at studying the behaviour of $\hat{n}_{N}(0)$ as $N\to\infty$. Using computations that are similar to the ones of Proposition 5.6 we prove the following Lemma.

Lemma 6.1**.**

Assume that $\hat{n}_{S}(0)\in\mathbb{R}_{*}$ and $\hat{n}_{0}(0)$ are given. Let $\{\hat{n}_{k}(0)\}_{k=1}^{N}$ be the sequence of $\hat{n}_{k}(0)$ defined as

[TABLE]

where

[TABLE]

and where

[TABLE]

The vector $(\hat{n}_{S}(0),\hat{n}_{0}(0),\hat{n}_{1}(0),\dots\hat{n}_{N}(0))$ satisfies

[TABLE]

As a consequence, enforcing that $\hat{n}_{S}(0)$ satisfies $\left(\frac{e^{-NE}-1}{1-e^{-E}}+\mu+e^{\sigma}\right)\hat{n}_{S}(0)=1+e^{\sigma}\hat{M}(0)$ we can prove the following Lemma.

Lemma 6.2**.**

Assume that $\hat{n}_{0}(0)\in\mathbb{R}_{*}$ is given by

[TABLE]

with

[TABLE]

Assume that $\hat{n}_{S}(0)\in\mathbb{R}_{*}$ is given by

[TABLE]

where

[TABLE]

Assume that the vector $(\hat{n}_{0}(0),\hat{n}_{1}(0),\dots,\hat{n}_{N}(0))$ is defined as in Lemma 6.1. Then $(\hat{n}_{S}(0),\hat{n}_{0}(0),\hat{n}_{1}(0),\dots,\hat{n}_{N}(0))$ is the unique solution to

[TABLE]

Lemma 6.3**.**

Let $F_{1}(N)$ and $F_{2}(N)$ be given by (6.3). Then we have

[TABLE]

and

[TABLE]

Proof.

In order to prove this Lemma it is enough to notice that Lemma 5.4 implies that $\varphi(\sigma)<\varphi_{2}$. ∎

Proposition 6.4**.**

Assume that $\hat{n}_{N}(0)$ is as in Lemma 6.2. Then we have that

• •

if $\Delta<\Delta_{c}$, then

[TABLE]

• •

If $\Delta>\Delta_{c}$, then

[TABLE]

where we recall that $G$ is given by (6.2).

Proof.

Lemma 6.3 together with Lemma 6.2 imply that

[TABLE]

Using the asymptotics for $\hat{n}_{1}$ we deduce that

[TABLE]

Hence $\hat{n}_{0}(0)\sim G\hat{n}_{S}$ where

[TABLE]

Using the asymptotics for $\hat{n}_{0}(0)$ as well as the fact that $n_{N}(0)$ is given by (6.2) we obtain that

[TABLE]

The fact that

[TABLE]

implies that

[TABLE]

As a consequence we have that $\hat{n}_{S}(0)\sim\overline{n}_{S}\hat{M}(0)$ as $N\to\infty$. This, combined with the fact that when $\Delta<\Delta_{c}$ we have that $\varphi(\sigma)<1$ while when $\Delta>\Delta_{c}$ we have that $\varphi(\sigma)>1$ to obtain the desired result. ∎

Theorem 6.5**.**

Let $p_{res}(\sigma)$ be the probability of response defined by (2.6). Then we have that

• •

if $\Delta<\Delta_{c}$, then

[TABLE]

where $C_{1}:=\left[\alpha e^{\Delta}A(0)\overline{n}_{S}\left(1-\varphi(\sigma)e^{-(E+\Delta)}\right)\right]^{-1}$.

• •

if instead $\Delta>\Delta_{c}$, then we have that

[TABLE]

where

[TABLE]

Proof.

Recall that

[TABLE]

Plugging $\hat{M}(0)$ in the asymptotics for $\hat{n}_{N}(0)$ and using the fact that $\frac{1}{\varphi_{2}}=\varphi(\sigma)e^{-(E+\Delta)}$, we obtain the desired conclusions. ∎

The constants $C_{1}$ and $C_{2}$ in Theorem 6.5 are the same constants that we obtained using matched asymptotics, i.e. (4.26) and (4.27). (In order to see that the constant $C_{2}$ is the same constant that appears in (4.27) notice that $G-A(0)=BA(0)$ where $B$ is given by (4.9)).

Remark 6.6**.**

Notice that Theorem 6.5 implies that the kinetic proofreading model has strong discrimination properties, in the sense of Definition 2.1, when $\Delta>\Delta_{c}$.

7 Some PDE limits

In this section we derive two PDEs that can be seen as an approximation of the discrete model (2.1) as $N\to\infty$ in different parameters regimes. These PDEs are particular types of renewal equations (see for instance [20]) and can be analysed with computations that are simpler compared to the ones of Section 6.

When $e^{\sigma}\approx\frac{1}{N}$ and $\Delta>E$ and both $E$ and $\Delta$ are of order one, we will derive the following PDE

[TABLE]

for a suitable $\beta>0$, $\delta>0$ and $L>1$ and initial condition

[TABLE]

The steady states of the equation above behave as $e^{-\lambda x}$ where $\lambda=\frac{\beta+\delta}{\alpha(e^{\Delta}-e^{E})}$. Since $\lambda$ depends non trivially on $\beta$ the model has strong discrimination properties. Notice that this is expected as in this regime we have that

[TABLE]

Instead when we assume that $e^{\sigma}\approx\frac{1}{N}$ and $\Delta\approx 1$ and $E\approx\frac{1}{N}$, we derive the following PDE

[TABLE]

If we ignore the loss term $-\delta f$, the steady states of this equation have the form $\overline{m}e^{-x}\left(1+Ce^{(1-\lambda)x}\right)$ for a suitable constants $\overline{m}$ and $C$ and for $\lambda=\frac{\beta}{\alpha(e^{\Delta}-1)}$. Notice that in order to have strong discrimination properties in this case we need to assume that

[TABLE]

This condition is expected, indeed in this regime we have that

[TABLE]

Notice that in order to obtain the PDEs above as an approximation of the discrete system (2.1) we have to assume that $\Delta>E$. The two PDEs that we obtain have a transport term that has a sign, in particular in both the cases we have a transport of complexes toward larger phosphorylation states, i.e. to larger $k$. However we can have situations in which detailed balance fails with $\Delta>\Delta_{c}$, but in which the PDE approximation is not valid because $E>\Delta>\Delta_{c}$. In these cases, the discrete model (2.1) cannot be approximated with a PDE and the transport towards small $k$ dominates the transport toward large $k$.

7.1 PDE limit when $E\approx 1$, $\Delta\approx 1$, $\Delta>E$, $e^{\sigma}\approx\frac{1}{N}$ and $\mu\approx\frac{1}{N}$

In this section we assume that the number of kinetic proofreading steps is $LN$ where $L>1$ and where $N$ tends to infinity. In order to derive equation (7.1) let us consider the behaviour of the solution to (2.1) when $k\approx N$ and when $k\approx 1$ under the assumption that $E\approx 1$, $\Delta\approx 1$, $\Delta>E$ and $e^{\sigma}\approx\frac{1}{N}$, $\mu\approx\frac{1}{N}$.

Behaviour of $n_{k}$ for $k\approx N$.

Let us assume that $e^{\sigma}N=\beta>0$ and that $\mu N=\delta>0$. Notice that (2.1) implies that

[TABLE]

where

[TABLE]

We now define $x=\frac{k}{N}$, $\tau=\frac{t}{N}$, $f(\tau,x)=Nn_{k}(t)$ and $m(\tau)=Nn_{S}(t)$.

We deduce that $f$ and $m$ satisfy the following system of equations

[TABLE]

Using the fact that $e^{\sigma}N=\beta$ and that $\mu N=\delta$ and the fact that $N\left[f\left(\tau,x-\frac{1}{N}\right)-f(\tau,x)\right]\approx\partial_{x}f(\tau,x)$ as $N\to\infty$ as well as the fact that $N\left[f\left(\tau,x+\frac{1}{N}\right)-f(\tau,x)+f\left(\tau,x-\frac{1}{N}\right)\right]\approx\frac{1}{N}\partial^{2}_{xx}f(\tau,x)\approx 0$ as $N\to\infty$ we deduce that

[TABLE]

Now notice that

[TABLE]

Therefore the equation for $m$ can be approximated as

[TABLE]

In particular this implies that $m(\tau)Ne^{-ExN}\rightarrow 0$ as $N\to\infty$ for any $x>0$. Therefore we obtain equation (7.1).

Behaviour of $n_{k}$ for $k\approx 1$.

In order to find the boundary condition at [math] for equation (7.1) we have to match the behaviour obtained for $k\approx N$ with the behaviour in the region $k\approx 1$. In this region we have that

[TABLE]

As in Section 4 we make the quasi steady state approximation, valid when $t\approx N$. As a consequence we deduce that in the region of $k\approx 1$ we have that

[TABLE]

With arguments that are analogous to the ones used in Section 4 and Section 5 we obtain that

[TABLE]

where

[TABLE]

We deduce that $n_{k}\approx B_{2}(n_{S})$ as $k\to\infty$ (when $k\gg 1$ and $N-k\gg 1$), hence

[TABLE]

In order to match the region where $k\approx 1$ with the region where $k\approx N$ we need to have the matching condition

[TABLE]

The boundary condition (7.2) follows.

7.2 PDE limit when $E\approx\frac{1}{N}$, $\Delta\approx 1$, $e^{\sigma}\approx\frac{1}{N}$ and $\mu\approx\frac{1}{N}$

As in Section 7.1, we assume that the number of kinetic proofreading steps is $LN$ where $L>1$ and where $N$ tends to infinity. We consider the behaviour of $n_{k}$ when $k\approx N$ and when $k\approx 1$ under the assumption that $E=\frac{1}{N}$, $\Delta\approx 1$, $e^{\sigma}=\frac{\beta}{N}$, $\mu=\frac{\delta}{N}$.

Behaviour of $n_{k}$ for $k\approx N$.

Defining $x=\frac{k}{N}$, $\tau=\frac{t}{N}$, $f(\tau,x)=Nn_{k}(t)$ and $m(\tau)=N^{2}n_{S}(t)$, we deduce that $f$ and $m$ satisfy the following equations

[TABLE]

Taking $N\to\infty$ in the equation for $m$ we deduce that

[TABLE]

Moreover taking the limit as $N\to\infty$ in the equation for $f$ we deduce that (7.4)-(7.5) holds. In order to obtain the boundary condition at $x=0$ we study the behaviour of $n_{k}$ when $k\approx 1$.

Behaviour of $n_{k}$ for $k\approx 1$.

In this region, since $e^{\sigma}\approx\frac{1}{N}$, $E\approx\frac{1}{N}$ we have that

[TABLE]

Since $t\approx N$ we can assume that $\frac{dn_{k}}{dt}\approx 0$ and that $\frac{dn_{0}}{dt}\approx 0$. We deduce that

[TABLE]

This implies that

[TABLE]

The matching condition with the region where $k\approx N$ implies that

[TABLE]

where we used the fact that $n_{S}\approx\frac{1}{N^{2}}$. This implies that (7.6) holds.

8 Some models without detailed balance and lack of strong discrimination properties

In this section we analyse some examples of kinetic proofreading networks that do not satisfy the property of detailed balance, and with the same amount of detailed balance $\Delta$ as the model analysed in Section 4, and that do not have strong discrimination properties. The networks that we study have the same architecture of the linear kinetic proofreading network described in Section 2.1. We select the chemical rates in such a way that also for these networks equality (2.2) holds for every cycle with $\Delta\neq 0$. In other words, in these networks the Wegscheider conditions associated to the cycles fails exactly in the same manner as in the network analysed in the previous sections. The only difference between the networks that we analyse in this section and the network analysed in the previous sections is the choice of the reaction in which the detailed balance property fails. In other words, the choice of the reaction whose rate is affected by the parameter $\Delta$.

In Subsection 8.4 we will also consider the case of $\Delta\to\infty$. We show that if $\Delta\to\infty$, and $\alpha\to 0$ in such a way that $\alpha e^{\Delta}$ is of order one then we recover the one directional model analysed in [10].

8.1 Detachment rate depending on $\Delta$

We start by introducing a model in which the parameter $\Delta>0$ affects the detachment, i.e. the detachment reactions are accelerated with respect to the case in which $\Delta=0$. We show now that this model does not have strong discrimination properties.

We consider the chemical reaction network generated by the following chemical reactions

[TABLE]

Notice that for every $k$ the set of the reactions in (8.8) forms a cycle. Hence the circuit condition (1.1) implies that the network satisfies the property of detailed balance if and only if $\Delta=0$.

The system of equations associated with the network is the following

[TABLE]

We analyse the behaviour of the steady states $\tilde{n}=(\tilde{n}_{S},\tilde{n}_{0},\tilde{n}_{1},\dots,\tilde{n}_{N})$ of (8.1) as $k\to\infty$, more precisely we show that

[TABLE]

for a suitable constant $C$. Notice that in this case the system does not have strong discrimination properties.

We now indicate with heuristic arguments how the behaviour (8.3) can be obtained. First of all notice that we can always find a particular solution to

[TABLE]

that is such that $\tilde{n}^{(p)}_{S}=1$ and $\tilde{n}_{0}^{(p)}=\tilde{n}_{1}^{(p)}=0$. We now study the homogeneous solutions corresponding to (8.1), i.e. the solutions to the equation

[TABLE]

We use the change of variables $\tilde{n}_{k}^{(h)}=e^{W_{k}}$ and deduce that

[TABLE]

Notice that we have two possibilities, either $W_{k+1}-W_{k}\rightarrow\infty$ as $k\to\infty$, hence $W_{k-1}-W_{k}\rightarrow-\infty$ or $W_{k-1}-W_{k}\rightarrow\infty$ as $k\to\infty$, hence $W_{k+1}-W_{k}\rightarrow-\infty$. In the first case we have that $\alpha e^{E}e^{W_{k+1}-W_{k}}\sim e^{\sigma+k\Delta}$ as $k\to\infty$, therefore $W_{k}\sim\frac{k^{2}}{2}\Delta$ as $k\to\infty$. When instead we assume that $W_{k-1}-W_{k}\rightarrow\infty$ we deduce that $\alpha e^{E}e^{W_{k-1}-W_{k}}\sim e^{\sigma+k\Delta}$ as $k\to\infty$, hence $W_{k}\sim-\frac{k^{2}}{2}\Delta$ as $k\to\infty$.

Let $\psi_{k}^{+},\psi_{k}^{-}$ be two solutions to (8.4) such that $\psi_{k}^{+}\sim e^{\frac{k^{2}}{2}\Delta}$ as $k\to\infty$ and $\psi_{k}^{-}\sim e^{-\frac{k^{2}}{2}\Delta}$ as $k\to\infty$. Let us define $\tilde{n}_{k}$ as

[TABLE]

Without loss of generality we can assume that $\psi_{0}^{+}=\psi_{1}^{+}=1$, so that $\psi^{+}_{k}$ is unique, and we select $a,b,\psi_{0}^{-},\psi_{1}^{+}$ in such a way that

[TABLE]

where

[TABLE]

This implies that $n_{k}^{(p)}\sim-a\frac{\psi_{k}^{+}}{\tilde{n}_{S}}$ as $k\to\infty$. More precisely,

[TABLE]

Therefore (8.3) holds.

8.2 Attachment rate depending on $\Delta$

We consider a model in which the parameter $\Delta>0$ affects the attachment rate, hence the attachment reactions are accelerated. We show that this model can have strong discrimination properties under certain assumptions on the parameters.

We study the chemical reaction network generated by the following chemical reactions

[TABLE]

The system of ODEs corresponding to this model is the following

[TABLE]

We analyse the steady states $\tilde{n}=(\tilde{n}_{S},\tilde{n}_{k})_{k=0}^{N}$ of the system of ODEs (8.2). Using arguments that are analogues to the ones used in Section 4 we deduce that the solutions of the form

[TABLE]

where

[TABLE]

and where $B$ and $C$ are suitable constants. Moreover, the function $f$ defined as

[TABLE]

is such that $f(0)<1$ while $\lim_{\Delta\to\infty}f(\Delta)>1$. As a consequence there exists a $\Delta_{c}$ such that if $\Delta<\Delta_{c}$ then the model does not have strong discrimination properties while if $\Delta>\Delta_{c}$ the model has strong discrimination properties, i.e. the steady state is given by $e^{-k\lambda(\sigma,E)}$ where $\lambda(\sigma,E)$ depends on $\sigma$.

8.3 Dephosphorylation rate depending on $\Delta$

We assume now that $\Delta$ affects the dephosphorylation reaction. In this case we can obtain strong discrimination only if $\alpha$, the phosphorylation rate, is sufficiently large. Consider the chemical reaction network generated by the following chemical reactions

[TABLE]

We formulate the system of equations

[TABLE]

The steady states $\tilde{n}=(\tilde{n}_{S},\tilde{n}_{k})_{k=0}^{N}$ of the system of ODEs (8.3) are given by

[TABLE]

where $B,C$ are suitable constants depending on the parameters and where

[TABLE]

In this case we have that if $\frac{e^{\sigma-E}}{\alpha(1-e^{-E})}<1$, then there exists a critical $\Delta_{c}$

[TABLE]

such that the asymptotic behaviour of $\tilde{n}_{k}$ as $k\to\infty$ depends on $\sigma$ only when $\Delta>\Delta_{c}$, i.e. $\tilde{n}_{k}\sim e^{-k\lambda(\sigma,E)}$ as $k\to\infty$ where $\lambda(\sigma,E)$ depends on $\sigma$ only when $\Delta>\Delta_{c}$.

8.4 $\Delta\to\infty$ and $\alpha e^{\Delta}\approx 1$

In this section we derive the one directional model that we analysed in [10] as a limit case of the model studied in this paper.

We consider the model with infinite states, i.e. we assume that the model is described by (5). We now study the steady states, that satisfy (4.3)-(4.4). As explained in Section 5 the behaviour of the steady states is expected to be the following one

[TABLE]

where $B$ is given by (4.9) and where $\varphi(\sigma)$ is given by (4.7). We assume that $\alpha e^{\Delta}=\gamma\approx 1$, then

[TABLE]

Since we have that

[TABLE]

we deduce that $\tilde{n}_{k}\sim e^{-\lambda(\sigma,E)k}$ where $\lambda(\sigma,E)=E$ if $\log\left(1+\frac{e^{\sigma}}{\gamma}\right)>E$ while $\lambda(\sigma,E)=\log\left(1+\frac{e^{\sigma}}{\gamma}\right)$ when $\log\left(1+\frac{e^{\sigma}}{\gamma}\right)<E$. In the second case we have strong discrimination properties because $\lambda(\sigma,E)$ depends on $\sigma$.

In order to understand the condition for discrimination $\log\left(1+\frac{e^{\sigma}}{\gamma}\right)<E$, recall that

[TABLE]

Assume now that $\Delta\to\infty$ and that $\Delta-\Delta_{c}=x>0$. Under this assumption, we expect the limiting model to have strong discrimination properties. This assumption implies that the parameters satisfy the assumption $\log\left(1+\frac{e^{\sigma}}{\gamma}\right)<E$, indeed

[TABLE]

In the model that we consider in [10] we have that $E=\infty$. Indeed in that paper we assume that when the ligands attach to the receptors then they form a complex that has always state $0.$ In that case is always true that $\log\left(1+\frac{e^{\sigma}}{\gamma}\right)<E=\infty$, hence we always have discrimination.

9 Conclusions and open problems

In this paper we prove that a minimal amount of lack of detailed balance is necessary in order to have strong discrimination properties. In other words we prove that the kinetic proofreading systems must be out of equilibrium systems. This absence of equilibrium can be measured in terms of the frozen concentrations of some substances, like ATP or ADP, that must be sufficiently far from equilibrium.

Given that the lack of detailed balance is a property on the chemical rates of the reactions that belong to the cycles, one could think that similar amounts of lack of detailed balance would result in a similar behaviour of the network. We prove that this is not the case. We give examples of models that have identical amount of lack of detailed balance and that do not have the same discrimination properties.

We also have shown that, in some parameters regimes, the chemical network studied in this paper can be approximated by means of PDEs for which the mathematical analysis of the discrimination properties is more amenable.

It could be interesting to study the behaviour of the network when the lack of detailed balance $\Delta_{k}$ in the phosphorylation reaction $(k)\rightarrow(k+1)$ depends on $k$. It would be relevant to understand if also in this case a critical amount of detailed balance is needed in order to have strong discrimination.

We could also consider more complex models in which the phosphorylation reaction is modelled in more detail. In particular we can assume that phosphorylation is a chain of reactions that involves $ATP$, $ADP$ and phosphate groups, but also other substances that are out of equilibrium, for instance some enzymes. It would be relevant to understand if including these substances in the system, hence changing the topology of the network, can improve the discrimination properties of the chemical network.

In this paper we analyse in detail the discrimination property of the kinetic proofreading mechanism and its relation with detailed balance. Understanding if a biological function can be achieved by passive systems (with detailed balance) or if it requires active processes is a relevant question for several biological systems. For instance, in [8] we analyse whether the property of adaptation, that takes place in many relevant signalling systems, can be achieved or not by equilibrium chemical reaction networks.

Acknowledgements The authors gratefully acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG) through the collaborative research centre ”The mathematics of emerging effects” (CRC 1060, Project-ID 211504053) and Germany’s Excellence Strategy EXC2047/1-390685813. The funders had no role in study design, analysis, decision to publish, or preparation of the manuscript.

E. Franco: Institute for Applied Mathematics, University of Bonn,

Endenicher Allee 60, D-53115 Bonn, Germany

E-mail: franco@iam.uni-bonn.de

ORCID 0000-0002-5311-2124

J. J. L. Velázquez: Institute for Applied Mathematics, University of Bonn,

Endenicher Allee 60, D-53115 Bonn, Germany

E-mail: velazquez@iam.uni-bonn.de