Algebraic study of receptor-ligand systems: a dose-response analysis

The study of a receptor-ligand system generally relies on the analysis of its dose-response (or concentration-effect) curve, which quantifies the relation between ligand concentration and the biological effect (or cellular response) induced when binding its specific cell surface receptor. Mathematical models of receptor-ligand systems have been developed to compute a dose-response curve under the assumption that the biological effect is proportional to the number of ligand-bound receptors. Given a dose-response curve, two quantities (or metrics) have been defined to characterise the properties of the ligand-receptor system under consideration: amplitude and potency (or half-maximal effective concentration, and denoted by EC$_{50}$). Both the amplitude and the EC$_{50}$ are key quantities commonly used in pharmaco-dynamic modelling, yet a comprehensive mathematical investigation of the behaviour of these two metrics is still outstanding; for a large (and important) family of receptors, called cytokine receptors, we still do not know how amplitude and EC$_{50}$ depend on receptor copy numbers. Here we make use of algebraic approaches (Gr\"obner basis) to study these metrics for a large class of receptor-ligand models, with a focus on cytokine receptors. In particular, we introduce a method, making use of two motivating examples based on the interleukin-7 (IL-7) receptor, to compute analytic expressions for the amplitude and the EC$_{50}$. We then extend the method to a wider class of receptor-ligand systems, sequential receptor-ligand systems with extrinsic kinase, and provide some examples.


Introduction
The human body consists of more than 3 × 10 13 cells [1], each of them receiving, at any given time, hundreds of signals from extra-cellular molecules when these bind their specific membrane receptors.These signals are integrated, translated and read by a small number of intra-cellular molecules to generate appropriate cellular responses [2].Surface receptors specifically bind to extra-cellular molecules called ligands.The binding of a ligand to its receptor induces an intra-cellular cascade of signalling events which regulate a cell's fate, such as migration, proliferation, death, or differentiation [3,4].Receptor-ligand interactions are essential in cell-to-cell communication, as is the case for immune cell populations [5], and thus, a large body of literature has been devoted to the experimental and theoretical study of cell signalling dynamics [6,7,8,9,10,11,12,13,14].Exploiting the controlled environment of in vitro experiments, most cell signalling studies focus on the estimation of the affinity constant for a given receptor-ligand system, and the quantification of biochemical on and off rates for the binding and unbinding, respectively, of receptor and ligand molecules.Recent single-cell studies have shown that cells have heterogeneous expression levels of receptor copy numbers.Not only does the copy number depend on the cell type, but receptor copy numbers vary strongly between isogenic cells of one cell type [14,15,5].Given the heterogeneity of receptor copy numbers across and within cell types, it is timely to understand how a cell's response to a given ligand depends on the expression levels of its receptor.This quantification will be a first step to account for the variability of receptor expression levels when designing and studying receptor-ligand models (both from an experimental and mathematical perspective) [5,15,8,7].
The study of a receptor-ligand system generally relies on the analysis of its dose-response (or concentration-effect) curve, which describes the relation between ligand concentration and the biological effect (or cellular response) it generates when binding its specific receptor [3,16,11].Mathematical models of receptor-ligand systems have been developed to compute a dose-response curve, under the assumption that a biological effect is proportional to the number of ligandbound receptors [13,15,5,6].Given a dose-response curve, two quantities (or metrics) have been defined to characterise the properties of the ligand-receptor system under consideration.These metrics are: the amplitude, which is defined as the difference between the maximal and minimal response, and the half-maximal effective concentration (or EC50), which is the concentration of ligand required to induce an effect corresponding to 50% of the amplitude [3,16,11].The amplitude is a measure of the efficacy of the ligand, and the EC50, a measure of the potency (or sensitivity) of the ligand (for a given receptor) [3,16,17].Both amplitude and EC50 are key quantities commonly used in pharmaco-dynamic modelling, yet a comprehensive mathematical investigation of the behaviour of these two metrics is still outstanding for most receptorligand systems.For instance, for a large (and important) family of receptors, called cytokine receptors [18,19,20], we still do not know how amplitude and EC50 depend on receptor copy numbers (for a given concentration of ligand) [14,15].
In this paper we bridge this gap by deriving closed-form expressions for a class of cytokine-receptor models.We further highlight how tools from computational algebra can be used to facilitate the calculation of both the amplitude and the EC50 for this family of models.
Previous work has shown that the estimation of the amplitude and the EC50 from experimental data is often possible, although strong inductive biases might be introduced [16,3].Usually one starts with a data set where the number (or concentration) of receptor-ligand signalling complexes formed (see Section 2.2) is measured for different values of the ligand concentration.Then, the estimation of the amplitude and the EC50 is turned into a regression problem by assuming a functional relationship in the data set and fitting a parametric curve.A simple first approach is to plot experimental values (corresponding to a measurable variable which quantifies cellular response) as a function of ligand concentration.The amplitude and the EC50 are then read directly from a curve formed by interpolation of the data points.Since the EC50 is likely to fall between two data points, a geometrical method [21] can be used for an accurate determination.Nowadays many software packages can compute the amplitude and the EC50 from the data set making use of statistical methods, which consist in finding the "best-fit" equation to the dose-response curve.The most common shape of the dose-response curve is a sigmoid, and thus, can be fitted with the famous Hill equation [22,23].However, other functions are also possible, such as a logistic equation [24,25], a log-logistic equation [26,27], or the Emax model [28,29].An asymmetrical sigmoid equation is sometimes needed for better precision [24,27].The amplitude and the EC50 are parameters of these equations and can thus, be directly inferred from the fitting process.When a data set does not follow the strictly increasing pattern of these Hill-like functions, then more complex functions, such as bell-shaped curves [30], or multi-phasic curves [31] can be used.It is important to note that even though these empirical regression methods allow one to quantify the two key receptor-ligand metrics, amplitude and EC50, they do not offer any mechanistic insights for the receptor-ligand system under consideration.To this end, mathematical models can be used to describe the receptorligand system at a molecular level; that is, mathematical models consider the biochemical reactions which initiate a cellular response [32,6,11].The challenge in such models is finding analytical, ideally closed-form, expressions for the amplitude and the EC50.Due to the non-linear nature of the biochemical reactions involved, this poses a significant and practical challenge.
Cytokine-receptor systems are of great relevance in immunology [18,19,20], and here we want to address this challenge in the context of this family of receptors [33,20].The advantages of having analytical (or closed-form) expressions of the amplitude and the EC50 for a large class of receptor-ligand systems are many: i) they allow to quantify their dependence on receptor copy numbers, ii) they facilitate mathematical model validation and parameter exploration, and iii) they reduce computational cost.To the best of our knowledge such expressions have been obtained in a few instances: closed or open bi-molecular receptor-ligand systems [34], monomeric receptors [35], or ternary complexes [36].More complicated receptor-ligand models have been studied with chemical reaction network theory (CRNT) [37,38,39,38], but CNRT has thus far, been focused on the analysis of the steady state of the system (i.e., existence and number of steady states and their stability).Yet, we believe CRNT is an essential and useful framework to start any mathematical investigation of the amplitude and the EC50.
Another aspect which can be effectively addressed by mechanistic mathematical modelling is the effect of internal or external perturbations to the state of a cell.For example, in single-cell experiments or even repetitions of bulk experiments [15,14], the experimental conditions can never be replicated exactly.This leads to noise not only in the measured quantities, but also in the reaction mechanisms themselves.This variation can be captured in mathematical models which encode parameters such as affinity constants or total copy number of constituent molecular species.An analytical study of the dependency of pharmacologically relevant quantities, such as amplitude and EC50, on the reaction parameters can facilitate in silico drug design [40].While amplitude and EC50 are widely employed to characterise biological phenomena, the manner in which they depend on the parameters of the receptor-ligand model is not fully understood.Thus, improved understanding of these relationships could provide novel biological and computational insights.
Motivated by the previous challenges and making use of methods from CRNT and algebraic geometry, such as the Gröbner basis, in this paper we propose a new method to obtain analytic expressions of the amplitude and the EC50 for a large class of receptor-ligand models, with a focus on cytokine receptors.The paper is organised as follows.In Section 2 we introduce the mathematical background and essential notions of CRNT used in the following sections.With IL-7 cytokine receptor as a paradigm, in Section 3 we propose a general method to calculate the amplitude and the EC50 of the dose-response curve for a class of receptor-ligand systems.In Section 4 we generalise the previous results to a wider class of receptor-ligand systems, sequential receptor-ligand systems with extrinsic kinase, and provide a few biological examples of these systems.Finally, we discuss and summarise our results in Section 5. We have included an appendix to provide additional details of our methods (perturbation theory) and our algebraic computations.

