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 the states in the network a Gibbs free energy (that will be written in $k_B T$ units). Without loss of generality, we assume that the free energy $E_S$ of free ligands is the basal energy, and we normalize it to be equal to zero. We then assume that the energy difference produced by each phosphorylation event is $E \gt 0$ . We denote with $\sigma \in \mathbb R$ the binding energy between the receptor and the ligand. This allows us to assume that the Gibbs free energy associated to a complex ligand-receptor at phosphorylation state $j \in \{ 0, \ldots , N \}$ (i.e. a complex ligand receptor with attached $j$ phosphate groups) is $E_j \,:\!=\, \sigma + j E$ .

As we will explain below, the assumptions we made on the Gibbs free energies allow us to write the chemical rates associated to the chemical reactions as follows. We assume that the attachment of a ligand to the receptor, i.e. the reactions of the form

\begin{equation*} {R_1^{(j)} \,:\, S \overset {e^{ -j E } }{\longrightarrow } C_j} \end{equation*}

take place at rate $k_{R_1^{(j)}} =e^{- j E }$ for every $j \in \{ 0, \ldots , N \}$ . The detachment rate depends on the binding energy between the ligand and the receptor, i.e. the reactions of the form

\begin{equation*} {-R_1^{(j)}\,:\, C_{j} \overset {e^{\sigma } }{\longrightarrow } S } \end{equation*}

take place at rate $k_{- R_1^{(j)}} = e^{\sigma }$ for every $j \in \{0, \ldots , N \}$ where $\sigma \in \mathbb R$ . The reactions

\begin{equation*} {R^{(j)}_2\,:\, C_j \overset {\alpha e^\Delta }{\longrightarrow } C_{j+1} } \end{equation*}

take place at rate $k_{R_2^{(j)}}= \alpha e^{\Delta }$ for every $j \in \{ 0, \ldots , N \}$ . On the other hand, we assume that the reaction

\begin{equation*} {-R^{(j)}_2\,:\, C_j \overset {\alpha e^{E } }{\longrightarrow } C_{j-1} } \end{equation*}

takes place at rate $k_{-R^{(j)}_2}=\alpha e^{E_j- E_{j-1} }= \alpha e^E$ for every $j \in \{0, \ldots , N \}$ .

Figure 1.

This is the sketch of the chemical reactions considered in this paper and the corresponding chemical rates.

See Figure 1 for a visual representation of the model. We now explain our choice of chemical rates and their interpretation. To this end, it is convenient to assume that $\Delta =0$ and that $b =0$ . Indeed, when $\Delta =0$ , we have that

\begin{equation*} 1= \frac { k_{R_1^{(j)}}k_{R_2^{(j)}} k_{-R_1^{(j+1)}} }{ k_{- R_1^{(j)} }k_{- R_2^{(j)} }k_{ R_1^{(j +1)} }} \end{equation*}

for every cycle

\begin{equation*} S \overset {R_1^{(j)}} \longrightarrow C_j \overset {R_2^{(j)}}\longrightarrow C_{j+1} \overset {- R_1^{(+1)}}\longrightarrow S. \end{equation*}

The Wegscheider criterium then guarantees that detailed balance holds when $\Delta =0$ . It is well known that in this case the ratio between the chemical rate of the direct reaction and the chemical rate of the reverse reaction can be written in terms of the Gibbs free energies. More precisely, for the specific system under consideration, we have that

(2.1) \begin{equation} \frac {k_{R_1^{(j)}} }{k_{- R_1^{(j)}}} = e^{ - E_j } , \quad \frac {k_{R_2^{(j)}} }{k_{- R_2^{(j)}}} = e^{ E_{j} - E_{j+1} } =e^{-E} , \quad j \in \{ 0, \ldots , N \}. \end{equation}

Notice that our choice of chemical rates satisfies (2.1) when $\Delta =0$ . The parameter $\alpha$ is the constant of proportionality between $k_{ R_2^{(j)}}$ and $k_{- R_2^{(j)}} e^{ E_j - E_{j+1} }$ . Finally, since we want to study a system with lack of detailed balance and since we assume that the energy is consumed during the phosphorylation reactions $R_2^{\kern1pt j}$ , we define the chemical rates $k_{R_2}^{(j) }$ as $e^\Delta k^{DB}_{{R_2}^{(j)}}$ , where $ k^{DB}_{{R_2}^{(j)}}$ is the rate of the phosphorylation reaction when detailed balance holds, hence $ k^{DB}_{{R_2}^{(j)}} =\alpha$ .

We call the parameter $\sigma$ that determines the detachment rate, binding energy because it is the energy associated with the dephosphorylated complex ligand-receptor $C_0$ . In the chemical literature, the detachment rate is usually written as $ e^{ - E_A}$ where $E_A$ is the activation energy of the unbinding reaction. Due to the normalization $E_S =0$ that we make in our model and the assumption that the rate of the reaction $S \rightarrow C_0$ is equal to $1$ , it turns out that, for the specific model under consideration, the energy $\sigma$ completely determines the detachment rate, which is $e^\sigma$ .

The assumptions above on the chemical rates guarantee that for every cycle $\omega _k=[R_1^{(j)}, R_2^{(j)}, -R_1^{(j+1)}]$ of the form

\begin{equation*} S \overset {R_1^{(j)}} \longrightarrow C_j \overset {R_2^{(j)}}\longrightarrow C_{j+1} \overset {- R_1^{(j+1)}}\longrightarrow S \end{equation*}

it holds that

