Autophosphorylation and the dynamics of the activation of Lck
AAutophosphorylation and the dynamics of theactivation of Lck
Lisa Maria KreusserDepartment for Applied Mathematics and Theoretical PhysicsUniversity of CambridgeWilberforce RoadCambridge CB3 0WA, UKandAlan D. RendallInstitut f¨ur MathematikJohannes Gutenberg-Universit¨atStaudingerweg 9D-55099 MainzGermany
Abstract
Lck (lymphocyte-specific protein tyrosine kinase) is an enzyme whichplays a number of important roles in the function of immune cells. Itbelongs to the Src family of kinases which are known to undergo au-tophosphorylation. It turns out that this leads to a remarkable variety ofdynamical behaviour which can occur during their activation. We provethat in the presence of autophosphorylation one phenomenon, bistability,already occurs in a mathematical model for a protein with a single phos-phorylation site. We further show that a certain model of Lck exhibitsoscillations. Finally we discuss the relations of these results to modelsin the literature which involve Lck and describe specific biological pro-cesses, such as the early stages of T cell activation and the stimulation ofT cell responses resulting from the suppression of PD-1 signalling whichis important in immune checkpoint therapy for cancer.
Phosphorylation and dephosphorylation, the processes in which proteins aremodified by the addition or removal of phosphate groups, play an importantrole in biology. The activity of an enzyme is influenced by its phosphorylationstate and these processes provide a way of switching enzymes on and off quickly.The enzymes which catalyse phosphorylation and dephosphorylation are called1 a r X i v : . [ q - b i o . M N ] J a n inases and phosphatases, respectively. The phosphorylation of a protein X isusually catalysed by another protein Y. It may also be catalysed by X itself,a process called autophosphorylation. This can happen either in trans (onemolecule of X catalyses the phosphorylation of a site on another molecule of X)or in cis (a molecule of X catalyses the phosphorylation of a site on that samemolecule). Here we are concerned with the kinase Lck [5], which can undergoboth autophosphorylation in trans and phosphorylation by another kinase Csk.Lck belongs to the Src family of kinases [26] which have many properties incommon, in particular those related to their phosphorylation.In what follows we are interested in understanding the way in which theactivity of Lck is controlled, an issue which is important for analysing how thefunction of immune cells is regulated. More specifically, we want to do so bystudying mathematical models for phosphorylation processes. There has beena lot of work on models for cases where there is a clear distinction betweensubstrates and enzymes. A standard example is the multiple futile cycle wherebounds for the maximal number of steady states were obtained in [29] and [8]and for the maximal number of stable steady states in [7]. Much less is knownin the case of autophosphorylation. To our knowledge the earliest papers onmathematical modelling of Src family kinases are by Fuß et al. [10], [11]. In thefirst of these papers the authors consider a system coupling Src (with autophos-phorylation included) to Csk and the phosphatase PTP α . They then introducea simplification by assuming the concentration of Csk to be constant, and findtwo fold bifurcations in simulations. In particular, this system appears to ex-hibit bistability. In [11] sustained oscillations and infinite period bifurcationswere observed in a slight extension of the model of [10]. These dynamical fea-tures occurred in a context where the basic system describing phosphorylationand dephosphorylation of Src is embedded in feedback loops. In fact it wasfound in [14] that complicated dynamical behaviour is possible even withoutthe feedback loops. More recently the dynamics of a model for autophospho-rylation of a protein with only one phosphorylation site was studied in [6]. Inthat case also two fold bifurcations were observed. The model considered thereis one-dimensional and thus relatively easy to analyse. The bistability found in[6] contrasts with the situation in the multiple futile cycle where in the case ofa single phosphorylation site there is only one steady state.In Section 2 a model for autophosphorylation is introduced which is of centralimportance in what follows and it is shown that in a certain Michaelis-Mentenlimit it can be reduced to a one-dimensional model. Section 3 contains ananalysis of some properties of solutions of this reduced model. In particular itis shown that this system can exhibit more than one stable steady state. Thissection provides a rigorous treatment of some features found in the simulationsof [6]. The property of bistability is lifted to the original model. The main resultsare Theorems 1-3. The model of Section 2 without external kinase only exhibitsbistability under the condition that phosphorylation has an activating effect onthe enzyme. The corresponding case with inhibition exhibits no multistability.The aim of Section 4 is to show that in the case of an inhibitory phosphorylationmultistability can be restored by modelling the external kinase explicitly. The2ain result is Theorem 4. Here, in contrast to the results of Section 3, themultistability is not present in the Michaelis-Menten limit.Section 5 is concerned with a model for Lck which can be reduced bytimescale separation to a two-dimensional one. The original model inheritscertain patterns of behaviour such as bistability, Hopf bifurcations and homo-clinic orbits from the two-dimensional one. It is proved that the two-dimensionalmodel does exhibit these phenomena as a consequence of the occurrence of aBogdanov-Takens bifurcation. The main result is Theorem 5. In Section 6 themodels analysed in the present paper are compared with ones which occur asparts of more comprehensive models in the literature describing some concretebiological situations. Section 7 presents some ideas on possible further develop-ments of the results of this paper. Consider a protein with one phosphorylation site. We denote the unphospho-rylated form of this protein by X and the phosphorylated form by Y. SupposeX is able to catalyse its own phosphorylation in trans . The simplest model forthis reaction is 2X → X+Y. If Y is also able to catalyse the phosphorylationof X then this can be modelled by the reaction X+Y → x , y , e , f , d and c , respectively. The evolution equations are of the form˙ x = − k x + k c − k ex + k d − k xy, (1)˙ d = k ex − ( k + k ) d, (2)˙ e = − k ex + ( k + k ) d, (3)˙ c = k f y − ( k + k ) c, (4)˙ f = − k f y + ( k + k ) c, (5)˙ y = k x − k f y + k c + k d + k xy, (6)where the dot stands for the derivative with respect to t and the k i are positivereaction constants. There are three conserved quantities defined by the totalamounts of the substrate and the two enzymes E and F . These are A = x + c + d + y , B = c + f and C = d + e . A situation where the amountsof both enzymes and the rates of both autophosphorylation reactions are smallcan be described using a Michaelis-Menten reduction. To do this, introduce newvariables by means of the relations k = (cid:15) ˜ k , k = (cid:15) ˜ k , c = (cid:15) ˜ c , f = (cid:15) ˜ f , d = (cid:15) ˜ d , e = (cid:15) ˜ e and τ = (cid:15)t . Substituting these relations into the above equations and3ropping the tildes gives x (cid:48) = − k x + k c − k ex + k d − k xy, (7) y (cid:48) = k x − k f y + k c + k d + k xy, (8) (cid:15)d (cid:48) = k ex − ( k + k ) d, (9) (cid:15)e (cid:48) = − k ex + ( k + k ) d, (10) (cid:15)c (cid:48) = k f y − ( k + k ) c, (11) (cid:15)f (cid:48) = − k f y + ( k + k ) c, (12)where the prime stands for the derivative with respect to τ . If we set (cid:15) = 0in these equations the last four become algebraic. Combining these with theconservation laws and doing the usual algebra for Michaelis-Menten reductionleads to the relations c = ByK M + y and d = CxK M + x where K M = k + k k and K M = k + k k . It follows that x (cid:48) = − k x + Bk yK M + y − Ck xK M + x − k xy (13)while y satisfies an analogous equation. These two equations are equivalentbecause x + y is a conserved quantity for (cid:15) = 0. Thus the whole dynamics iscontained in the single equation (13) in that case. When C = 0 (no externalkinase) the equation for y reduces (up to a difference of notation) to the equation(1) in [6]. To make it clear that this is an equation for a single unknown it isnecessary to use the conserved quantity A = x + y . Thus for C = 0 the evolutionequation for y is y (cid:48) = k ( A − y ) + k ( A − y ) y − Bk yK M + y . (14) In [6] the authors describe certain aspects of the dynamics of solutions of equa-tion (14). Here we complement their analysis by giving rigorous proofs of someof these. Steady states of this equation are zeroes of the polynomial p ( y ) = [( k − k ) y + ( − k A + k A ) y + k A ]( K M + y ) − Bk y = − ( α − y + [ − K M ( α −
1) + A ( α − y + [ K M A ( α −
2) + A − Bk − k ] y + K M A (15)where α = k k . Positive steady states of the evolution equations for x and y are in one to one correspondence with roots of this polynomial in the interval(0 , A ). Note that p (0) > p ( A ) <
0. If k − k > p must haveone root with x < x > A . Thus it has exactly one root in thebiologically relevant region. When k − k < , A ). Since no root can cross the endpoints of the interval the number of4oots counting multiplicity is odd for any values of the parameters. In biologicalterms, bistability is only possible when phosphorylation activates the enzyme.In the case of Lck there are two phosphorylation sites of central importance forthe regulation of the kinase activity, Y394 and Y505, whose phosphorylationis activatory and inhibitory, respectively. Thus if we wanted to use this modelto describe Lck with mutations targeting one of its phosphorylation sites thento have a chance of bistability it is the inhibitory site Y505 which should beknocked out. This type of modification of Lck has been studied experimentallyin [3]. It was discovered that the mutated protein exhibits carcinogenic effects.It will now be shown that there is a region in parameter space where threepositive steady states exist. Theorem 1. If α > , A k < Bk , k is sufficiently large and K M is suf-ficiently small for fixed values of the other parameters then the equation (14)has three hyperbolic steady states, of which two are asymptotically stable and theother unstable. Proof.
When three steady states exist they must be simple zeros of p and itfollows that when ordered by the value of y the first and third steady states arestable while the second is unstable. Each of these steady states is hyperbolic.Thus to complete the proof of the theorem it suffices to prove the existence ofthree steady states under the given assumptions. The condition for a steadystate can be written in the form q ( y ) = q ( y ) where q ( y ) = ( A − y )[ k A +( − k + k ) y ] and q ( y ) = Bk yK M + y . Note that q ( y ) < Bk for all y ≥
0. If α > q has a local maximum when y = y = ( α − A α − . Assume that α > y >
0. Evaluating at the maximum gives q ( y ) = k αA α − . By choosing k large enough while keeping all other parameters fixed we can ensure that thismaximum is greater that Bk . It follows that q ( y ) > q ( y ). Then choosing K M small enough while keeping all other parameters fixed and using the factthat A k < Bk ensures that there is some y < y with q ( y ) < q ( y ).This implies that there are two roots of p which are less than y and these aresimple. Under these conditions p has three positive roots in the interval (0 , A )and so there exist three positive steady states. (cid:4) This theorem and its proof are illustrated by Fig. 1 where we show q , q and p for parameters A = 1, B = 2, k = 1, k = 1, k = 8, α = k k and K M = 0 . Lemma 1.
The polynomial p ( x ) = ax + bx + cx + d has a triple root if andonly if b = 27 a d and c = 27 ad . Proof. If x ∗ is a triple root then p ( x ∗ ) = p (cid:48) ( x ∗ ) = p (cid:48)(cid:48) ( x ∗ ) = 0. From the last ofthese equations we can conclude that x ∗ = − b a . Substituting this in the othertwo equations gives b = 3 ac and b = a ( bc − ad ). Combining the last twoequations gives b = 27 a d and c = 27 ad . Suppose conversely that b = 27 a d and c = 27 ad . Then bc = 9 ad and abc = 9 a d . Thus ( abc − a d ) = 27 a d b = a ( bc − ad ). Using b = 27 a d then implies that b = 3 ac . With all this information it can be checked directly that x ∗ = − b a isa triple root of p . (cid:4) Theorem 2.
The three steady states in Theorem 1 arise in a generic cuspbifurcation.
Proof.
To prove this it will be shown that the parameters can be chosen so thatthe polynomial p satisfies the conditions of Lemma 1. Assume that K M < A .Since b = [ − K M ( α −
1) + A ( α − we see that b = ( A − K M ) α + . . . for α large and b = ( K M − A ) + . . . for α →
0. Now A − K M > K M − A <
0. Thus if we consider b as a function of α with the otherparameters fixed it is an increasing function which takes on all values in theinterval [( K M − A ) , ∞ ). On the other hand a d = ( α − K M A and so a d = K M A α + . . . for α large and a d = K M A + . . . for α →
0. It followsthat there exists an α ∗ for which b = 27 a d . In this way the first conditionof Lemma 1 has been achieved. Since a d is non-negative there the same mustbe true of b and it follows that α ∗ >
2. Hence ad is negative and so in orderto achieve the second condition of the Lemma 1 it is enough to show that c can be given any prescribed negative value by choosing k appropriately whilefixing the other parameters. Note that a , b and d do not depend on k so thatthe first condition remains satisfied. Since α ∗ > c is positive for α = α ∗ and k sufficiently small. By increasing k it can then be made to haveany desired negative value. Thus it can be ensured that the second condition issatisfied. Note that the point x ∗ at which the bifurcation takes place does liein the biologically relevant region (0 , A ) since there is one steady state in thatregion and x ∗ is, neglecting multiplicity, the only one.Next we note that the derivative of the mapping ( α, K M , A, k ) (cid:55)→ ( a, b, c, d )6s always invertible for α >
2. Thus by the inverse function theorem we seethat by varying the parameters arbitrarily we can vary the coefficients of thepolynomial p arbitrarily in a neighbourhood of the values for the triple root.Thus we can choose two parameters so that the point with the triple root isembedded in a generic cusp bifurcation as defined in [16]. More specifically wecan choose a mapping ( β , β ) (cid:55)→ ( α, K M , A, k ) such that, after translatingthe coordinate y so that the bifurcation is at the origin, we have ( a, b, c, d ) =(1 , , β , β ). (cid:4) Consider now the rescaled mass action system (7)-(12) in the case C = 0.In this case we can discard the equations for d and e . Moreover, we can usethe conservation laws to discard the equations for x and c and replace thesequantities in the right hand sides of the equations for y and f . The result is y (cid:48) = k ( A − B − y + f ) − k f y + k ( B − f ) + k xy, (16) (cid:15)f (cid:48) = − k f y + ( k + k )( B − f ) . (17)We now want study the limit (cid:15) → Theorem 3.
There is a choice of parameters such that the system (1)-(6) with d = e = 0 , C = 0 and fixed values of A and B imposed has three steady states, ofwhich two are asymptotically stable and the other a hyperbolic saddle. The threesteady states arise in a generic cusp bifurcation. For arbitrary values of theparameters each solution converges to a steady state as t → ∞ . In particular,this system has no periodic solutions. Proof.
It suffices to prove corresponding results for the system (16)-(17). Thetheorem can be proved using the results of Theorems 1 and 2 and geometric sin-gular perturbation theory (GSPT) [15]. The important condition to be checkedis that of normal hyperbolicity. It says that on the critical manifold, which is thezero set of the right hand side in the equation for f (cid:48) , the derivative of that righthand side with respect to f should be non-zero. This is indeed the case since thederivative is − k f − k − k <
0. It can be concluded that for each hyperbolicsteady state of the Michaelis-Menten system there is a nearby steady state ofthe mass action system which is hyperbolic within the invariant manifold ofconstant A and B . In addition, when the steady state of the Michaelis-Mentensystem is stable the same is true of the corresponding steady state of the massaction system and when the steady state of the Michaelis-Menten system is un-stable the steady state of the mass action system is a saddle point whose stablemanifold is one-dimensional. To obtain the statement about the convergenceof general solutions to steady states we compute the linearization of (16)-(17)which is A = (cid:20) k ( f + y − A − B ) − k f k ( f + y − A − B ) − k y − k − (cid:15) − k f − (cid:15) − [ k y + ( k + k )] (cid:21) (18)It is always the case that A − y and B − f are positive on the region of biologicalinterest. Thus the system is competitive. Every solution of a competitive two7imensional system converges to a steady state [27] and this completes the proofof the theorem. (cid:4) To conclude this section we consider the limiting case of the system (14)obtained by setting k = 0. In this case only the phosphorylated form of theprotein is catalytically active. Bistability for a system of this type was consid-ered in [18]. If we continue to assume C = 0 then y = 0 is a steady state. Thusin order to get bistability we need to include that boundary steady state in thecounting. With this understanding we obtain an analogue of Theorem 1 for thiscase, where the condition on α is absent. The proof is strictly analogous to thatof Theorem 1. To see what happens to Theorem 2 in this case we need to replace p , which was got by division by k , by ˜ p = y [ k ( − y + A )( y + K M ) − Bk ]. Thispolynomial has a triple root at the origin when A = K M and AK M = Bk . We next consider the case where the phosphorylated kinase is completely inac-tive, which can be modelled by setting k = 0 in the model of the last section.This might be thought of as a model of the mutant of Lck where the activatorysite Y394 is knocked out. It should, however, be noted that in reality the cat-alytic activity of this mutant, although much reduced, is not actually zero [28].In that case we have k − k > k − k > C > k x + k x ( A − x ) + Ck xK M + x = Bk ( A − x ) K M + A − x . (19)Since the function on the left hand side of this equation is monotone increasingon [0 , A ] and is zero for x = 0 while the function on the right hand side ismonotone decreasing on [0 , A ] and is zero for x = A these two functions areequal at a unique point x ∈ (0 , A ). Thus there cannot be more than one steadystate in the Michaelis-Menten system with k < k . It turns out, however,that there can be more than one steady state in the corresponding mass actionsystem, even in the case k = 0.Positive solutions of the mass action system with k = 0 are in one to onecorrespondence with solutions of the following system obtained by using theconserved quantities to eliminate d , c and y .˙ x = − k x + k ( B − f ) − k ex + k ( C − e ) , (20)˙ e = − k ex + ( k + k )( C − e ) , (21)˙ f = − k f ( A − B − C − x + e + f ) + ( k + k )( B − f ) . (22)8efine a polynomial by p ( x ) = (cid:80) i =0 a i x i with coefficients a = k k k , (23) a = 2 k k k ( k + k ) + k k k k , (24) a = k k [ − k (( A + B ) k − ( k + 2 k ) C − k ( k + k ))+ k ( k + k ) + k k ( k + k )] − k k ( k + k ) k , (25) a = k k ( k + k ) { k [ − A + B ) k + ( k + 2 k ) C ] − k + k ) k } + k k [ k ( Bk − k C ) + k ( k + k ) ] , (26) a = k [ − k k ( A + B )( k + k ) − k k ( k + k ) B + k ( k A − ( k + k ) C − k ( k + k ))( Bk − k C )] − ( k + k ) k [ k ( k + k ) + k k C ] , (27) a = k k k ( k + k )[ B ( Ak − ( k + k ) C − k ( k + k ))+ A ( k B − k C )] − ( k + k ) k k ( k + k ) k C, (28) a = k k ( k + k ) AB. (29)Define x max to be the largest value of x satisfying the inequalities k k x + k k Cxk ( k x + k + k ) ≤ B, (30) x + k k x + (cid:18) k k (cid:19) k Cx ( k x + k + k ) ≤ A. (31)Note that x max depends continuously on the parameters. Lemma 2.
For given positive values of A , B and C positive steady state solu-tions of the system (1)-(6) with k = 0 are in one to one correspondence withroots of the polynomial p in the interval (0 , x max ) . Proof.
Note first that the equations for steady states of (1)-(6) are equivalent tothe equations for steady states of (20)-(22) and that these in turn are equivalentto the equations k ex = ( k + k )( C − e ) , (32) k f ( A − B − C − x + e + f ) = ( k + k )( B − f ) , (33) k x = k ( B − f ) − k ( C − e ) . (34)Now suppose that ( x, d, e, c, f, y ) is a positive steady state. It follows from(32) that e = ( k + k ) Ck x + k + k and combining this with (34) gives f = B − ( k /k ) x − k k Cxk ( k x + k + k ) . Thus we have solved for e and f in terms of x . Substituting thisinformation into (33) and rearranging gives the equation p ( x ) = 0. Supposeconversely that x is a root of p with 0 < x < x max . Define e = ( k + k ) Ck x + k + k , (35) f = B − ( k /k ) x − ( k /k )( C − e ) , (36) y = A − B − C − x + e + f. (37)9t follows from (30) that the quantity f defined by (36) is positive and from(31) that the quantity y defined by (37) is positive. It follows directly that (32)and (34) hold. The fact that x is a root of p implies that (33) holds and hencethat ( x, d, e, c, f, y ) is a positive solution of (1)-(6). (cid:4) Consider now the real roots of p . They depend continuously on the param-eters and their number is constant modulo two. Since p (0) > p cannot pass through zero. We claim that x max can also never be a root of p .For if x were equal to x max while satisfying the inequalities (31) then at leastone of them would become an equality. In the first case f as defined by (36)would be equal to zero. This contradicts equation (33). In the second case y defined by equation (37) would be equal to zero. But then it follows from (4)that c = 0 and from (6) that x = 0, a contradiction. It can be concluded thatthe sign of p ( x max ) is independent of the parameters. To determine what thesign is it suffices to evaluate it for some particular values of the parameters.Choose k i = 1 for all i , A = 5, B = 2 and C = 3. When x = 1 we see thatequality holds in (30) while the strict inequality holds in (31). Thus in thiscase x max = 1. Evaluating the coefficients in p gives a = 1, a = 5, a = 8, a = − a = − a = − a = 40. Hence p ( x max ) = − <
0. Itfollows from the intermediate value theorem that p has a least one root in eachof the intervals (0 , x max ) and ( x max , ∞ ). The number of sign changes of thecoefficients in the polynomial is even and at most four. Thus Descartes’ rule ofsigns implies that the number of positive roots is zero, two or four. The casewith no positive roots has already been ruled out. Thus there are two or fourand at least one of them must be greater that x max . With the parameter valuesin the example there are only two changes of sign, only two positive roots andwe know that precisely one is less than x max . By continuity the number of rootsin (0 , x max ) counting multiplicity is odd. It can only be one or three and wehave already seen an example of parameters where it is one. In that case thesystem (1)-(6) admits precisely one positive steady state.We will show that there also exist parameter values such that the systemhas three positive steady states. One approach would be to show that thereare parameters for which there is a triple root in the desired interval and thenperturb. In fact we will show directly that there are parameters for which thereare three roots in that interval, since that approach is simpler. Theorem 4.
There is a choice of parameters for which the system (1)-(6) with k = 0 has three positive steady states. Proof.
Due to Lemma 2 it suffices to find parameter values for which thepolynomial has three roots in the interval (0 , x max ). It follows from the precedingdiscussion that it is enough to show that the interval contains at least two roots.It turns out that is suffices to choose A = 6, B = 20, C = 2, k i = 1 for all i (cid:54) = 5and k sufficiently large. With these choices it follows that we get the followingasymptotics for k → ∞ . a = k + . . . , a = k + . . . , a = − k + . . .a = 18 k + . . . , a = − k + . . . where the terms not written explicitly are o ( k ), as are the coefficients a and a . It follows that k − x − p ( x ) = q ( x )+ o (1),where q ( x ) = x + x − x + 18 x −
4. In the limit the inequalities defining the10dmissible interval become x <
18 and x + x <
2. By continuity it suffices toshow that q has two roots in the interval (0 , q (0) < q (1 / > q (1) < (cid:4) The results of the previous sections were related to situations in which one of thetwo key regulatory phosphorylation sites in a Src family kinase such as Lck ismutated. In the present section we move to the case where both sites are present.The starting point for the discussion is the model introduced in [14]. Therefour phosphorylation states of the kinase are included in the description. Thefirst, denoted by S i , is that where the inhibitory site is phosphorylated while theactivatory site is not. This form of the kinase shows no catalytic activity. S , S a and S a are the forms where neither site is phosphorylated, only the activatorysite is phosphorylated and both sites are phosphorylated, respectively. All ofthese are catalytically active to some extent and can catalyse the transition S → S a . The transitions S → S i and S a → S a are catalysed by Csk.The transitions S i → S and S a → S a are catalysed by one phosphatase andthe transitions S a → S and S a → S i are catalysed by another phosphatase.Experimental results obtained in [12] indicate that some modifications of theseassumptions may be needed to obtain a biologically correct model. In particular,it was found that Y505 in Lck undergoes autophosphorylation in trans , albeitwith a much lower rate than Y394. The variants of the model of [14] whichwould be needed to take this into account will not be considered further in thepresent paper - the aim here is rather to see the variety of dynamical behaviourwhich this type of system can produce.Let us introduce the following neutral notation for the quantities involvedin the model, denoting the concentrations of S , S i , S a and S a by x , x , x and x respectively. Then X = x + x + x + x is a conserved quantity. Thereaction rate for the autophosphorylation is bilinear, the dephosphorylation of S a is given by Michaelis-Menten kinetics and the other reactions are assumedto be linear. Using the notations of [14] for the reaction constants gives thesystem ˙ x = − k x + k x + k x β + x − k x ( δx + x + x ) , (38)˙ x = k x − k x + k x , (39)˙ x = k x ( δx + x + x ) − k x β + x + k x − k x , (40)˙ x = k x − ( k + k ) x . (41)Before considering this system in the general case note that setting x , x ,11 and k to zero reduces this system to˙ x = k x β + x − k x ( δx + x ) , (42)˙ x = k x ( δx + x ) − k x β + x . (43)Either of the variables can be eliminated using the conserved quantity givingan equation which is, up to a difference in notation, exactly the equation ofDoherty et al. discussed in previous sections. Only setting k and k to zero in(38)-(41) gives a partially decoupled system which is the product of the systemof [6] with a hyperbolic saddle. It follows immediately from Theorem 1 thatfor suitable values of the parameters the system (38)-(41) admits at least threepositive steady states, of which two are stable and hyperbolic and the third isa hyperbolic saddle.We now return to the general system (38)-(41). In [14] the authors find aremarkable variety of dynamic behaviour in the system above which, after fixinga value of the conserved quantity, is of dimension three. They remark that thereis a limiting case which gives rise to a system of dimension two which alreadyexhibits a lot of this dynamics. To investigate this possibility we define a newvariable by y = x + x and use it to replace x . In addition we introduce rescaledparameters satisfying ˜ k = (cid:15)k and ˜ k = (cid:15)k . Making these substitutions anddiscarding the tildes leads to the system˙ x = − k x + k x + k y − x β + y − x − k x ( δx + y ) , (44)˙ x = k x − k x + k x , (45)˙ y = k x ( δx + y ) − k y − x β + y − x − k x , (46) (cid:15) ˙ x = k ( y − x ) − ( k + (cid:15)k ) x . (47)This is a fast-slow system with one fast and three slow variables. We have theconserved quantity X = x + x + y . In the limiting case (cid:15) = 0 equation (47)reduces to y − x = ξx , where ξ = k k . It follows that x = ξ y . Substitutingthis into (44) and (45) and using the conserved quantity gives the followingsystem of two equations:˙ x = − k x + k x + k ξ ( X − x − x ) β ( ξ + 1) + ξ ( X − x − x ) − k x [( X − x − x ) + δx ] , (48)˙ x = k x − k x + k ξ + 1 ( X − x − x ) . (49)In the terminology of GSPT this is the restriction of the system to thecritical manifold. This critical manifold is normally hyperbolic and stable sincethe partial derivative of the right hand side of (47) with respect to x evaluatedat (cid:15) = 0 is negative. This allows us to transport information about stability and12ifurcations from steady states of the two-dimensional system to steady statesof the full system. Positive steady states of the full system with a given valueof X are in one to one correspondence with positive steady states of (48)-(49)with x + x < X . At steady state equation (49) can be used to express x interms of x and substituting this into (48) shows that x is a root of a cubicpolynomial which is not identically zero. Thus the system (48)-(49) has at mostthree steady states.It will be shown that the system (48)-(49) admits periodic solutions whicharise in a Hopf bifurcation and homoclinic orbits. In order to do this it sufficesto show that this system admits a generic Bogdanov-Takens bifurcation [16].By saying that the bifurcation is generic we mean that it satisfies the condi-tions BT.0, BT.1, BT.2 and BT.3 of [16]. Then the desired results follow fromTheorem 8.5 of [16] and the analysis of the normal form of the bifurcation pre-ceding that theorem. Let J ( x , x ) be the linearization of the system (48)-(49)about the point ( x , x ). Finding a bifurcation point where the condition BT.0is satisfied means finding a point ( x , x ) and a choice of the parameters of thesystem so that J ( x , x ) has a double zero eigenvalue but is not itself zero. Ifthe right hand sides of equations (48) and (49) are denoted by f and f thismeans solving the system of four equations given by the vanishing of f , f ,tr J and det J . The general strategy is to choose all but four of the parametersand use the four equations to solve for the rest. An obstacle to this is that thequantities resulting from this process might fail to be positive. This obstaclewas overcome by trial and error.The equations for steady states can be written in the following form. k = [ k x − k x + k x ( X − x − x + δx )][ β ( ξ + 1) + ξ ( X − x − x )] ξ ( X − x − x ) , (50) k = ( ξ + 1)( − k x + k x ) X − x − x . (51)The linearization is J = (cid:20) − k − φ ( β ) − k ( X − x − x + 2 δx ) k − φ ( β ) + k x η − ω (cid:21) (52)Here we have introduced the auxiliary quantities η = k − k ξ +1 , ω = k + k ξ +1 and φ ( β ) = k ξ ( ξ +1) β [ β ( ξ +1)+ ξ ( X − x − x )] . Suppose that we have a Bogdanov-Takensbifurcation. Since the trace is zero we have that the first element in the firstrow of the Jacobian must be equal to ω . Hence φ ( β ) + ω + k = − k [ X − − δ ) x − x ] . (53)It follows that k − φ ( β ) + k x = ω + k + k + k [ X − (1 − δ ) x − x ] . (54)13ince the determinant is zero we have ω + η [ ω + k + k + k ( X − x − x + 2 δx )] = 0 . (55)Choose X = , x = 1, x = , k = 8, k = 1, δ = , ξ = 1. It follows from(51) that k = 8. Putting this into (55) gives k = . It then follows from(53) that φ ( β ) = . On the other hand (50) implies that k = (8 β + 1).Combining this with the definition of φ shows that β = . Finally we compute k = . The conclusion is that with the given choices there is exactlyone solution for the remaining parameters ( k , k , β, k ) such that the systemsatisfies the condition BT.0 for a Bogdanov-Takens bifurcation at the chosenpoint with coordinates ( x , x ) = (cid:0) , (cid:1) . At this point the linearization is of theform J = (cid:20)
12 48 − − (cid:21) . (56)When talking about a Bogdanov-Takens bifurcation we need a system dependingon two parameters. In our example we choose these to be δ and k and considerall other parameters in the system as fixedAs will now be explained, a calculation shows that conditions BT.1, BT.2and BT.3 are also satisfied so that this is a generic Bogdanov-Takens bifurcation.For this purpose it is convenient to transform to coordinates y = − x + and y = x + 4 x − J . Then ˙ y i = ( J y ) i + Q i ( y ) + O ( | y | ), where J is in Jordan form and the Q i are quadratic. In the notationof [16] the elements of Q and Q are denoted by a ij and b ij , respectively. Inthe present example it turns out that a = 0 and in that case BT.1 and BT.2are the conditions that b (cid:54) = 0 and b (cid:54) = 0. A lengthy calculation shows that b = − W + 168 k > b = − W + 7 k >
0, where W = βk (8 β +1) . Herewe use the values of the parameters at the bifurcation point. The condition BT.3is that the linearization J T of the mapping ( x , x , δ, k ) (cid:55)→ ( f , f , tr J, det J )at the bifurcation point is non-singular. This matrix is J T =
12 48 − k − − −
12 0 0 − W + k − W + k − k W − k W − k k − . (57)and det J T = − k − W k (cid:54) = 0.When a generic Bogdanov-Takens bifurcation is present in a dynamical sys-tem then there are always generic Hopf bifurcations nearby. The periodic solu-tions which arise in these Hopf bifurcations are hyperbolic and may be stable(supercritical case) or unstable (subcritical case). We may correspondingly callthe Bogdanov-Takens bifurcation super- or subcritical and it turns out thatthese two cases are distinguished by the relative sign of b and b . In thepresent case the signs of these two coefficients are equal and the bifurcation issubcritical. Hence the periodic solutions are unstable. In comparison with thephase portrait given in [16], which corresponds to the supercritical case, the14irection of the flow is reversed. For the parameter values for which the bifur-cation takes place the cubic polynomial for x has a double root at x = 1 andmust therefore have a factor ( x − . Carrying out this factorization allows athird root to be calculated explicitly. The result is p ( x ) = 3728 (2457 x − x − . (58)The additional root is x = and at the corresponding steady state x = . At that point the trace of the linearization is negative and the determinantpositive. Hence this steady state is stable.Some of these results will now be collected in a theorem. Theorem 5.
There are parameter values for which the system (48)-(49) has ageneric Bogdanov-Takens bifurcation. In particular, there are nearby parametervalues for which it has an unstable periodic solution and ones for which it hasa homoclinic orbit. In the case where there is an unstable periodic solution withparameter values sufficiently close to those at the bifurcation point there are alsotwo stable steady states and one saddle point.
The structural stability of the bifurcation and the fact that the limit usedto obtain this system is normally hyperbolic implies that these features can belifted to the system (38)-(41). In more detail, note first that this system isequivalent by rescaling to the system (44)-(47). Moreover we can concentrateon a fixed value of the conserved quantity X . Thus it remains to consider alimit from a three-dimensional system to a two-dimensional one. Restrictingto the slow manifold we get a regular limit of two-dimensional systems. For (cid:15) = 0 the mapping ( x , x , δ, k ) (cid:55)→ ( f , f , tr J, det J ) has full rank and a zeroat the bifurcation point. It follows from the implicit function theorem that ithas a unique zero near the bifurcation point for (cid:15) small. This is a point whereBT.0 is satisfied. By continuity BT.1, BT.2 and BT.3 remain satisfied for (cid:15) sufficiently small and the bifurcation remains subcritical. Thus the featureslisted in Theorem 5 are also seen in the system on the slow manifold. Thisimplies immediately that there is a heteroclinic orbit in the full system. Thehyperbolic periodic solution in the slow manifold is also hyperbolic as a solutionof the full system. It is of saddle type with there being both solutions whichconverge to it for t → + ∞ and solutions which converge to it for t → −∞ . Ifit could be shown that the limiting system admits a stable periodic solution forsome values of the parameters it could be concluded that the full system doesso too. We have not been able to prove the existence of stable periodic solutionsof the limiting system. To see how such a stable solution might occur, considerthe predator-prey model of Bazykin discussed in [16]. It has two subcriticalBogdanov-Takens bifurcations and a stable periodic solution in a distant partof the phase space.It turns out that it is possible to find an extension of the explicit Bogdanov-Takens point to an explicit two-parameter family of steady states, includinga one-parameter family of points where the eigenvalue condition for a Hopf15ifurcation is satisfied. The parameters are the determinant σ and the trace τ of J at the steady state. This family is obtained by fixing the same quantitiesas in the original case, including the coordinates ( x , x ) of the steady state andcomputing the parameters k = 324 + 4 σ + 36 τ , β = 132 + 5 σ + 24 τ σ + 1248 τ (59)and k = (384 + 5 σ + 45 τ )(1536 + 20 σ + 180 τ )21(1404 + 15 σ + 156 τ ) . (60)The Hopf points are those with σ > τ = 0. These formulae are helpfulin finding parameter values for which the system admits an unstable periodicsolution. A solution of this type is illustrated in Fig. 2 in red for the case σ = 1, τ = − .
02. In Fig. 2(a) we show neighbouring solutions of the unstable periodicsolution which move away from it (inward and outward spirals). A larger partof the phase space is shown in Fig. 2(b) where the Bogdanov-Takens point andthe stable steady state are shown. (a) Unstable periodic solution. (b) All steady states.
Figure 2: Phase diagram for σ = 1, τ = − . In this section the results of this paper will be put into a wider context bycomparing the models analysed here with some more complicated models in theliterature which arise when studying specific biological phenomena. One of themost exciting recent developments in medicine are immune checkpoint therapiesfor cancer [22]. Immune checkpoint molecules such as CTLA4 and PD-1 canresult in the deactivation of T cells under certain circumstances and this isexploited by cancer cells to evade attacks by the immune system. Antibodies tothe immune checkpoint molecules can prevent this and thus be used in cancertherapies. This type of therapy has had remarkable success in curing somecancers. On the other hand, although these therapies could in principle work16or all types of cancer, in practice they only work for some cancers, notablymetastatic melanoma, and even in the most favourable cases for only a certainpercentage of patients. It is important to obtain a better understanding of themolecular mechanisms of these therapies, so as to explain in which cases theyare effective and, hopefully, to improve them so as to increase the range of theirefficacy.Up to now the most effective types of immune checkpoint monotherapy arethose involving PD-1. In this context it is important to understand how theactivation of PD-1 leads to suppression of T cell activity. This has been studiedexperimentally in [13]. The main conclusion of that work is that the inhibitionof T cell activity caused by PD-1 is due less to a decrease in signalling via theT cell receptor than to a decrease of the second signal coming from CD28. Thisleads to a certain mechanistic model of how the influence of PD-1 is exerted.In an effort to obtain a better understanding of the mechanisms involved amathematical model was introduced in [4]. Simulations of that model gaveresults agreeing well with the results of [13] and at the same time suggestingan additional path by which PD-1 can influence T cell signalling. In the pathhighlighted in [13] activation of Lck plays an important role. The suggestionin [4] is that this change in the activation state of Lck could have an indirectinfluence via phosphorylation of molecules downstream of the T cell receptorand CD28. The model of [4] consists of several modules. One of these describesthe activation of Lck and plays a central role.In the context of their model of Lck regulation the authors of [4] cite amodel given in [23]. The latter includes complexes which are intermediatesin the autophosphorylation reactions. This would correspond in the case withone phosphorylation site to replacing the reaction 2X → X+Y by the reactions2X → X → X+Y, where X is the complex formed when two molecules of X bindto each other. It also includes certain complexes of Lck with Csk which areanalogous to XF in the basic model introduced in section 2. According to[23] the inclusion of these complexes was necessary to obtain a good agreementbetween the results of simulations and the experimental data of [12]. The modelof [23] includes no phosphatases and so it is clear that in that case the evolutionmust converge to the state where only the unique maximally phosphorylatedstate is present. The non-trivial characteristics of the evolution have to do withthe way in which the solution approaches that state.One difference of the model of [4] compared to that of [14] is that it includesfive forms of Lck rather than four. The model contains two different forms ofdoubly phosphorylated Lck which are supposed to differ by the order in whichthe two sites were phosphorylated. The issue of the order of phosphorylationis mentioned in [12] but we are not aware of any justification for including thefifth form in the model. It is stated in [4] that the model includes autophos-phorylation of Lck but in the equations the dependence on the concentrationsof the different forms of Lck is everywhere linear and this does not seem to beconsistent.In [25] the author discusses the model of [14] and the alternative where theMichaelis-Menten term in that model is replaced by a linear one. When that17implification is made bistability is eliminated. The author presents unpub-lished data of Acuto and Nika which addresses the issue of bistability in Lckexperimentally. The idea is that if bistability was present the distribution ofthe measurements of certain quantities in a population should be bimodal, i.e.the graph should exhibit two maxima. In these data most (but not all) of thegraphs have a unique maximum and this is taken as evidence that there is nobistability in the system. However no detailed justification for this conclusionis given. The significance of this conclusion is that if bistability were presentin the biological system this would mean that the simplified model would notbe sufficient. In [25] the simplified model is used. The advantage is that thereare less parameters in the simplified model and that their values can be morestrongly constrained by experimental data.Another biological phenomenon where Lck plays a central role is that ofT cell activation. It will now be discussed how Lck has been modelled in theliterature on that subject. One of the first and most important steps in Tcell activation is the phosphorylation of the ITAMs (immunoreceptor tyrosine-based activation motifs) of the T cell receptor complex. The most importantkinase carrying out this process is Lck. In one successful model of early T cellactivation [9] Lck is not one of the chemical species included in the model. Inthe process of ITAM phosphorylation Lck is treated as an external kinase whoseactivity is represented by a reaction constant. It was proved in [21] that thismodel can exhibit more than one steady state. The model of [9] is a radicalsimplification of a more extensive one introduced in [2]. In the latter modelactivated Lck is one of the chemical species included. It takes part in manyreactions where it binds to a complex X containing the T cell receptor and someother molecules and then phosphorylates some element of the complex. Thekinetics of these reactions is extended Michaelis-Menten. Other forms of Lckplay a role in mechanisms represented in this model but they do not appearexplicitly. Another model implementing some of the same mechanisms waspresented in [17]. It includes four forms of Lck arising from phosphorylation atY394 and the serine S59. The serine phosphorylation may have an importantrole to play in T cell activation but will not be considered further here. Thetyrosine phosphorylation is supposed to occur by autophosphorylation in trans but the Lck molecules responsible for the catalysis are supposed to belong to adifferent population to those being phosphorylated. The former population istreated as external and so no nonlinearity arises from this process. The modelof [17] exhibits bistability. In this paper we proved that the model of [6] of an enzyme with a single sitesubject to autophosphorylation in trans can exhibit bistability. This improveson the simulations in [6] showing this type of behaviour for specific parametervalues by identifying a large part of parameter space where it occurs. We alsoshow that in the context of this model multiple steady states can only occur18hen the phosphorylation increases the activity of the enzyme. It is shown thatin a case where phosphorylation decreases the activity of the enzyme multiplesteady states can also occur but this requires a more complicated model with anexternal kinase which is operating well away from the Michaelis-Menten limit.We related the models studied in this paper to other models involving Lckwhich have been applied in the literature to describe particular biological phe-nomena. It is of interest to consider the possible biological meaning of the resultsof this paper. Switches arising through bistability are a well-known phenomenonin biology and the bistability found in the regulation of Lck might be of impor-tance for immunology as a mechanism by which the activity of immune cells isswitched off or on in certain circumstances. As discussed in the last section itseems unclear on the basis of experimental evidence whether bistability due tothe properties of Lck occurs in biologically interesting circumstances. We arenot aware that oscillations in the activation of Lck have been observed experi-mentally. The biological significance of those periodic solutions whose existencewe proved is limited by their instability.This paper is a preliminary exploration of dynamical features of models forLck involving autophosphorylation. At this point it is appropriate to thinkabout what biological issues could be illuminated by continuing these investiga-tions. A question of great biological and medical interest, already mentioned inthe last section, is that of the mechanism by which ligation of the receptor PD-1leads to the suppression of the activity of T cells. (For a recent review of thistopic see [20].) Normally the activation of a T cell requires both a signal fromthe T cell receptor and a second signal from CD28. Both of these receptors getphosphorylated. (In the case of the T cell receptor it is rather the associatedproteins CD3 and the ζ -chain which are phosphorylated). A question which isapparently still controversial is whether the main effect of PD-1 activation isdesphosphorylation of the T cell receptor or that of CD28. The conclusion of[13] is that it is CD28 but this has been disputed in [19], where it has beensuggested that this finding of [13] was an artefact of using a cell-free system andthat in reality it is dephosphorylation of the T cell receptor which is the mostimportant consequence of the activation of PD-1. This indicates that betterunderstanding of these phenomena is necessary. The authors of [4] claim thattheir model can reproduce the results of [13]. Could that model, or a relatedone, reproduce the results of [19]?There is a wide consensus that, whatever the targets of dephosphorylationresulting from the activation of PD-1, the phosphatase which carries it out isSHP-2. There is one caveat here since it was observed in [24] that PD-1 canhave an inhibitory effect on T cells in the absence of SHP-2. This issue de-serves further investigation. Another interesting question is that of the way inwhich Lck interacts with SHP-2 and PD-1. When PD-1 is fully activated itis phosphorylated at two sites. These provide binding sites for SHP-2. Thephosphorylation of PD-1 is catalysed primarily by Lck [13]. SHP-2 can dephos-phorylate PD-1 and thus promote its own unbinding. This effect is opposedby Lck. Here there is an incoherent feed-forward loop [1]. On the one handLck causes phosphorylation of PD-1 by a direct route and on the other hand19t causes its dephosphorylation by an indirect route. These interactions are de-scribed by one of the modules in the model of [4]. They are sufficiently complexthat it would be desirable to carry out a deeper mathematical analysis of theirdynamics. Acknowledgments
LMK acknowledges support from the European Union Horizon 2020 researchand innovation programmes under the Marie Sk(cid:32)lodowska-Curie grant agreementNo. 777826 (NoMADS), the Cantab Capital Institute for the Mathematics ofInformation and Magdalene College, Cambridge (Nevile Research Fellowship).
References [1] Alon, U. 2006 An introduction to systems biology: design principles ofbiological circuits. Chapman and Hall, Boca Raton.[2] Altan-Bonnet, G. and Germain, R. N. 2005 Modelling T cell antigen dis-crimination based on feedback control of digital ERK responses. PLoS Biol.3(11), e356.[3] Amrein, K. E. and Sefton, B. M. 1988 Mutation of a site of tyrosine phos-phorylation in the lymphocyte-specific tyrosine protein kinase, p56lck, re-veals its oncogenic potential in fibroblasts. Proc. Natl. Acad. Sci. USA 85,4247–4251.[4] Arulraj, T. and Barik, D. 2018 Mathematical modeling identifies Lck asa potential mediator for PD-1 induced inhibition of early TCR signaling.PLoS ONE 13:e0206232.[5] Bommhardt, U., Schraven, B. and Simeoni, L. 2019 Beyond TCR signaling:emerging functions of Lck in cancer and immunotherapy. Int. J. Mol. Sci.20, 3500.[6] Doherty, K., Meere, M. and Piiroinen, P. T. 2015 Some mathematical mod-els of intermolecular phosphorylation. J. Theor. Biol. 370, 27–38.[7] Feliu, E, Rendall, A. D. and Wiuf, C. 2020 A proof of unlimited multista-bility for phosphorylation cycles. Nonlinearity 33, 5629–5658.[8] Flockerzi, D., Holstein, K. and Conradi, C. 2014 n -site phosphorylationsystems with 2 n −−