Mathematical background
In this section we briefly summarise the relevant notions of chemical reaction network theory and formally define amplitude, EC50, and signalling function.A very short introduction to the use of Gröbner bases is also given.

A brief introduction to chemical reaction network theory
In this paper we view a chemical reaction network (CRN), N , as a multi-set N = {S, C, R}, where S is the set of species, C the set of complexes, and R the set of reactions.We note that in the context of CRN, a "complex" is a linear combination of species and need not be a "biological functional unit", which we refer to as a biological complex.We denote, whenever useful, a biological complex formed by species X and Y as X : Y , where the colon denotes the physical bond between X and Y .The order of species in the biological complex is irrelevant, i.e., X : Y = Y : X.
Example 2.1 (Heterodimeric receptor tyrosine kinase).A simple heterodimeric receptor tyrosine kinase (RTK) model has a species set S = {X1, X2, Y1}, a complex set C = {X1 + X2, Y1, Y2}, and a reaction set R = {X1 + X2 → Y1, Y1 → X1 + X2, Y1 → Y2, Y2 → Y1}.Ligand binding induces dimerisation of these receptors resulting in auto-phosphorylation of their cytoplasmic domains (tyrosine autophosphorylation sites) [41].X1 and X2 are the two components of the heterodimeric RTK.The biological complexes Y1 = X1 : X2 and Y2 = L : Y1 are the heterodimeric receptor with intrinsic kinase activity and the heterodimeric receptor bound to the ligand, respectively.In this paper the ligand concentration (L) is taken to be an input parameter and, hence, it does not feature as a separate chemical species in the species set S.
We can associate a reaction graph to every CRN N , by letting the vertex set be C and the (directed) edge set R. There exists a class of important CRNs defined by their network reversibility.

Definition 2.2 (Network reversibility).
Let N be a CRN with its associated reaction graph G(C, R).An edge between Ci and Cj ∈ C exists if Ci → Cj ∈ R. If for every edge Ci → Cj ∈ R, the edge Cj → Ci ∈ R also exists, then the network is called reversible.If for every edge, Ci → Cj ∈ R, a directed path exists going back from Cj to Ci, then the network is called weakly reversible.All reversible networks are also weakly reversible.
A general reaction from complex Ci to complex Cj can be written as where the sum is over the set of species (X1, X2, . . ., Xn), and αi = (αi1, ..., αin) T and αj = (αj1, ..., αjn) T are nonnegative integer vectors.The corresponding reaction vector is given by r = αj − αi.For a CRN with n species and m reactions we can now define the n × m matrix of all reaction vectors, Γ, such that Γ = (r1, . . ., rm).This matrix is called the stoichiometric matrix.
To derive dynamical properties from the static description so far provided, we make use of the law of mass action kinetics [42].First, we assign a rate constant k ∈ R>0 to each and every reaction in the network.Second, we denote the concentration of species Xi by xi.With this notation, we then associate a monomial to every complex Ci = k α ik X k , as follows where n is the number of species in the network.We define the reactant complex of a reaction as the complex on the left hand side of reaction (1).The reaction rate of a reaction is the monomial of its reactant complex multiplied by the rate constant.The flux vector, R(x), is the m × 1 column vector of all reaction rates.The ordinary differential equations (ODEs) governing the dynamics of the reaction network are given by where Γ is the stoichiometric matrix (defined above).We note that the reaction rate of the i th reactant complex is the i th row in R(x), and similarly, the stoichiometry of i th reaction is given by the i th column of Γ.
From (3) we can also deduce the conserved quantities of the reaction network.That is, if a vector exists, z ∈ Z n , such that d(z T x)/dt = z T ΓR(x) = 0, the quantity z T x is conserved.Consequently, the left kernel of Γ defines a basis for the space of conserved quantities.In this way, conservations induce linear relations between the variables.Informally we say that a molecular species Xi is conserved if its total number of molecules, Ni, is constant.Ni is determined by the initial conditions.
Example 2.4 (Heterodimeric RTK continued).The dynamical system associated with the heterodimeric RTK model is given by with ki as the reaction constants of the forward chemical reactions ( ) and qi the reaction constants of the backward reactions ( ).A basis for the conservation equations is given by the linear relations X1 +Y1 +Y2 = N1 and X2 +Y1 +Y2 = N2.These imply that the total amount of the species X1 (X2) is conserved by adding the amounts of the bound states of the molecule (Y1 and Y2) to the amount of free molecule X1 (X2).
We can now define the biologically relevant steady states of a CRN.
A useful connection between the static network structure (defined earlier) and the existence (and stability) of unique biologically relevant steady states can be made via deficiency theory [37].
Definition 2.6 (Deficiency).Let N be a CRN with connected components in the reaction graph and η = dim span(r1, . . ., rm) be the dimension of the span of the reaction vectors.The deficiency of N is then given by The notion of network deficiency leads to one of the fundamental theorems of CRNT, the Deficiency Zero Theorem [37], which connects the network structure to the dynamics of a CRN.Theorem 2.7 (Deficiency zero theorem).Let N be a weakly reversible CRN with δ = 0. Then the network has a unique biologically relevant steady state for every set of initial conditions, and this steady state is asymptotically stable.
With certain additional conditions on the reaction rates (see Refs. [43,44]), biologically relevant steady states are detailed balanced.This means that for every reaction of the form (1), the steady states satisfy where K = k/q, the ratio of the rate constants of the forward and backward reactions, is called the affinity constant of the reaction.
Example 2.8 (Heterodimeric RTK continued).The heterodimeric RTK model has 3 complexes, 1 connected component and the dimension of the span of the reaction vectors is 2; hence, δ = 3 − 1 − 2 = 0. Since the network is reversible, we know from Theorem 2.7 that there exists exactly one stable positive steady state for each set of initial conditions.One can show that in fact y1 = (k1/q1)x1x2 and y2 = (k2/q2)y1.