(2.2) \begin{equation} {e^{\Delta } = \frac { k_{R_1^{(j)}}k_{R_2^{(j)}} k_{-R_1^{(j+1)}} }{ k_{- R_1^{(j)} }k_{- R_2^{(j)} }k_{ R_1^{(j +1)} }}. } \end{equation}

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 $C_j \rightarrow C_{j+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.

Finally, we assume that the ligand can jump to the absorbing state $\emptyset$ from any other state. Once the ligand reaches the absorbing state $\emptyset$ , it cannot interact again with the receptor. The rationale behind this is that we assume that the pMHC and the T-cells interact for a characteristic time that we call signalling lifetime, and we denote with $\mu ^{-1}$ . We stress that the fact that the ligand reaches the state $\emptyset$ does not necessarily mean that the ligand is degraded. The ligand can reach the state $\emptyset$ also, for instance, due to the fact that the pMHC stops interacting with the T-cells or due to the fact that the activation time of the receptors is limited due to any kind of mechanism, including ligand degradation, T-cell motion, or the detachment of the peptide from the MHC.

The existence of a characteristic interaction time between the ligands and the receptors is crucial in kinetic proofreading systems. If we assume that the interaction time between the T-cells and ligands is unbounded (and hence $\mu =0$ ), then all the ligands will eventually produce a response, even if the response would require a very long time. Therefore, the strong specificity property cannot hold when $\mu =0$ .

Moreover, we assume that the signalling lifetime $\mu ^{-1}$ is such that

(2.3) \begin{equation} \mu \sim e^{- b N } \text{ as } N \to \infty \end{equation}

where $b \gt 0$ . This means that the characteristic time at which the interaction between the T-cells and the ligands takes place scales as $T \sim e^{bN }$ as $N \to \infty$ . This is a mathematical assumption that we make to get a precise scaling limit in which we can prove that strong specificity properties hold. We do not claim that there exists a physical dependence between the characteristic interaction time and the number of kinetic proofreading steps, but the numerical relation between $\mu$ and $N$ allows to obtain a mathematical limit in which the strong specificity property can be defined in precise mathematical terms.

Let $n_k (t)$ be the probability that a ligand, characterised by the binding energy $\sigma$ , reaches the state $k \in \{ 0, \ldots , 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), \ldots , n_N(t) )^T \in \mathbb R_*^{N+1}$ is described by the following system of ODEs

(2.4)

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

\begin{equation*} S_N(E)\,:\!=\,\sum _{k=1}^{N} e^{- k E } = \frac {1-e^{-NE}}{e^E-1 } . \end{equation*}

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

\begin{equation*} {R(t) \,:\!=\, \alpha e^\Delta \int _0^t n_N(s ) ds .} \end{equation*}

The probability $p_{deg} (t)$ that a ligand stops interacting with the receptor before producing a response in the time interval $(0, t]$ is defined as

\begin{equation*} {p_{deg} (t) \,:\!=\, \mu \int _0^t M(s) ds } \end{equation*}

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

(2.5) \begin{equation} M(t)\,:\!=\,n_S(t) +\sum _{k=0}^{N} n_k (t) . \end{equation}

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

(2.6) \begin{equation} \frac {d M }{ dt } = - \alpha e^\Delta n_{N} - \mu M . \end{equation}

Finally, notice that, as expected, we have that

\begin{equation*} p_{deg}(t) + R(t) + M(t)=1, \quad \forall t \geq 0. \end{equation*}

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

(2.7) \begin{equation} p_{res} (\sigma ) \,:\!=\, \int _0^\infty R(t) dt = \alpha e^\Delta \int _0^\infty n_{N} (t) dt \end{equation}

as the number of kinetic proofreading steps $N \to \infty$ . In other words, $p_{res} (\sigma )$ is the probability that the lifetime of the ligand is longer than the time to reach the phosphorylation state $C_N$ (when both times are exponentially distributed).

2.2. Strong discrimination properties for $\Delta \gt \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.4) has strong discrimination properties with respect to the energy $\sigma$ if

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

where $ \lambda _\Delta (\sigma , E )$ is a parameter that depends non trivially on $\sigma$ and there exists a solution $\sigma _c$ to

(2.9) \begin{equation} \lambda _\Delta (\sigma _c, E)= b. \end{equation}

We recall that the parameter $b$ here characterizes the ligand signalling lifetime as the number of kinetic proofreading steps $N$ goes to infinity (see 2.3).

In this paper, we prove that for the kinetic proofreading networks introduced in Section 2.1, we have that the probability of response always satisfies (2.8), but $\lambda _\Delta (\sigma , E)$ depends non trivially on $\sigma$ only when $\Delta \gt \Delta _c(\sigma )$ ; indeed, we will prove that

(2.10) \begin{equation} \lambda _\Delta (\sigma , E )=\begin{cases} E &\text{ if } \Delta \lt \Delta _c (\sigma ) \\ \psi _\Delta (\sigma , E) &\text{ if } \Delta \gt \Delta _c(\sigma ) \end{cases} \end{equation}

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

(2.11) \begin{equation} \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 ] \end{equation}

and where

(2.12) \begin{equation} \Delta _c= \Delta _c (\sigma ) \,:\!=\, \log \left ( 1+ \frac {e^\sigma }{\alpha (e^E-1 ) } \right ). \end{equation}

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 \lt \Delta _c(\sigma )$ and depends on $\sigma$ for $\Delta$ and $\sigma$ such that $\Delta \gt \Delta _c(\sigma )$ is the most important result of this paper. It implies that the kinetic proofreading networks described by (2.4) can have strong discrimination properties, in the sense of Definition 2.1, only when $\Delta$ and $\sigma$ are such that $\Delta \gt \Delta _c(\sigma )$ .

2.3. Necessary conditions on $b$ for strong discrimination when $\Delta \gt \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.9) 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 $0$ around a critical binding energy $\sigma _c$ . Indeed, (2.8) implies that the probability $p_{res} (\sigma )$ that one ligand characterized by the binding energy $\sigma$ satisfies

(2.13) \begin{equation} p_{res} (\sigma ) \sim \left ( 1 + Ce^{N(\lambda (\sigma , E)- b)}\right )^{-1} \text{ as } N \to \infty \end{equation}

where $C$ is a positive constant that depends on the parameters. When $\Delta \lt \Delta _c(\sigma )$ , we have that $\lambda (\sigma , E )=E$ . As a consequence, if $\lambda (\sigma , E )=E \lt b$ , then $\lim _{N\to \infty } p_{res} (\sigma ) =1$ , while when $\lambda (\sigma , E )=E \gt 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 \gt \Delta _c = \Delta _c( \sigma _c)$ . By the continuity of the function $\Delta _c(\sigma )$ defined by (2.12), there exists a neighbourhood $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\gt 0$ such that (2.9) holds. Precise conditions that the parameter $b$ must satisfy are given in (4.28). In particular, we must have that $b\lt 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 2, we must have $\psi (\sigma _c, E ) = b \lt E$ .

Figure 2.

In this figure, we plot the line $\psi (\sigma , E)=E$ for $\alpha =1$ , $E=\ln (2)$ . Notice that $\psi (\sigma , E )=E$ if and only if $\Delta = \Delta _c(\sigma )$ . The line $\psi (\sigma , E)=E$ separates the region in which we have strong discrimination properties from the region where we do not have strong discrimination.

The condition $b \lt E$ , needed to have strong discrimination, has an interpretation: the characteristic lifetime 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.9) 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

\begin{equation*} |\lambda (\sigma , E ) - \lambda (\sigma _c, E ) |= |\lambda (\sigma , E) -b | \lesssim \frac {1}{N}. \end{equation*}