Signalling function: amplitude and half-maximal effective concentration
In this paper we want to closely investigate pharmacological properties of receptor-ligand systems, rather than the steady state structure of the models.In particular, we want to study the dose-response (or concentration-effect) curve of the system, which describes the relation between ligand concentration and the biological effect (or cellular response) it generates when binding its specific cell surface receptor.As mentioned in Section 1 a lot of effort has been devoted to explore the steady state structure of chemical reaction networks.In this paper we make use of algebraic methods to explore the dose-response of receptor-ligand systems.To do so we start with the definition of signalling complex.We note that in most biological instances the signalling complex is formed by all the sub-unit chains that make up the full receptor, intra-cellular kinases and the specific ligand [2,4,8,10,14,15,17].
Definition 2.9.The signalling complex of a receptor-ligand system is the biological complex which induces a biological response.
This leads to the following definitions of the signalling function and dose-response curve.Definition 2.10.We define the signalling function, σ : R+ → R+, L → σ(L), as the univariate function which assigns to a given value of ligand concentration, L, the number (or concentration) of signalling complexes formed at steady state.The dose-response curve is the corresponding plot of the signalling function.
We note that in what follows we will not distinguish between number (or concentration) of signalling complexes since one can be obtained from the other if we know the volume of the system and Avogadro's number.
The specific choice of σ will depend on the receptor-ligand system under consideration.In this paper we focus on the class of cytokine receptors and the signalling function will be defined in Section 3. In our examples the signalling function will be a product of the steady state values (numbers) of sub-unit chains that make up the full receptor, intra-cellular kinases, affinity constants of the reactions involved, and ligand concentration.This together with equations ( 2) and ( 3) indicate that the signalling function will always be algebraic.Next, we define a central object of study in this paper; namely, the amplitude of the signalling function, often referred as efficacy in the pharmacology literature [3].Definition 2.11.The amplitude of the signalling function, A, is the difference between the maximum and the minimum of σ; that is, A ≡ max(σ) − min(σ).
We note that when min(σ) = 0, which is the case considered in this paper (min(σ) = σ(0) = 0), the amplitude is given by the maximum of the signalling function.If, in addition, the dose-response curve attains its maximum at large concentrations (for instance, when the dose-response curve is sigmoid), we have The amplitude provides information about the magnitude of the intra-cellular response to the stimulus, L. The larger the amplitude is, the larger the response variability will be.The amplitude is always bounded by the number of molecules available.However, this bound is often not tight [45].To quantify the sensitivity of the model to the stimulus, i.e., the potency of the ligand L, we introduce the half-maximal effective concentration, EC50.

2
= min(σ) + A 2 .We say that the EC50 is inversely proportional to ligand potency; namely, the lower the EC50, the higher the potency of the ligand.Figure 1 illustrates the amplitude and the EC50 of a sigmoid dose-response curve (A) when its minimum is zero: increasing the amplitude shifts up the maximum of the curve and results in greater efficacy (B), and decreasing the EC50 shifts the dose-response curve to the left and increases the potency of the ligand (C).We now review some algebraic and analytic tools which will enable us to compute the EC50 and the amplitude.

Gröbner bases
Since we assume the law of mass action, the models studied in this paper are systems of polynomial equations, and thus, we can use the techniques developed in the field of computational algebra and algebraic geometry [46].Such methods have also been successfully applied to many topics in chemical reaction network theory, see e.g., Refs.[47,48,49].In particular, we make use of Gröbner bases.Informally speaking, a Gröbner basis is a non-linear generalisation of the concept of a basis in linear algebra and, therefore when a Gröbner basis for a polynomial system is calculated, many properties of the system can be investigated, such as the number of solutions and the dimensionality of the space of solutions.Strictly speaking, however, a Gröbner basis is not a basis as it is not unique and it depends on the lexicographical (lex) monomial ordering chosen.For more details we refer the reader to Ref. [46].
A lex Gröbner basis is a triangular polynomial system; that is, for a polynomial system (ideal) in Q[x1, . . ., xn] we obtain a polynomial system of the form We note that when the solution space is positive dimensional, then g1, . . ., gn are identically zero.For a given Gröbner basis with zero-dimensional solution space we can now iteratively, and often numerically, solve the constituent polynomials to obtain the solutions (in C n ) for the polynomial system.We can also find all real and, further, positive solutions, if there are any [46].
1.The system is in steady state.
2. The ligand is in excess (we consider ligand concentration, L, as a parameter instead of a dynamic variable).

A unique biologically relevant solution exists for any given set of rate constants and initial conditions.
The IL-7R models we have chosen are simple enough to illustrate our method, and thus, to derive analytic expressions for the amplitude and the EC50, yet complex enough to show its limitations.

Two motivating examples: IL-7 cytokine receptor as a paradigm
We now consider the cytokine interleukin-7 (IL-7) and its receptor (IL-7R) [50,9,8,10,15,12] as a motivating receptorligand system.IL-7 is a cytokine involved in T cell development, survival and homeostasis [51].Its receptor, IL-7R, is displayed on the surface of T cells and is composed of two trans-membrane chains: the common gamma chain (denoted by γ) and the specific high affinity chain IL-7Rα (denoted by α) [8,12,51,13].Cytokine receptors do not contain intrinsic kinase domains, and thus, make use of Janus family tyrosine kinases (JAKs) and signal in part by the activation of signal transducer and activator of transcription (STAT) proteins [52].In the case of the gamma chain, it binds to the intra-cellular extrinsic Janus kinase molecule, JAK3.Binding of IL-7 to the dimeric JAK3-bound IL-7 receptor, defined as α : γ : JAK3, initiates a series of biochemical reactions from the membrane of the cell to its nucleus, which in turn lead to a cellular response.For the IL-7R system the STAT protein preferentially activated is STAT5 [52], so that the amount of phosphorylated STAT5 can be used as the experimental measure of the intra-cellular response generated by the IL-7 stimulus.The IL-7R receptor is illustrated in Figure 2a, where the hatched area determines the intra-cellular environment.
The first model we consider is shown in Figure 2c.As discussed in Ref. [9], the gamma chain is shared by other cytokine receptors.This model does not include the competition for the gamma chain between different cytokine receptors, therefore later in this section we introduce a second model to account for this competition.In this section we will provide an (algebraic) analytic treatment of both models.We consider the formation of "dummy" receptors, α : γ, which are formed of the IL-7R devoid of JAK3 and, therefore, they cannot signal (see Figure 2b).We further assume no allostery; that is, the affinity constants of the biochemical reactions involved in the formation of the dummy complex, L : α : γ, are the same as the affinity constants involved in the formation of the signalling complex, L : α : γ : JAK3.

The IL-7 receptor-ligand system: two receptor chains and a kinase
We first consider a model in which the IL-7R is formed sequentially, one molecule at a time; the γ chain binds to the kinase, JAK3, then the α chain binds to the complex formed by γ and JAK3.Finally, the ligand, IL-7, binds to the signalling receptor composed of γ, α and JAK3.The model also includes the formation of "dummy" receptors, which do not involve the kinase JAK3. Figure 2c illustrates the sequential formation of the signalling and dummy complexes.The reaction scheme for this model is as follows where for i = 1, 2, 3, Ki is the affinity constant of the appropriate reaction.One can show that this system has deficiency zero and is reversible (see Section 2.1).Therefore, for every set of rate constants and initial conditions, there exists exactly one positive steady state.Moreover, this positive steady state is in detailed balance.We remind the reader that in this paper we assume mass action kinetics to determine reaction rates.We denote the concentration of γ, α, JAK3 and IL-7 by x, y, z, and L, respectively.The reaction rate for the forward/backward reaction ( / ) is given by ki and qi, respectively, for i = 1, 2, 3. We note that Ki = ki/qi.The concentrations of the product complexes of the forward reactions are denoted by ci in order of appearance (see Figure 2c).We can now write down the ordinary differential equations (ODEs) associated to the system of reactions (6):  A suitable basis for the conservation equations is that is, single chain molecules are conserved since we do not consider the generation or degradation of molecules.The constants Nx, Ny and Nz represent the total copy number of γ, α and JAK3 molecules per cell, respectively.Detailed balance leads to the following steady state equations: Substituting the steady state equations into the conservation equations, we obtain a system of polynomials Analytic computation of the steady state.The polynomial system (10) can be solved numerically for a particular set of parameter values.However, an analytic solution will provide greater insight and will allow us to derive expressions for the amplitude and the EC50.We make use of Macaulay2 [53] to compute a lex Gröbner basis for this model, which will lead to a triangular set of polynomials 1 , as follows: Equation (11c) gives where the last equality follows from equation (11a).Solving the system (11) and selecting the biologically relevant solution, we obtain an analytic expression for the number of free (unbound) JAK3, α and γ molecules at steady state where we have introduced and We study the dose-response curve of this model given by the number of signalling complexes, L : γ : α : JAK, per cell at steady state and as a function of L. The signalling function, σ(L), is given by .
Analytic computation of the amplitude.A simple inspection of the behaviour of (12) shows that the doseresponse curve is a sigmoid, such that σ(0) = 0. Therefore the amplitude A is given by the asymptotic behaviour of the signalling function as follows: We will prove this result rigorously for a more general class of models in Section 4. We first notice that z is independent of L. We now compute the product xy (at steady state) as follows xy = (Nx − Ny)y + y 2 1 + K1z .
From equation (11b) we can replace the polynomial in y of degree two by an expression linear in y: .
Thus, we obtain the following analytic expression for the signalling function: Since K 3 L 1+K 3 L → 1 when L → +∞, we need to study the expression Ny − y in this limit.We have where Finally, noticing that Nx + Ny − |Nx − Ny| 2 = min(Nx, Ny), we obtain the amplitude where z is the analytic expression obtained in (12).This result indicates that the amplitude of this model is the total number of the limiting trans-membrane chain modulated by a factor, valued in the interval [0, 1], which only depends on K1, Nx and Nz.
Analytic computation of the EC 50 .We now determine the EC50 by finding the value of L50 such that where x50, y50 and z50 are the steady state expressions found in (12) evaluated at L = L50.Two expressions satisfy this equation but only one provides a relevant biological solution with L, x, y, z > 0. The relevant analytic expression of the EC50 is given by with M = min(Nx, Ny).The details of the computation can be found in Appendix B. This result shows that the EC50 value for this system is independent of the kinase, since the parameters K1 and Nz are absent in the previous expression.
Alternatively, we now propose a more algebraic method to derive the analytic expression of the EC50.We compute a lex Gröbner basis for the augmented system of polynomials consisting of the steady state equations (10) and where this time x, y, z, and L are variables.The resulting triangular system describes directly the EC50 and x, y, z at L =EC50: Solving (23a) and selecting the solution for which y and x, given by equations (23c) and (23d), respectively, are positive yields the final result, in agreement with (21).