Moreover, we have that if $\sigma \gt \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 \lt \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.9), 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 5, we compare the probability of response $p_{res} (\sigma )$ that we obtain selecting the parameters in such a way that $\Delta \gt \Delta _c$ with the probability of response when $\Delta \lt \Delta _c$ . Notice that in the latter case, the probability is almost constantly equal to $0$ , while in the first case, we have a sharp transition from response to non response around a critical binding energy $\sigma _c$ .

Due to the fact that $\Delta _c(\sigma ) \rightarrow 0$ as $ \sigma \to - \infty$ , it is not possible to prove the existence of a minimal lack of detailed balance in order to have strong specificity for all the binding energies $\sigma$ . However, it is possible to define a critical lack of detailed balance $\overline \Delta _c$ if the order of magnitude of $\sigma$ is prescribed. More precisely, suppose that we want to study the ability of a kinetic system to distinguish between ligands characterized by binding energies that are in the interval $I\,:\!=\, [\sigma _m , \sigma _M ]$ . Let us define $\overline \Delta _c \,:\!=\, \max _{\sigma \in I } \Delta _c(\sigma )$ . If $\Delta \gt \overline \Delta _c$ , then (2.10) implies that $\lambda _\Delta (\sigma , E)$ depends non trivially on $\sigma$ on the interval $ I$ . Assume that $b= \lambda _\Delta (\sigma _c, E )$ for some $\sigma _c \in I$ . As a consequence, since $p_{res}(\sigma )$ satisfies (2.8), we deduce that if $\Delta \gt \overline \Delta _c$ then the strong discrimination property holds. In particular, we will see a sharp transition from $1$ to $0$ in the probability of transition $p_{res} (\sigma )$ in region of order $1/N$ around the critical binding energy $\sigma _c$ .

If, instead, we define $ \underline { \Delta _c} \,:\!=\, \min _{\sigma \in I } \Delta _c(\sigma )$ and consider $\Delta \lt \underline { \Delta _c}$ , then the system does not have strong discrimination properties as $\lambda _\Delta (\sigma , E)= E$ for every $\sigma \in I$ . Hence, depending on whether $b \gt E$ or $E \gt b$ , we will have either that $p_{res}(\sigma ) \approx 1$ for every $\sigma \in I$ or that $p_{res}(\sigma ) \approx 0$ for every $\sigma \in I$ .

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

\begin{equation*} \lim _{\Delta \to \infty } \psi (\sigma , E )= E. \end{equation*}

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.

2.4. Some numerical solutions

In Figure 3, we compare the probability of response that we compute numerically when $N=3$ and when the lack of detailed balance is very small with the probability of response that we compute numerically for $N = 15$ for the same parameters. We notice that when $N =3$ , we have a transition of the probability of response from a positive value to $0$ . In this case, we have that when $\sigma$ is small, the probability of response does not reach the value $1$ , but remains between $1$ and $0.8$ . The reason is that when $N \approx 1$ , the parameter $\mu$ is not negligible; hence, even when $\sigma$ is small, the probability of reaching the state $\emptyset$ is positive. Notice that the situation changes when we select $N = 15$ . In this case, the probability of response approaches the value one for small binding energies.

On the other hand, we notice that in both cases ( $N=3$ and $N=15$ ), the transition in the probability of response takes place in a region that has thickness of order $1$ . In particular, the thickness of the region of transition does not depend on $N$ and remains of order $1$ as $N$ increases.

Figure 3.

The continuous line in blue is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b$ is selected in order to have the transition in $\sigma =1$ , $\Delta =0.01$ and $N=15$ . The dashed red line is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b= \log (2)$ , $\Delta =0.01$ and $N=3$ .

In Figure 4, we compare the probability of response that we compute numerically when $N=3$ and when the lack of detailed balance is large with the probability of response that we compute numerically for $N=7$ , $N = 15$ , and $N=45$ for the same parameters. As before, when $N =3$ and $N=7$ , the probability of response does not approach the value $1$ for small $\sigma$ . This is due to the fact that, also in this case, we have that when $N \approx 1$ the parameter $\mu$ is of order $1$ . However, in contrast with the situation in Figure 3, the probability of response approaches a step function as $N \to \infty$ . In particular, for the parameters considered in Figure 4, we have that the thickness of the region of transition from response to non response decreases as $N$ increases. This is the mechanism of increasing specificity for larger $N$ that we are studying in this paper.

Figure 4.

The dotted line in blue is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b$ is selected in order to have the transition in $\sigma =1$ , $\Delta =2$ and $N=3$ . The violet line with dots and dashes “ $\cdot - \cdot$ ” is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b$ is selected in order to have the transition in $\sigma =1$ , $\Delta =2$ and $N=7$ . The dashed red line is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b$ is selected in order to have the transition in $\sigma =1$ , $\Delta =2$ and $N=15$ . The continuous line in yellow is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $\Delta =4$ , $b$ is selected in order to have the transition in $\sigma =1$ , and $N=45$ .

Finally, in Figure 5, we plot the asymptotics of the probability of response (2.13) for different parameter regimes. In particular, we compare the asymptotics obtained when $\Delta$ is subcritical with the asymptotics when $\Delta$ is supercritical. For the selected interval of binding energies $\sigma \in [1, 3]$ , we have a sharp transition from response to non response when $\Delta$ is supercritical; we do not have a transition when $\Delta$ is subcritical.

Figure 5.

The dashed line in orange is the probability of response $p_{resp} (\sigma )$ when $\alpha =1$ , $E=\log (3)$ , $b= \log (2)$ , $\Delta =1/10$ and $N=20$ . In this regime, we have that $\Delta$ and $\sigma$ are such that $\Delta \lt \Delta _c(\sigma )$ , hence we are in the subcritical case and since $b \lt E$ we have $p_{res} \approx 0$ . The continuous line in blue is a plot of the probability of response when $\alpha , E, b , N$ are as before, but $\Delta =2$ . In this regime $\Delta$ and $\sigma$ are such that $\Delta \gt \Delta _c(\sigma )$ , hence we are in the supercritical case.

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 the molecules of $ATP$ , $ADP$ and the phosphate groups explicitly. 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

(3.1) \begin{equation} S \overset {e^{- k E} }{\underset {e^\sigma } \rightleftarrows } C_k \overset {\alpha e^\Delta }{\underset {\alpha e^E }\rightleftarrows } C_{k+1} \overset {e^\sigma }{\underset {e^{- (k+1)E}} \rightleftarrows } S. \end{equation}

Notice that if $\Delta =0$ , then we have that the Wegscheider condition (1.2) 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