3.1.2
The IL-7 receptor-ligand system: an additional sub-unit receptor chain  Decoy complexes can be formed to prevent the formation of signalling or dummy complexes.The mathematical notation used is annotated below each molecule or complex.
The previous model described the IL-7 receptor system without any consideration for the fact that the γ chain is shared with other cytokine receptors [9].We now account for this competition by including in the previous model an additional receptor chain, R, which can bind to the γ chain, or the complex γ : JAK3, to form decoy receptor complexes (see Figure 3a and Figure 3b, where the hatched area indicates the cytoplasmic region).The resulting reaction scheme (summarised in Figure 3c) is given by We use w to describe the concentration of the additional chain R. Similarly to the previous model, we write the system of ODEs describing the time evolution for each complex and then derive (a basis for) the conservation and steady state equations.Combining them, we obtain the following polynomial system: where Nw is the additional conserved quantity.Again, we compute a lex Gröbner basis for this set of polynomials to obtain the following triangular system: where Solving (25a) gives the number of free JAK3 molecules per cell at steady state, z; solving (25b) gives the number of free (unbound) α chains per cell, y; and substituting y and z into (25c) and (25d) gives the remaining steady states.We obtain the following implicit steady state expressions for the number of free (unbound) chains: The problem now reduces to finding the positive real roots of (25b).As (25b) is a polynomial of degree three, we could, in principle, find an exact analytic solution.However, such a solution might not be very informative.Instead, we show how perturbation theory can be used to obtain the amplitude of the dose-response.In this model, the signalling complex is still L : α : γ : JAK3.The signalling function is given by In Section 4.3 we will show that, for this model, the maximum of σ is attained in the limit L → +∞.Hence, we have Combining (25a), written as Nx − Nz + z = Nx 1+K 1 z , and (25b), we obtain a reduced expression for the product xy which allows us to rewrite the amplitude as We note that z is independent of L and, therefore, to compute the amplitude we only need to find the behaviour of y as L → +∞.
Perturbation theory to determine y as L → +∞.We now apply the method described in Ref. [54] and summarised in Appendix A. Let = 1 L and define the polynomial P as follows: where P2 is the polynomial (25b).We added a factor of 2 to remove any negative powers of in P2.We obtain the polynomial where We now replace y by p ω( ) with ω(0) = 0 according to theorem A.3.We obtain where The smallest exponents in the previous equation are We note that 0 is not in E because we multiplied P2 by 2 .Applying the graphical algorithm detailed in Appendix A, we find the proper values (0, 0) and (1, 2) (see Figure 4).We investigate these two branches.Branch (0,0).We make use of the notation in Appendix A, to define T (1) (ω) ≡ 0 P (ω 0 ).
The least common denominator of {2,1,0,0} is q1 = 1.Therefore in accordance with the notation of Appendix A = β, and the polynomial R is the polynomial P .It means that we have y = ω and we can directly carry out a regular perturbation expansion.
Let us write the asymptotic expansion y = y0 + y1 + y2 2 + . . .and replace it in P (y).Since P (y) = 0, by the fundamental theorem of perturbation theory (Theorem A.2) we obtain a system of equations in y0, y1, . .., which can be solved.The first equation of the system is given by We are are only interested in non-negative values of y0, since we want y to be biologically relevant.We also require ω(0) = y(0) = y0 = 0 from Theorem A.3.Thus, solving equation (33), we obtain y0 = Ny − Nx if Ny > Nx and y0 = 0 otherwise.Assuming y0 = 0 (i.e., Nx ≥ Ny), we solve the next order equation We select the positive root of this polynomial and obtain an expression for y1 when y0 = 0. Thus, we have Equation (35) shows that y1 > 0 when Nx ≥ Ny.Hence, y ≈ y1 converges to zero from above and, therefore, represents a biologically relevant solution.We can conclude that Branch (1,2).On this branch, and again following the notation of Appendix 47, we define The least common denominator of {2,1+1,0+1,0+1} is q2 = 1, so R β is the same polynomial as T (2) .Since y = ω , we have Ny − y ∼ →0 Ny.Furthermore, when replacing ω by an asymptotic expansion ω0 + ω1 + . . . in T (2) and applying the fundamental theorem of perturbation theory (Theorem A.2), we obtain the same equation for w0 as for y1 in the previous branch (see equation ( 34)): We have y1 = ω0.In other words, at large, but finite L = 1/ , the convergence behaviour of the two branches is identical.This agrees with Theorem 2.7 which states that there is only one positive solution for each set of reaction constants and initial conditions.In conclusion, we find that Ny − y = min(Nx, Ny), which gives the following expression for the amplitude with z defined in (26).As the steady state concentration of JAK3, z, is the same in the IL-7R model with or without the extra chain R, the amplitude of both models has the exact same expression.
Computation of the EC 50 .Since we did not compute analytic expressions for each steady state concentration, the EC50 expression has to be obtained by computing a Gröbner basis of the polynomial system (24) augmented by the polynomial considering x, y, z and L as variables, with M = min(Nx, Ny).The lex Gröbner basis obtained for this system is: where we wrote: The polynomial (40a) is expected to be independent of the ligand concentration, L. The EC50 expression is the real positive root of polynomial (40b) at which x, y and w (obtained via polynomials (40e), (40c) and (40d), respectively) are positive.The polynomial (40b) reflects the parameter dependency of the EC50: since the parameters K1 and Nz are not present in its coefficients, we can affirm that the EC50 is, once again, independent of the kinase.Thus, we reduced the problem of computing the EC50 to solving a univariate polynomial (equation ( 40b)).In comparison, before any algebraic manipulation was possible, the polynomial system (24) had to be solved multiple times to obtain the dose-response curve (a sigmoid), which was then fitted with a Hill equation.Finally, the EC50 was computed from the fitted parameters of the Hill curve.Alternatively, if one wanted to apply the Gröbner basis-free method of Section 3.1.1,one would have to solve the polynomial (25b) in y (which is possible in theory), find its positive real solution (which is potentially hard), substitute the expression of y into σ(L), and then solve for the EC50.

Summary of proposed algebraic method to study the signalling function
From the two previous examples, we devise a general algorithm to compute analytic expressions of the steady state, the amplitude and the EC50 for some receptor-ligand systems when ligand is in excess.

Key steps
1) Write the mass action kinetics set of ODEs for the system under consideration.
2) Obtain the polynomial system by combining the steady state and conservation equations.
3) Compute a lex Gröbner basis of the polynomial system obtained in step 2. 4) Define the signalling function σ(L).
6) Compute a lex Gröbner basis of the polynomial system obtained in step 2 augmented by the equation with the ligand concentration, L, considered an additional variable.This additional equation corresponds to definition 2.12 of the EC50.7) Find the positive roots of the univariate polynomial in L of the Gröbner basis obtained in step 6.The root which allows the other variables to be positive is the EC50.
One of the crucial parts of the proposed algebraic algorithm is the amplitude computation.Usually, we have the simplification that min(σ) = σ(0) = 0, however, finding max(σ) can be challenging.For certain classes of models we have which greatly reduces the calculation.We can now either solve the Gröbner basis from step 3 directly, to obtain analytic expressions of the steady state concentrations of the single chains components, or use perturbation theory as outlined in Section 3.1.2.In the final step, if an exact expression for the EC50 cannot be computed, i.e., the univariate polynomial in L has a large degree, one already reduces the cost of the EC50 computation compared to the naive approach.In summary, in this section we compute the lex Gröbner bases with the computer algebra package Macaulay2 [53] and provide a Macaulay2 code example in Appendix C.

Analytical study of general sequential receptor-ligand systems
In spite of the general applicability of the method outlined in the previous section, we still have to make the assumption that the computed limit of the signalling function coincides with its amplitude.In this section we show that this is indeed the case for a wider class of receptor-ligand systems.An analytic closed-form expression for the amplitude follows with little extra work.The EC50 can then be studied making use of the extended Gröbner basis introduced in Section 3.2.We start by giving an abstract generalisation of the example from Section 3.1.1.
Note that the IL-7R model of Section 3.1.1 is one example of an SRLK model and the definition of signalling function given in Section 2.2 is equivalent.We now introduce the notion of a limiting component.If the signalling function attains its maximum for large values of the ligand concentration, then, since by definition σ(0) = 0, the amplitude of such model is given by In this section we present some general results for limL→+∞ δ(L) and limL→+∞ σ(L) applicable to SRLK models.The proofs of the lemmas and theorems can be found in Appendix D.