(3.2) \begin{equation} C_k + (ATP) \overset {R_1^{(k)}}{\underset {- R_{1}^{(k)}} \rightleftarrows } C_{ k +1} + (ADP), \quad (ATP) \overset {R_2}{\underset {-R_{2}} \rightleftarrows } (ADP) +(P), \quad C_k \overset {R_3^{(k)}}{\underset {-R_{3}^{(k)}} \rightleftarrows } (S) + C_k * (P), \end{equation}

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

\begin{equation*} C_k *(P) = \overbrace {(P) + (P) + \ldots + (P)}^{ k \text{ times} } \end{equation*}

We assume that the chemical rates of the reactions $ \{ R_1^{(k)},- R_{1}^{(k)} \}_{k=0}^{N}$ are independent on $k =0, \ldots 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

\begin{equation*} K_{R_1^{(k)}}=K_1 =\alpha e^{E_T} , \ K_{-R_1^{(k)}}=K_{-1 }= \alpha e^{ E+ E_D }, \quad K_{R_2}=e^{E_T},\ K_{-R_2} = e^{ E_D+E_P }, \end{equation*}

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

\begin{equation*} K_{R_3^{(k)}}= K_3=e^\sigma ,\ K_{-R_3^{(k)}}= K_{-3}= e^{ k (E_P- E)}, \quad k \in \{ 0, \ldots , N\}. \end{equation*}

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

\begin{align*} (0,0,0,0,0, 0) & \overset {R_1^{(k)}} \rightarrow ({-}1,1,-1,1,0, 0) \overset { R_{-2}} \rightarrow ({-}1,1,0,0,-1, 0) \overset {R_{-3}^{(k)}}\rightarrow (0,1,0,0,-(k+1), -1 )\\& \overset {R_{3}^{(k+1)}} \rightarrow (0,0,0,0,0, 0). \end{align*}

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, \ldots , N$ are cycles of the enlarged chemical reaction network. Notice that by construction, we have that the circuit condition (1.2) holds along every cycle $\overline \omega _k$ ; indeed, the rates were selected in such a way that

\begin{equation*} K_1 K_{-2 } K_{- 3 }^{(k)} K_3 = K_{-1} K_{2 } K_{ 3 } K_{-3}^{(k+1)}, \quad \forall k \in \{0, \ldots , N\} . \end{equation*}

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:

(3.3)

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

\begin{equation*} \textbf {E} =(\sigma , \sigma + E, \sigma +2E, \ldots , \sigma + (N-1) E, \sigma + N E, E_T , E_D, E_P , 1 ) \end{equation*}

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

Moreover, notice that we have some conserved quantities in our system. In particular, we have that if $n(t)=(n_0, n_1, \ldots , 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

\begin{equation*} m_1^T n (t) = m_1^T n(0), \quad m_2^T n (t) = m_2^T n(0) \end{equation*}

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

\begin{equation*} n_T(t) +2 n_{D} (t)+ n_P (t)+ \sum _{k=0}^N k n_k(t) = n_T (0)+2 n_{D} (0) + n_P (0) + \sum _{k=0}^N k n_k (0) \end{equation*}

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

\begin{equation*} \sum _{k=0}^N n_{k} (t) + n_S (t) = \sum _{k=0}^N n_{k} (0) + n_S (0). \end{equation*}

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.3) satisfies the property of detailed balance, all the steady states of (3.3) are of the form

\begin{equation*} N = e^{- \textbf {E} + \mu _1 m_1 + \mu _2 m_2 } \end{equation*}

for suitable constants $\mu _1, \mu _2$ , see Lemma 3.10 in [Reference Franco and Velázquez10].

Therefore, the concentrations of $ATP$ , $ADP$ and $P$ at any steady state $N \in \mathbb R_+^{N +4 }$ of (3.3), denoted respectively by $N_{T}, N_D, N_p$ , satisfy the following equality

\begin{equation*} \frac {N_T }{N_P N_D } = e^{ E_D + E_P -E_T }. \end{equation*}

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.3) is in contact with a reservoir 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

(3.4) \begin{equation} \Delta = E_T - E_D - E_P + \log \left ( \frac {\overline n_T }{ \overline n_P \overline n_D } \right ). \end{equation}

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

(3.5) \begin{equation} { C_k \overset {\overline R_1^{(k)} }{\underset {\overline R_{-1}^{(k)} } \rightleftarrows }C_{k+1}, \quad C_k \overset { \overline R_3^{(k)} }{\underset {\overline R_{-3}^{(k)} } \rightleftarrows } (S) \quad k \in \{ 0, \ldots , N\} } \end{equation}

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

\begin{equation*} \overline {K}_1 = \alpha e^{E_T}\overline n_T , \quad \overline K_{-1}= \alpha e^{E+ E_D }\overline n_D, \quad \overline K_3= e^\sigma , \quad \overline K_{-3}^{(k)} = e^{k (E_P- E) } \overline n_P^k. \end{equation*}

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

\begin{equation*} \overline n_P = e^{- E_P} , \quad \overline n_D=e^{- E_D}, \quad \overline n_T = e^\Delta e^{- E_T} \end{equation*}

then the chemical reaction network (3.5) 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.5) and the chemical reaction network (3.1) coincide.

Notice that the concentrations 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.3), taking into account these fluxes. The equations that we obtain are the following:

(3.6) \begin{align} \frac {d n_T }{dt } &= - n_T \left (e^{E_T} + \alpha e^{E_T} \sum _{k=0}^{N-1} n_k \right ) + \alpha e^{ E+ E_D } n_D \sum _{k=1}^{N}n_{k} + e^{ E_D+E_P } n_P n_D + J^{ext}_T(n) \nonumber\\[3pt] \frac {d n_D }{dt } &= - n_D \left ( e^{ E_D+E_P } n_P + \alpha e^{ E+ E_D } \sum _{k=1}^{N} n_{k} \right ) + \alpha e^{E_T} n_T \sum _{k=0}^{N-1}n_{k} + e^{E_T} n_T + J^{ext}_D(n) \nonumber \\[3pt] \frac {d n_P }{dt } &= - e^{ E_D+E_P } n_P n_D - n_S \sum _{k=0}^{N} k e^{k (E_p - E)} n_P^k + e^{E_T} n_T + e^\sigma \sum _{k=0}^{N} k n_k + J^{ext}_P(n). \end{align}

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

\begin{align*} J^{ext}_T(N) &= e^{\Delta } -1 + \alpha e^{- \sigma } \left ( \sum _{k= 0}^{N-1}e^{- k E} - e^\Delta \sum _{k= 1}^{N}e^{- k E} \right ) = \left ( e^{\Delta } -1 \right ) \left ( 1 + \alpha e^{- \sigma } \frac {1-e^{- NE} }{ 1-e^{-E}} \right ) \gt 0 \\[3pt] J^{ext}_D(N) &= 1 - e^{\Delta } + \alpha e^{- \sigma } e^E \sum _{k=1}^{N} e^{- k E} - e^\Delta \sum _{k=0}^{N-1} e^{- k E} = (1 - e^{\Delta } )\left ( 1 + \alpha e^{- \sigma } \frac {1-e^{- NE} }{ 1-e^{-E}} \right ) \lt 0 \\[3pt] J^{ext}_P(N) &= 1 - e^\Delta + n_S \sum _{k=0}^{N} k e^{- k E} -\sum _{k=0}^{N} k e^{- k E} = 1- e^\Delta \lt 0. \end{align*}

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 maintain the system at steady state.

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.4), i.e.

(4.2) \begin{align} \frac {d n_0^{(out)} }{dt} &= n_S - e^{\sigma } n_0^{(out)} + \alpha e^E n_1^{(out)} - \alpha e^{\Delta } n_0^{(out)} \nonumber \\[3pt] \frac {d n_k^{(out)} }{dt} &= n_S e^{- k E} - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) n_k^{(out)} + \alpha e^E n_{k+1}^{(out)} + \alpha e^{\Delta } n_{k-1}^{(out)}, \quad k \in \mathbb N_0. \end{align}

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

\begin{equation*} \frac {d}{dt} M^{(out)}(t)=\frac {d}{dt} n_S(t) + \frac {d}{dt} \sum _{k=0}^\infty n_k^{(out)}(t) = 0. \end{equation*}

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

(4.3) \begin{align} 0 &= \tilde {n}_S - e^{\sigma } \tilde {n}_0^{(out)} + \alpha e^E \tilde {n}_1^{(out)} - \alpha e^{\Delta } \tilde {n}_0^{(out)}\\[-10pt]\nonumber \end{align}

(4.4) \begin{align} 0 &= \tilde {n}_S e^{- k E} - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) \tilde {n}_k^{(out)} + \alpha e^E n_{k+1}^{(out)} + \alpha e^{\Delta } \tilde {n}_{k-1}^{(out)}, \quad k \in \mathbb N_0. \end{align}

In order to study the solutions to this equation, we consider the case of $\Delta = \Delta _c$ and $\Delta \neq \Delta _c$ separately.

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 \gt \Delta _c$ or $\Delta \lt \Delta _c$ . For completeness, we also study the case $\Delta = \Delta _c$ .

Explicit solution to (4.3) in the subcritical and supercritical case $\boldsymbol{\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.

(4.5) \begin{align} 0 &= - e^{\sigma } \tilde {n}_0^{(out)} + \alpha e^E \tilde {n}_1^{(out)} - \alpha e^{\Delta } \tilde {n}_0^{(out)} \nonumber\\[3pt] 0 &= - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) \tilde {n}_k^{(out)} + \alpha e^E \tilde {n}_{k+1}^{(out)} + \alpha e^{\Delta } \tilde {n}_{k-1}^{(out)}, \quad k \in \mathbb N_0. \end{align}

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

\begin{equation*} n_k^{part} = e^{- k E } A(0) \overline n_S M(t) = e^{- k E } A(0) \tilde {n}_S(t) , \quad k \geq 1 \end{equation*}

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

\begin{equation*} A(0)\,:\!=\, \frac {1}{e^\sigma - \alpha (e^\Delta -1) ( e^E-1 )}. \end{equation*}

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

(4.6) \begin{equation} \theta ^2 (0)- \frac {\Omega (0)}{\alpha e^E} \theta (0) + e^{\Delta - E} =0, \end{equation}

where

\begin{equation*} \Omega (0)\,:\!=\, e^\sigma + \alpha (e^E + e^\Delta ). \end{equation*}

As a consequence, we have that the solutions $\theta _{1}(0)$ and $\theta _2(0)$ to the equation (4.6) are

\begin{equation*} \theta _{12}(0)= \frac {1}{2 \alpha e^E } \left ( \Omega (0) \pm \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{\Delta + E} } \right ). \end{equation*}

Now notice that

\begin{equation*} \theta _1 (0)= \frac {1}{2 \alpha e^E } \left ( \Omega (0) - \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{\Delta + E} } \right ) \lt 1 \text{ while } \theta _2 (0)= \frac {1}{2 \alpha e^E } \left ( \Omega (0) + \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{\Delta + E} } \right ) \gt 1. \end{equation*}

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 \lt 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

\begin{equation*} \tilde {n}_k^{(out)}(t) = e^{- k E } A(0) \overline n_S M(t) + F(t) \theta _1^k(0),\ k \in \mathbb N. \end{equation*}

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

\begin{equation*} \tilde {n}_0^{(out)}(t) = A(0) \tilde {n}_S(t) + F(t), \quad t\gt 0. \end{equation*}

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

\begin{align*} \tilde {n}_0^{(out) } (t) &= \frac { A (0) \tilde {n}_S(t)+ \alpha e^E \tilde {n}_1^{(out)}(t) }{e^\sigma + \alpha e^\Delta } = \frac { A(0) \tilde {n}_S(t) + \alpha e^E \tilde {n}_1^{(out)} (t)}{e^\sigma + \alpha e^\Delta } \\ &= \frac { \tilde {n}_S(t) A(0) + \alpha A(0) \tilde {n}_S(t) + \alpha e^E F(t) \theta _1(0) }{e^\sigma + \alpha e^\Delta }. \end{align*}

We deduce that

\begin{equation*} F(t)= \frac {A(0) \tilde {n}_S(t) (1+\alpha - e^\sigma - \alpha e^\Delta ) }{ e^\sigma + \alpha e^\Delta - \alpha \varphi (\sigma )} \end{equation*}

where we have that

(4.7) \begin{equation} \varphi (\sigma )\,:\!=\, \frac {1}{2 \alpha } \left ( \Omega (0) - \sqrt { \Omega (0)^2- 4 \alpha ^2 e^{\Delta + E }} \right ). \end{equation}

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

(4.8) \begin{equation} n_k^{(out)} (t) \sim \tilde {n}_k^{(out)} (t) = e^{- k E } A(0) \tilde {n}_S(t) \left [ 1 + B \varphi (\sigma )^k \right ],\ N- k \gg 1 , \quad N \to \infty \end{equation}

for

(4.9) \begin{equation} B= \frac { 1+\alpha - e^\sigma - \alpha e^\Delta }{ e^\sigma + \alpha e^\Delta - \alpha \varphi (\sigma ) }. \end{equation}

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 \lt \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 \lt \Delta _c$ . A direct computation shows that in this case we have that $\varphi (\sigma ) \lt 1$ . As a consequence, we deduce that in the subcritical case (i.e. when $\Delta \lt \Delta _c$ ), we have that