Asymptotic study of the steady states
While it is difficult to find closed-form expressions of the steady states for general receptor-ligand systems, in what follows we show that considerable progress can be made for the specific case of SRLK models.In this section we describe the behaviour of the concentrations, xi, in the limit L → +∞.The proofs of our results can be found in Appendix D. First, we recall the definition and a property of algebraic Definition 4.5.A univariate function y = f (x) is said to be algebraic if it satisfies the polynomial equation: where the Ri(x) are rational functions of x, i.e., are of the form p(x) q(x) , where p and q are polynomial functions and q(x) = 0 for all x ∈ R. Remark 4.6.Note that the polynomial ( †) has m solutions.These solutions are called the branches of an algebraic function and one often specifies a particular branch.
Since we are interested in the limit behaviour, the following lemma proves useful.Lemma 4.7.Any bounded, continuous solution of ( †) defined on R has a finite limit at +∞ (and −∞).
With this background in place, we can now proceed to study SRLK models in detail.We start by showing that in steady state the signalling and the dummy functions have a positive limit when L tends to +∞.
Lemma 4.8.The signalling and the dummy functions of an SRLK model satisfying the experimental hypotheses admit a finite limit when L → +∞ and this limit is positive.
An equivalent result holds for the steady state concentration of the kinase.Lemma 4.9.In an SRLK model under the experimental hypotheses, the concentration of the extrinsic intra-cellular kinase Z admits a positive finite limit, cz > 0, when L → +∞.
In the particular case of no allostery, we can write an explicit expression of the limit of z, cz.Lemma 4.10.Consider an SRLK model which satisfies the experimental hypotheses.If we assume no allostery, then the steady state value of the extrinsic intra-cellular kinase, z, is given by where By Lemma 4.10 z is independent of L (thus, cz = z) and only depends on K1, N1 and Nz.Note that this result is equivalent to the one obtained in Section 3.1.1for the IL-7R model.Finally, we study the behaviour of the concentration xi in the limit L → +∞.We first give bounds to the asymptotic dependency of xi on L. Lemma 4.11.Let us consider an SRLK model which satisfies the experimental hypotheses.Then no concentration xi behaves proportionally to L q , q > 0 or 1 L p , p > 1 when L → +∞.We can now state the main theorem of this section.where ci 0 and ci are positive constants.

Asymptotic study of the signalling and dummy functions
The previous section presented numerous small results which give insight into the steady state behaviour of SRLK receptor-ligand systems.We are now in a position to combine these results to state and prove our main theorem, which gives closed-form formulae for the limits of the signalling and dummy functions.
and the limit of the dummy function is Under the assumption of no allostery, these expressions can be further simplified.
Corollary 4.15.Consider the system of Theorem 4.14 and assume there is no allostery.Denote the limiting components by Xi 1 , . . ., Xi r and Ni 0 ≡ Ni 1 = . . .= Ni r , their corresponding total number.The limit of the signalling function is and the limit of the dummy function is with z given by equation (43) in Lemma 4.10.
From the previous expressions we observe that the limit of the signalling and dummy functions are equal to the total copy number of the limiting component, Ni 0 , multiplied by a term which is bounded between 0 and 1.This term only depends on the affinity constant K0 and the steady state concentration of the kinase.In order to relate the above limits back to biologically meaningful quantities, all there is left to show is that the explicit expression of the limit of σ is in fact the amplitude of the system.Since σ(0) = 0, let us first note that the amplitude is equal to the maximum of σ.Under the no allostery assumption, we can show mathematically that this maximum is the limit of σ when L → +∞.To this end, the following lemma is needed.The supremum here is attained and is a maximum.Thus, the amplitude for a SRLK receptor-ligand system when there is no allostery is the limit of σ described in Corollary 4.15.This result is the generalisation of the example discussed in Section 3.1.1.We note that the amplitude of the IL-7R model of Section 3.1.1can be recovered by setting Ni 0 = min(Nx, Ny).We now have also rigorously shown that the limit of the signalling function is indeed the amplitude.The EC50 can now be found as outlined in Section 3.1.1.
(a) Decoy complexes with kinase.The sub-unit chain Y i can also bind to the intermediate dummy complexes X 1 : X 2 : . . .: X i , forming decoy complexes without kinase.(c) Scheme of the sequential formation of the signalling and dummy complexes.At each step their formation can be interrupted by the binding of a sub-unit chain, Y i , to the intermediate complex, forming a decoy complex.
Definition 4.18.An extended SRLK model is said to be under the assumption of no allostery if for each i = 1, . . ., n, Ki = K i and Ky i = K y i .
With these expanded definitions, we can extend the results previously obtained for the SRLK receptor-ligand systems.
Theorem 4.19.The theorems and lemmas previously true for the SRLK models are true for the extended SRLK models under the same (extended) hypotheses.

A few examples of (extended) SRLK models
In spite of some presumably strong modelling assumptions, the (extended) SRLK models can describe a broad range of cytokine-receptor systems.The IL-7R models described in Section 3.  A number of interleukin receptors share different molecular components.For instance, cytokine receptors of the common gamma family (comprising the receptors for IL-2,4,7,9,15 and 21 [9]) share the common gamma chain, γ.In addition the IL-2 and IL-15 receptors share the sub-unit chain, IL-2Rβ.The competition for these sub-unit chains can be mathematically described with an extended SRLK model, as follows.This example is illustrated in Figure 7a.(i) models the competition for IL-12Rβ1 between the IL-12 and the IL-23 receptors.We assume that IL-23R and IL-12Rβ2 are already bound to their associated extrinsic kinase JAK2; (ii) models the competition for IL-12Rβ2 between the IL-12 and IL-35 receptors.We consider the complexes IL-12Rβ1:TYK2 and gp130:JAK1 already pre-formed; (iii), models the competition for gp130 between the IL-35 and the IL-27 receptors.We consider the complexes IL-12Rβ2:JAK2 and IL-27R:JAK2 already pre-formed.
A further extended SRLK example is that of the IL-12 family of receptors, which share multiple components [55], and each of which is composed of two trans-membrane sub-unit chains.The IL-12 receptor is composed of the sub-unit chains IL-12Rβ1 and IL-12Rβ2.The IL-23 receptor signals via the IL-23R chain and the IL-12Rβ2 chain.The IL-27R (also known as WSX-1) and glycoprotein 130 (gp130) form the IL-27 receptor.Finally, IL-12Rβ2 and gp130 form the IL-35 receptor.The sub-unit chains gp130, IL-12Rβ1 and IL-12Rβ2 bind to a kinase from the JAK family (JAK1, TYK2 and JAK2 respectively).This competition can be described with extended SRLK models as follows.
Example 4.22 (Extended SRLK models: IL-12R family).We provide three examples of extended SRLK systems which characterise the competition for receptor sub-units between receptors of the IL-12 family (see Figure 7b).Above we have made use of the notation X * to denote the pre-formed complex composed of the receptor chain X and its intra-cellular extrinsic kinase (TYK2 for IL-12Rβ1, JAK1 for gp130 and JAK2 for all the others).