(4.10) \begin{equation} \tilde {n}_k^{(out)} (t) \sim A(0) \overline n_S M(t) e^{- k E } \text{ as } k \to \infty . \end{equation}

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

4.1.2. Asymptotic behaviour of $n_k^{(out)}$ when $\Delta \gt \Delta _c$ and $k \gg 1$ , $N-k \gg 1$

On the contrary, in the supercritical case in which $\Delta \gt \Delta _c$ , we have that $\varphi (\sigma ) \gt 1$ , hence

(4.11) \begin{equation} \tilde {n}_k^{(out)} (t) \sim A(0) \overline n_S M(t) B e^{ - k \psi (\sigma , E) } \text{ as } k \to \infty , \end{equation}

where we recall that $B$ is defined by (4.9) and $\psi$ is given by (2.11), 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) B e^{ - 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 \gt \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 \lt \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 \gt \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

\begin{equation*} \tilde {n}_k^{part} (t) = \frac { k e^{- k E } \tilde {n}_S(t)}{ \alpha (e^{\Delta + E} -1 )}, \quad k \in \mathbb N \end{equation*}

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

\begin{equation*} 0 =- \alpha \left (1+e^{E+ \Delta } \right ) \tilde {n}_k^{(out)} + \alpha e^E \tilde {n}_{k+1}^{(out)} + \alpha e^{\Delta } \tilde {n}_{k-1}^{(out)}, \quad k \in \mathbb N_0. \end{equation*}

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

\begin{equation*} \theta (0)^2 - (e^\Delta + e^{- E} ) \theta (0) + e^{\Delta - E } =0. \end{equation*}

Notice that the two solutions $\theta _{1}(0)$ and $\theta _2(0)$ to the second-order equation are given by

\begin{equation*} \theta _2(0)=e^\Delta , \quad \theta _1(0)=e^{- E}. \end{equation*}

As before, the only admissible solution is $\theta _1(0)\lt 1$ . Therefore, the solution to equation (4.4) is

(4.12) \begin{equation} \tilde {n}_k^{(out)}(t) =e^{-k E } \left ( \frac {k \tilde {n}_S(t) }{ \alpha (e^{\Delta + E} -1 )}+ F(t) \right ), \quad k \in \mathbb N \end{equation}

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

\begin{equation*} 0=\tilde {n}_S (t) -(e^\sigma +\alpha e^\Delta ) \tilde {n}_0^{(out)} (t) + \alpha e^E \tilde {n}_1^{(out)} (t) = \tilde {n}_S (t) -(e^\sigma +\alpha e^\Delta ) \tilde {n}_0^{(out)} (t)+ \alpha \left ( \frac {\tilde {n}_S (t)}{ \alpha (e^{\Delta + E} -1 )}+ F(t) \right ) \end{equation*}

which implies that

\begin{equation*} \tilde {n}_0^{(out)}(t)=\frac {e^\Delta }{ \alpha (e^{\Delta +E} -1 )(e^\Delta -1 )} \tilde {n}_S(t) . \end{equation*}

Asymptotic behaviour of $\boldsymbol{n_k^{(out)}}$ as $\boldsymbol{k \to \infty}$ when $\boldsymbol{\Delta = \Delta _c}$ . Equality (4.12) implies that

(4.13) \begin{equation} \tilde {n}_k^{(out)} (t) \sim \frac { k e^{-k E } \tilde {n}_S(t) }{ \alpha (e^{\Delta + E} -1 )}\text{ as } k \to \infty . \end{equation}

The quasi steady state approximation implies that $ n_k^{(out)} (t) \approx \frac { k e^{-k E } \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 \lt \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.4), i.e.

(4.14) \begin{align} \frac {d n_k }{dt} &= n_S e^{- k E} - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) n_k + \alpha e^E n_{k+1} + \alpha e^{\Delta } n_{k-1} - \mu n_k, \quad k \in \{0, \ldots , N \} \nonumber\\ n_{ N+1}&=0. \end{align}

Notice that we are rewriting the equation for $n_N$ in (2.4), artificially introducing the state $N+1$ in the model and imposing that $n_{N+1}=0$ . This allows us to have the same equation to describe the evolution of $n_k$ for $k \lt 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^{- b N }$ as $N \to \infty$ , this quasi steady state approximation for the solution to equation (4.14) satisfies the following system of equations

(4.15) \begin{align} 0 &= \tilde {n}_S e^{- k E} - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) \tilde {n}_k + \alpha e^{E } \tilde {n}_{k+1} + \alpha e^{\Delta } \tilde {n}_{k-1}, \quad k \in \{0, \ldots , N \} \nonumber\\ \tilde {n}_{N+1}&=0. \end{align}

Notice that the condition $ \tilde {n}_{N+1} =0$ plays the role of an absorbing boundary condition; indeed, when a complex reaches the state $N+1$ , it is lost. This boundary condition, together with the n 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 the subcritical case from the supercritical case separately.

4.2.1. Subcritical case $\Delta \lt \Delta _c$

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

(4.16) \begin{equation} h_m (t) \,:\!=\, e^{(N-m)E} \tilde {n}_{N-m }(t) =e^{k E} \tilde {n}_k (t) \end{equation}

where we recall that $N-k=m$ . Therefore, we now study the following equation for $h_m$

(4.17) \begin{align} 0 &= \tilde {n}_S(t) - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) h_m (t)+ \alpha h_{m-1} (t) + \alpha e^{\Delta + E} h_{m+1} (t), \quad m \in \mathbb N \nonumber\\ h_{-1}(t)&=0, \nonumber\\ \lim _{m \to \infty } h_m(t) &= A(0) \tilde {n}_S(t) . \end{align}

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

\begin{equation*} h^{(h)}_m =\Theta ^m, \quad \forall m \in \mathbb N. \end{equation*}

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

\begin{equation*} \Theta ^2 - \frac { \Omega (0)}{\alpha e^{E+\Delta }} \Theta + e^{-(E+ \Delta ) } =0 . \end{equation*}

This equation has two solutions $\Theta _{1}$ and $\Theta _2$ given by

\begin{equation*} \Theta _{1} = \frac {1}{2 \alpha e^{E+\Delta } } \left ( \Omega (0) - \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) \text{ and } \Theta _{2} = \frac {1}{2 \alpha e^{E+\Delta } } \left ( \Omega (0) + \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right )\! . \end{equation*}

Notice that since we expect the limit as $m \to \infty$ of $ h_m$ to be finite for every fixed time $t\gt 0$ , then $\Theta ^m$ is an admissible homogeneous solution only when $\Theta \leq 1$ . Moreover, since $\Delta \lt \Delta _c$ , then we have that

\begin{equation*} \Theta _2 = \frac {1}{2 \alpha e^{E+\Delta } } \left ( \Omega (0) + \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) \gt 1. \end{equation*}

Hence, there exists a unique admissible homogeneous solution, which is

\begin{equation*} \Theta _1= \frac {1}{2 \alpha e^{E+\Delta } } \left ( \Omega (0) - \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) = \varphi (\sigma ) e^{- (E+\Delta ) } \lt 1. \end{equation*}

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

\begin{align*} h_m(t) = h_m^{(p)} (t)+ h_m^{(h)} (t)= A(0) \tilde {n}_S (t)+ F(t) \Theta ^m = A(0) \tilde {n}_S(t) + F(t) \varphi (\sigma )^m e^{- m (E+ \Delta ) }, \quad m \in \mathbb N \cup \{ -1\}. \end{align*}

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

\begin{equation*} F(t) = - A(0) \varphi (\sigma ) e^{- (E+\Delta ) } \tilde {n}_S(t) . \end{equation*}

Summarizing, we have that

\begin{equation*} h_m = A(0) \tilde {n}_S ( 1 - \varphi (\sigma )^{m+1} e^{- (m+1) (E+ \Delta ) })\quad m \in \mathbb N \cup \{ -1\}. \end{equation*}

Hence, we obtain that

(4.18) \begin{equation} \tilde {n}_{N} (t) \sim e^{- N E } h_0= e^{- N E } A(0) ( 1 - \varphi (\sigma ) e^{- (E+ \Delta ) }) \tilde {n}_S(t) \text{ as } N \to \infty . \end{equation}

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

\begin{equation*} n_{N} (t) \sim e^{- N E } h_0= e^{- N E } A(0) ( 1 - \varphi (\sigma ) e^{- (E+ \Delta ) }) \tilde {n}_S(t) \text{ as } N \to \infty , \ t \approx N. \end{equation*}

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 \lt \Delta _c$ , does not depend on $\sigma$ and where $F(t)$ is a suitable function of time.

4.2.2. Supercritical case $\Delta \gt \Delta _c$

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

\begin{equation*} \ell _k = e^{ k \lambda (\sigma ) } \tilde {n}_k = e^{ k E - k \log (\varphi (\sigma )) } \tilde {n}_k, \end{equation*}

where $\lambda$ is defined as (2.10). Performing the change of variable in equation (4.15), we deduce that $\ell _k$ must satisfy the following equation

\begin{align*} 0 &= \tilde {n}_S(t) - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta) e^{ k \log (\varphi (\sigma )) } \ell _k(t) + \alpha e^{ (k +1)\log (\varphi (\sigma )) } \ell _{k+1} (t) + \alpha e^{\Delta + E + ( k-1) \log (\varphi (\sigma )) } \ell _{k-1}(t) \\ 0& =\ell _{N+1}(t) \end{align*}

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

\begin{align*} 0 & \sim - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) \varphi (\sigma ) \ell _k (t)+ \alpha \varphi (\sigma )^2 \ell _{k+1}(t) + \alpha e^{\Delta + E } \ell _{k-1} (t), \text{ as } k \to \infty \\ 0&= \ell _ {N+1}(t) \end{align*}

Moreover, we define $h_m = \ell _{N-m }$ . Then $h_m$ satisfies the following equation

\begin{align*} 0 & = - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) \varphi (\sigma ) h_m(t) + \alpha \varphi (\sigma )^2 h_{m-1} (t) + \alpha e^{\Delta + E } h_{m+1} (t), \quad m \in \mathbb N \\ h_{-1}(t)&=0 \\ \lim _{m \to \infty } h_m (t)&=B A (0) \tilde {n}_S(t) \end{align*}

where the last equality comes from the matching condition; indeed, 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

\begin{equation*} \alpha e^{\Delta + E } \Theta ^2 - \Omega (0) \varphi (\sigma ) \Theta + \alpha (\varphi (\sigma ))^2 =0. \end{equation*}

We deduce that

\begin{equation*} \Theta _{12} = \frac {\varphi (\sigma )}{2\alpha e^{\Delta + E }} \left ( \Omega (0) \pm \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) \end{equation*}

In particular, we obtain that

\begin{equation*} \Theta _1 = \frac {\varphi (\sigma )}{2\alpha e^{\Delta + E }} \left ( \Omega (0) - \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) = \varphi (\sigma )^2 e^{- (\Delta + E ) } \end{equation*}

and

\begin{equation*} \Theta _2 = \frac {\varphi (\sigma )}{2\alpha e^{\Delta + E }} \left ( \Omega (0) + \sqrt {\Omega (0)^2 - 4 \alpha ^2 e^{E+\Delta }} \right ) = 1 \end{equation*}

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

\begin{equation*} h_m(t) = F_1(t) + F_2(t) e^{ 2m \log \varphi (\sigma ) - m (\Delta + E ) }, \quad \forall m \in \mathbb N \cup \{ -1 \}. \end{equation*}

Matching the behaviour in the inner with the behaviour in the outer region implies that $F_1 (t)= B A(0) n_S(t)$ . On the other hand, since $h_{-1} (t)=0$ , we also have that

\begin{equation*} F_2(t) = - F_1 (t) e^{ -(\Delta + E)} (\varphi (\sigma ))^2= - B A(0) n_S(t)e^{- ( \Delta + E )}(\varphi (\sigma ))^2, \end{equation*}

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

\begin{equation*} h_m(t) = B A(0) \left ( 1 - e^{ 2(m+1) \log \varphi (\sigma ) - (m+1) (\Delta + E ) } \right ) \tilde {n}_S (t), \quad \forall m \in \mathbb N \cup \{ -1 \} \end{equation*}

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

(4.19) \begin{equation} n_{N}(t) \sim e^{- NE + N \log (\varphi (\sigma )) } h_0 (t)= e^{- N E + N \log (\varphi (\sigma )) } A(0) B ( 1 - ( \varphi (\sigma ))^2 e^{ - (\Delta + E ) } ) \tilde {n}_S (t), \text{ as } N \to \infty . \end{equation}

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^{k E}}{k} \tilde {n}_k$ in equation (4.15) and deduce that $\ell _k$ satisfies