Conclusion
In the first part of this paper we propose a method to compute analytic expressions for two relevant pharmaco-dynamic metrics, the amplitude and the EC50 for receptor-ligand systems, based on two (simple) IL-7 receptor models.Our method starts with the computation of a Gröbner basis for the polynomial system of the receptor-ligand system in steady state.As shown in our IL-7R models from Section 3.1.1and Section 3.1.2,the derivation of the amplitude is easier when the maximum of the dose-response curve is attained at large ligand concentration (for instance when the dose-response curve is a sigmoid).In that case, the amplitude is the limit of the signalling function when the ligand concentration tends to infinity.When the model is simple enough, as is the case of the first IL-7R model, the polynomial system, simplified by the computation of the Gröbner basis, can be solved iteratively to obtain an analytic expression for the steady state.
From these expressions, it is then relatively straightforward to compute the amplitude (i.e., the limit of the signalling function at large values of the ligand concentration) and the EC50.For more complex models, such as our second IL-7R model, getting such steady state expressions can be more challenging.However, perturbation theory can be used to derive the expression for the amplitude.Computing another Gröbner basis can dramatically simplify the calculation of the EC50, and in turn display how it depends on the parameters of the model.Analytic expressions for the amplitude and the EC50 offer mechanistic insight for the receptor-ligand systems under consideration, allow one to quantify the parameter dependency of these two key variables, and can facilitate model validation and parameter exploration.Indeed, for both IL-7R models, we noticed that the affinity constant of the association of the gamma chain to the kinase JAK3, K1, was the only constant involved in the expression for the amplitude.As a consequence, and if conducting parameter inference to fit the model to experimental data, K1 would be the only parameter that could be inferred by comparison of the theoretical to the experimental amplitude.On the contrary, this constant was absent from both EC50 expressions and thus, its value would be impossible to infer by only comparing the experimental to the theoretical EC50.Our exact analysis also showed that both models have the same amplitude.Finally, the application of our method no longer requires the numerical computation of the dose-response curve, finding its maximum to then obtain the amplitude and fit the curve to derive the EC50.This reduces dramatically the computational cost and numerical errors.However, our method requires models simple enough to be able to compute a lex Gröbner basis, which is known to be computationally intensive [46,56].Additionally, computing the amplitude when the maximum response is not the asymptotic behaviour of the dose-response curve can be tricky.For instance the computation of the maximum for bell-shaped dose-response curves (which has been done for simple models in Refs.[35,36]) may involve the computation of the derivative of the signalling function.This computation can be laborious even with the use of symbolic software.Finally, our method often requires additional mathematical tools or knowledge, such as perturbation theory in Section 3.1.2,which makes it rather a challenge to be used by those who are not mathematically trained.In spite of the (sometimes, complicated) calculations that our method requires, we believe that analytic expressions of the pharmacological metrics characterising simple receptor-ligand systems may provide significant advantages when studying such biological systems.
In the second part of this paper, we introduced a family of receptor-ligand systems, called SLRK, in which the signalling complex, composed of a kinase, a ligand and n trans-membrane sub-unit chains, is built sequentially.These models could also form dummy and decoy complexes, similarly to the IL-7R models which the SRLK family encompasses.By manipulating the polynomials describing the SRLK models, we are able to derive an analytic expression of the amplitude under the no allostery assumption.We also show that the maximum of the dose response curve for both our IL-7R models was indeed the amplitude of the models.Despite relatively strong assumptions, we believe that the SRLK approach can be used to model a broad range of biochemical systems, such as receptor competition in interleukin signalling.The analytic expressions obtained for the amplitude could improve our understanding of biological mechanisms requiring a fine tuning of cytokine signalling such as cancer treatment [57] or cytokine storm control [58,59].We showed in Section 4.4 how our SRLK models can account for the competition for the gamma chain between the IL-2 family of receptors and the competition for receptor components between the IL-12 family of receptors.However, many receptors signal through different configurations.IL-35, for instance, can signal through homodimerisation of gp130 or IL-12Rβ2 [60].It has been shown that IL-6, a cytokine implied in cytokine storms [61,58], signals through an hexameric structure composed of two IL-6Rα chains and two gp130 molecules [62].Furthermore, it seems that the ligand IL-6 first binds to the IL-6Rα chain before any association with gp130 [62].Thus, one could imagine other general receptor models that may involve any of the following: 1) homo-oligomerisation (when two trans-membrane chains Xi are identical), 2) other orders of signalling complex formation (non-sequential orders or for instance, if the ligand is not the final sub-unit to be bound), 3) thermodynamic cycles (when there are several ways to form the signalling complex), 4) multiple kinases (including kinases binding to other sub-unit chains, such as JAK1 which binds to IL-7Rα [12]), or 5) a more detailed JAK-STAT pathway (most cytokine receptors activate multiple STAT molecules, whose copy numbers tune the immune response elicited [52]).
With this paper we hope to have initiated, or renewed, an interest for the algebraic analysis of receptor-ligand systems.Finally, we believe the results presented in this paper are a first step to account for the variability of receptor expression levels when designing and studying receptor-ligand models (both from an experimental and mathematical perspective) [5,15,8,7].