(4.20) \begin{align} 0 &= \tilde {n}_S (t) - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) k \ell _k (t) + \alpha (k+1) \ell _{k+1} (t) + \alpha e^{\Delta +E} (k-1) \ell _{k-1} (t), \quad k \in \{ 0, \ldots , N \} \nonumber\\ 0 &=\ell _{N+1} (t). \end{align}

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

(4.21) \begin{align} 0 &= - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) k \ell _k (t) + \alpha (k+1) \ell _{k+1} (t) + \alpha e^{\Delta +E} (k-1) \ell _{k-1} (t) \text{ as } k \to \infty \nonumber\\ \ell _{N+1}(t) &=0. \end{align}

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

(4.22) \begin{align} 0 &= - ( e^{\sigma } + \alpha e^E + \alpha e^\Delta ) h_m(t) + \alpha h_{m-1} (t) + \alpha e^{\Delta +E} h_{m+1} (t)\nonumber \\ h_{-1}(t)&=0, \quad m \in \mathbb N \nonumber \\ \lim _{m \to \infty } h_m (t) &= \frac { \tilde {n}_S(t) }{\alpha (e^{\Delta + E } -1 ) } . \end{align}

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

\begin{equation*} \theta ^2 - \left ( \frac { 1+ e^{\Delta + E } }{e^{\Delta + E }}\right ) + \frac {1}{e^{\Delta + E }}=0 . \end{equation*}

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

\begin{equation*} h_m(t) = F_1(t) + F_2(t) e^{-m (\Delta + E) }, \end{equation*}

with $F_2=-F_1e^{- (\Delta + E) }$ and with

\begin{equation*} F_1(t) = \frac {\tilde {n}_S (t)}{\alpha (e^{\Delta + E} - 1 ) }, \end{equation*}

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

(4.23) \begin{equation} n_{N} (t) \sim \frac {\tilde {n}_S(t)}{\alpha e^{E+\Delta } } N e^{-N E } \text{ as } N \to \infty \text{ and } t \approx N. \end{equation}

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.7) 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 $\boldsymbol{\Delta \lt \Delta _c}$ . Using the asymptotics for $n_{N}$ in (4.18), we approximate the probability of response as follows

\begin{align*} p_{res} (\sigma ) = \alpha e^\Delta \int _0^\infty n_{N} (t) dt \approx \alpha e^{\Delta } e^{-NE} A(0) \overline n_S \left (1- \varphi (\sigma ) e^{-(E+ \Delta )} \right ) \int _0^\infty M (t) dt \text{ as } N \to \infty . \end{align*}

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

(4.24) \begin{equation} \frac {dM }{dt} = - \mu M - \alpha e^{\Delta } n_{N} \sim - ( \mu + \alpha e^{\Delta } e^{-NE} A(0) \overline n_S (1- \varphi (\sigma ) e^{-(E+ \Delta )} ) ) M (t) \text{ as } N \to \infty . \end{equation}

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

(4.25) \begin{equation} M(t) \sim e^{- \left ( \mu + \alpha e^{\Delta } e^{-NE} A(0) \overline n_S \left (1- \varphi (\sigma ) e^{-(E+ \Delta )} \right ) \right ) t } \text{ as } N \to \infty . \end{equation}

Using this approximation for $M$ , we deduce that

(4.26) \begin{equation} p_{res} (\sigma ) \approx \left ( 1 +\frac {\mu e^{NE}}{ \alpha e^{\Delta } A(0) \overline n_S \left (1- \varphi (\sigma ) e^{-(E+ \Delta )} \right ) } \right )^{-1} \text{ as } N \to \infty . \end{equation}

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

\begin{equation*} p_{res}(\sigma ) \approx ( 1 + C_1e^{N(E- b)} )^{-1} \text{ as } N \to \infty . \end{equation*}

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\gt b$ , $p_N \approx 1$ as $N \to \infty$ if $E\lt b$ . Notice that in this case, if we plot the probability $p_{res} (\sigma )$ , defined as (2.7), as a function of $\sigma$ , we do not have a sharp transition from $0$ 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 $\boldsymbol{\Delta = \Delta _c}$ . Arguing as in the subcritical case $\Delta \lt \Delta _c$ , we deduce that (4.23) implies

\begin{equation*} p_{res} (\sigma ) \approx \left ( 1 + \frac {e^{N(E-b)}}{ N \overline n_S } \right )^{-1} \text{ as } N \to \infty \end{equation*}

We conclude that $p_{res} (\sigma ) \approx 0$ as $N \to \infty$ when $E\gt 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 \lt \Delta _c$ .

Supercritical case $\boldsymbol{\Delta \gt \Delta _c}$ . Using the approximation for $n_{N}$ as $N \to \infty$ given in (4.19), we deduce that

(4.27) \begin{equation} p_{res} (\sigma ) \approx \left ( 1+C_2 e^{(E -b- \log (\varphi (\sigma ) ) N }\right )^{-1} \end{equation}

where $C_2= \left [ \alpha e^\Delta B A(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

\begin{equation*} \varphi (\sigma )= \frac {1}{2 \alpha } \left ( e^\sigma + \alpha (e^\Delta + e^E) - \sqrt {(e^\sigma + \alpha (e^\Delta + e^E) )^2 - 4 \alpha ^2 e^{\Delta + E} } \right ). \end{equation*}

Moreover, notice that

\begin{equation*} \varphi '(\sigma )= \frac {e^\sigma }{2 \alpha } \left ( 1 - \frac {e^\sigma + \alpha (e^\Delta + e^E) }{ \sqrt {(e^\sigma + \alpha (e^\Delta + e^E) )^2 - 4 \alpha ^2 e^{\Delta + E} } } \right ) \lt 0, \end{equation*}

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

(4.28) \begin{equation} \lim _{\sigma \to \infty } \lambda (\sigma , E) = E\gt b \gt \max \{ E - \Delta , 0 \} = \lim _{\sigma \to - \infty } \lambda (\sigma , E), \end{equation}

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 $0$ , 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

\begin{equation*} |\lambda (\sigma , E) - \lambda (\sigma _c, E) | \lesssim \frac {1}{N }. \end{equation*}

This region tends to $0$ as $N \to \infty$ . If instead we have that $\sigma$ is such that $\lambda (\sigma , E ) - \lambda (\sigma _c, E ) \gt \frac {1}{N }$ then we have $p_{res}(\sigma ) \approx 1$ as $N \to \infty$ and if $\lambda (\sigma _c, E) - \lambda (\sigma , E) \gt \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 \gt \Delta _c$ relies on the fact that we assume that the interaction between the ligands and the receptors takes place in a characteristic time, the signalling lifetime $\mu ^{-1 }$ . If we assume that the ligand can interact with the receptor for an arbitrarily long amount of time, the discrimination property is lost.