A Perturbation theory
A well known difficulty with the lex Göbner basis method, and polynomial equations in general, is that there is usually no analytic solution when the degree of a univariate polynomial is greater than four.This result is known as the Abel-Ruffini theorem [63].Therefore, in order to make progress, we need to resort to either numerical computations or analytic approximations.Since receptor-ligand systems are often characterised by a sigmoidal dose-response curve, at least to calculate the amplitude, the only quantity of interest is the limit of the signalling function at infinity.In order to calculate this limit (where possible analytically, otherwise numerically) we make use of perturbation theory for polynomial equations.
Greatly inspired by the Dover book written by Simmonds and Mann [54], this section reviews some notions of perturbation theory and justifies the steps of the method used to compute the analytic amplitude expression in Section 3.1.2.We start by defining an asymptotic expansion.The core of perturbation theory is the notion of asymptotic expansion and the following fundamental theorem.We are now ready to study the behaviour of the root of a univariate polynomial.Let n ∈ N * , where N * is the set of natural numbers without zero.We consider a univariate polynomial, P (x), of degree n, in the variable x, with coefficients which depend on the parameter .We are interested in the behaviour of the roots of P (x) when → 0. This polynomial can be re-written in the following form P (x) = (1 + b0 + c0 2 + . ..) + A1 α 1 (1 + b1 + c1 2 + . ..)x + . . .+ An αn (1 + bn + cn 2 + . ..)x n , where for each i αi is a rational number, bi, ci, . . .are real constants, (1 + bi + . ..) is a regular asymptotic expansion of the general from a0 + a1 + . . .+ aN N + O( N +1 ).
For such a polynomial, P (x), we have the following result.
Theorem A.3.Each root of a polynomial, such as equation ( 45) is of the form x( ) = p ω( ), ω(0) = 0, where ω is a continuous function of for sufficiently small and p ∈ Q.
The proof of this theorem (see Ref. [54]) gives a method to study the asymptotic behaviour of the roots of polynomial (45).
Proof.The concentration of kinase z being an algebraic function bounded on R between 0 and Nz, it admits a finite limit cz when L → +∞.We know that cz ≥ 0 because z is a concentration.We now prove that cz > 0. Since δ converges to a positive constant when L → +∞, we must have where ∆z = (1 + K1(N1 − Nz)) 2 + 4K1Nz.
Proof.We assumed no allostery so Ki = K i for all i = 1, . . ., n. Equation (42a) gives: The expression z1 is always positive while z2 is always negative.Hence z1 is the steady state kinase concentration, z.
Lemma 4.11.Let us consider an SRLK model which satisfies the experimental hypotheses.Then no concentration xi behaves proportionally to L q , q > 0 or 1 L p , p > 1 when L → +∞.
Proof.Lemma 4.9 affirms that z tends to a positive constant when L → +∞.In order for σ or δ to converge to a positive constant as stated in lemma 4.8, we need where c is a positive constant.Since the concentrations x1, . . ., xn are bounded functions (between 0 and their respective Ni), it is impossible to have for any i = 1 . . .n, xi ∼ L→+∞ ciL q with ci constant and q > 0. Since all the xi, with i = i1, i = i2, tend to a positive constant when L → +∞, we have where C is a positive constant.Thus, we obtain the behaviour of the left side of equation ( 63) Since the right side is given by it results in p1 = p2 = 1 2 .If there are r limiting components xi 1 , . . ., xi r , i1 < . . .< ir, then we have with ci 1 , . . ., ci r positive constants, and p1, . . ., pr > 0 such that p1 + . . .+ pr = 1.We proceed the same way as for the case of two limiting components and we obtain p1 = . . .= pr = 1 r .Ki .

D.2 Asymptotic study of the signalling and dummy functions
Using the limit properties and since everything converges, we obtain: Consequently since z → cz > 0 from lemma 4.9, we obtain: We substitute this limit into the expression of σ and δ and obtain the desired expressions.

Figure 1 :
Figure 1: Sigmoid dose-response curve: number of signalling complexes formed, σ(L), as a function of the concentration of ligand L (arbitrary units).A) The maximum value defines the amplitude.The EC 50 is the concentration of ligand which corresponds to half the amplitude.B) Three dose-response curves with the same EC 50 value and different amplitudes.Increasing the amplitude shifts up the maximum of the curve and increases the efficacy.C) Three dose-response curves with the same amplitude and different EC 50 values.Decreasing the EC 50 shifts the dose-response curve to the left and increases the potency of the ligand.

Figure 2 :
Figure 2: First IL-7R model: (a) The IL-7 receptor is composed of the trans-membrane γ and α chains.The γ chain can bind the intra-cellular downstream kinase JAK3.When the ligand, IL-7, binds the full receptor, it phosphorylates STAT5.(b) The IL-7R model allows the formation of "dummy" complexes: IL-7 bound IL-7R complexes, devoid of JAK3, which are unable to induce intra-cellular signalling.(c) IL-7 bound IL-7R complexes with JAK3 are able to induce intra-cellular signalling, and thus, are called "signalling" complexes.IL-7R dummy and signalling complexes are formed sequentially.The mathematical notation used in this paper is shown below each molecule or complex.
a) Decoy complex with kinase.Second IL-7R model: sequential chemical reaction scheme.

Figure 3 :
Figure 3: IL-7R model with an additional receptor sub-unit.The signalling and dummy complexes are the same as in the first IL-7R model.This second model allows the formation of decoy complexes: (a) with the kinase JAK3, or (b) without the kinase.(c) The IL-7R dummy and signalling complexes are formed sequentially.Decoy complexes can be formed to prevent the formation of signalling or dummy complexes.The mathematical notation used is annotated below each molecule or complex.

Figure 4 :
Figure 4: The lines defined in set E and the proper values, (black dots), computed following the graphical algorithm described in Appendix A.

Theorem 4 . 14 .
Consider an SRLK model which satisfies the experimental hypotheses.Let us write Xi 1 , . . ., Xi r as the limiting components and Ni 0 ≡ Ni 1 = . . .= Ni r as their corresponding total number.The limit of the signalling function is given by

Figure 6 :
Figure 6: Extended SRLK model with n = 4 trans-membrane chains: (a) An additional sub-unit chain, Y i , can bind to each intermediate signalling complex Z : X 1 : . . .: X i , to form decoy complexes with kinase.(b)The sub-unit chain Y i can also bind to the intermediate dummy complexes X 1 : X 2 : . . .: X i , forming decoy complexes without kinase.(c) Scheme of the sequential formation of the signalling and dummy complexes.At each step their formation can be interrupted by the binding of a sub-unit chain, Y i , to the intermediate complex, forming a decoy complex.

1 . 1 and
Section 3.1.2are, respectively, an SRLK and an extended SRLK model.In this section, we provide examples of other interleukin-signalling systems which are part of the SRLK family.

Figure 7 :
Figure 7: (a) Illustration of example 4.21: IL-2R, IL-7R and IL-15R competing for the common gamma chain and IL-2Rβ.The IL-2R is composed of three sub-unit chains: the gamma chain, IL-2Rβ and IL-2Rα.IL-15R is composed of the gamma chain, IL-2Rβ and the specific chain IL-15Rα.The IL-7R is composed of the gamma chain and IL-7Rα.All these receptors signal through the Janus kinase JAK3.(b) Illustration of example 4.22:(i) models the competition for IL-12Rβ1 between the IL-12 and the IL-23 receptors.We assume that IL-23R and IL-12Rβ2 are already bound to their associated extrinsic kinase JAK2; (ii) models the competition for IL-12Rβ2 between the IL-12 and IL-35 receptors.We consider the complexes IL-12Rβ1:TYK2 and gp130:JAK1 already pre-formed; (iii), models the competition for gp130 between the IL-35 and the IL-27 receptors.We consider the complexes IL-12Rβ2:JAK2 and IL-27R:JAK2 already pre-formed.

Lemma 4 . 10 .
where c d is a positive constant.Since σ also admits a finite limit when L → +∞, it means that z cs is a positive constant.So z has to satisfyz ∼ L→+∞ cz,where cz = cs c d is a positive constant.Consider an SRLK model which satisfies the experimental hypotheses.If we assume no allostery, then the steady state value of z is given by

L p for p > 1 .
From equation (61) it follows that it is impossible to have any xi ∼ L→+∞ c i Theorem 4.12.We consider an SRLK model which satisfies the experimental hypotheses.If there exists a unique limiting component Xi 0 then xi 0 ∼ L→+∞ ci 0 L , and for all i = 1, . . ., n, i = i0, xi ∼ L→+∞ ci,where ci 0 and ci are positive constants.