Poster CN0314

Lattice Simulations of Protein-Tyrosine Kinase 
Autophosphorylation Reactions

Petteri Kettunen
University of Oulu, Department of Physical Sciences, Division of Biophysics
Linnanmaa, FIN-90570 Oulu, Finland & Department of Chemistry
West Virginia University, Morgantown, WV 26506-6045, U.S.A.

Marcelle Kaufman
Centre for Nonlinear Phenomena and Complex Systems
C.P. 231, Boulevard du Triomphe, B-1050 Bruxelles, Belgium
--- Internet address: [email protected] ---


Autocatalytic phosphorylation reactions are involved in the early steps of helper T cell activation. A model of autocatalytic phosphorylation-dephosphorylation reactions is simulated on a two-dimensional lattice. The model consists of two autophosphorylation reactions and an enzyme-catalyzed dephosphorylation reaction. Simulations are performed with particle diffusion implemented as a nearest-neighbor random walk. It is shown that nonlinear features, such as phase transitions and bistability, are not obtained with low diffusion, but are recovered by increasing the diffusion rate. Therefore a theoretical model of T-cell unresponsiveness that relies on a memory effect due to bistability is shown to remain relevant when spatial fluctuations of the concentrations are allowed.


Our computational study is based on a theoretical model of T-cell unresponsiveness, which considers bistability at the level of protein tyrosine kinase activity. In this mean-field model the phosphorylation-dephosphorylation reactions are independent of space. Here, the reactions are studied with probabilistic lattice simulations. It is shown that bistability can be obtained with realistic short-range interactions and with relatively high particle diffusion. Thus, the predictions of the mean-field model are shown to be relevant for biological systems, for which the spatial dependency cannot be omitted.

Key words

Autophosphorylation, bistability, enzyme kinetics, immunology, lattice simulations, Michaelis-Menten reaction, protein-tyrosine kinase, T-cell.

Title | Introduction | Materials and Methods | Results | Conclusion | Comments | Sessions


T cell activation is initiated by the aggregation of cell surface receptors (Fig. 1) with the activation of receptor-associated, autophosphorylating, protein tyrosine kinase (PTK) enzymes as one of the first consequences (Metzger, 1992; Rudd et al., 1994). Two plausible canditates for an increased tyrosine kinase activity following TRC aggregation are the src-like tyrosine kinases p59fyn and p56lck, which have been shown to be physically associated with the TCR complex (Veillette et al., 1988). Autophosphorylation markedly enhances the phosphorylative activity of these enzymes and represents an important autocatalytic step in the early stages of transmembrane cell signaling (Cooper and Howell, 1993). A theoretical model of density-dependent PTK phosphorylation-dephosphorylation reactions has been developed, whose mean-field equations have, in some parameter regions, given bistable solutions and irreversible transitions between two regimes of phosphorylative activity (Kaufman et al., 1992, 1996). It has been proposed that residual kinase activity after antigen removal may be linked to an inability to respond normally to further stimuli by the same antigen (Cho et al., 1993; Gajewski et al., 1994; Quill et al., 1992). This unresponsive state, which is also referred to as T-cell anergy (Schwartz, 1990), is not permanent but reversible, as shown in Fig. 2.

FIG1 Click to enlarge

Fig. 1: Schematic view of receptor density increase upon TCR-antigen interaction. (A) In the absence of antigen T-cell expresses a low and homogeneous TCR-distribution. (B) Physical contact between APC and T-cell induces TCR aggregation on the cell-cell contact area. (C) Planar view of T-cell surface before (left) and after (right) TCR aggregation.

FIG2 Click to enlarge

Fig. 2: Flow diagram of the induction of reversible unresponsiveness in T-cells.

Spatial constrains have been shown to play an important role in the dynamics of chemical reactions (Nicolis and Prigogine, 1987; Imbihl and Ertl, 1995). In particular, in systems where the reactive species are membrane bound but free to diffuse laterally the space-diffusion limitations might not allow all the particles to travel far enough and to interact with all the other particles in the system, as is tacitly assumed in a mean-field approach (Nicolis and Prigogine, 1987). In this work, a representative model of autophosphorylative surface reactions is simulated on a two-dimensional lattice, where the interactions between reactive species are spatially restricted. A numerical analysis is performed as a function of particle density and with particle transport implemented by a nearest-neighbor random walk. The model considers an enzyme-catalyzed dephosphorylation reaction, and a phenomenological simulation algorithm for a Michaelis-Menten type of reaction on a lattice is applied. The study is focused on nonlinear features of the steady state behavior such as phase transitions and bistability.

Title | Introduction | Materials and Methods | Results | Conclusion | Comments | Sessions
Materials and Methods

Protein tyrosine-kinase reaction model

The TCR-associated PTKs have a specific phosphorylation-dependent catalytic activity (Veillette and Fournel, 1990; Cooper and Howell, 1993). In our protein phosphorylation-dephosphorylation cascade shown in Fig. 3 (Kaufman et al., 1992, 1996; Kettunen and Kaufman, 1996), E and P represent two characteristic states of a prototype src-family PTK enzyme with different catalytic activities. E and P possess low and high catalytic activity, respectively. E is converted into P by intermolecular phosphorylation by another E and by P (Fig. 3. R1,R2). Since the phosphorylated species P has a higher phosphorylative activity than E, we set b>a. The dephosphorylation reaction from P to E in vivo is known to occur through intra-cellular phosphatases (Rudd et al., 1994). For simplicity, we do not decompose this reaction into its basic steps and assume that dephosphorylation follows globally a Michaelis-Menten scheme (Murray, 1989). Since the catalyzing enzyme is not consumed, the net enzymatic dephosphorylation is modeled by R3 (Fig. 3) where c=c(K,V) approximates the influence of all factors underlying a Michaelis-Menten type of reaction.

FIG3 Click to enlarge

Fig. 3: PTK phosphorylation-dephosphorylation cascade.

There are no fluxes of the reactive species, i.e. d[E]/dt+d[P]/dt=0 so that the total enzyme concentration obeys a conservation relation denoted by [E]+[P]=x. Therefore the mean-field rate equation for [P] can be written as:

d[P]/dt = a (x-[P])(x-[P]) +b [P] (x-[P]) -V [P]/(K+[P]),

where K and V are the characteristic constants of a Michaelis-Menten reaction. With appropriate parameter values (Kaufman et al., 1992, 1995; Kettunen and Kaufman, 1996), the stable fixed points of this reaction system form a familiar S-shaped curve as a function of total enzyme concentration x with one unstable and two stable steady states, as shown in Fig. 4 below. Low densities of x are characterized by dominating E-particles, while increasing x enhances the autocatalytic phosphorylation reactions (Fig. 3, R1, R2); consequently the steady state at high x is characterized by a relatively high [P] concentration. At intermediate x there is a domain of co-existence of two stable steady states. In this region different initial conditions may lead to different final states. Such bistability provides a mechanism for memory and could account for the induction of an unresponsive state in T-cell immunology (Kaufman et al., 1993, 1996).

FIG4 Click to enlarge

Fig. 4: Steady-state concentrations of E and P. Solid and dashed lines indicate stable and unstable steady states, respectively. Kinetic parameters: a=0.01; b=0.8; K=0.01; and V=0.02.

Lattice model and simulation algorithm

We consider the TCRs with their associated E- and P-enzymatic forms as hard-core particles distributed on a two-dimensional triangular lattice. Occupation of the lattice sites is subject to a strong exclusion principle: one lattice site can be occupied by one particle only. This implies realistic hard-core interactions between the TCRs. Any particle can interact only with its six nearest neighbors, which are located equidistantly. The initial particle density is determined by the surface density x, which describes the total enzyme concentration. An arbitrary site is occupied with probability x, and vacant with probability 1-x. The initial state of a particle is set to P with probability p0, otherwise it is set to E (Stauffer, 1985; Kaufman and Stauffer, 1993).

Using an appropriate scaling one can set all the concentration and parameter values between zero and unity. Thus, we can treat the concentrations and reaction rates as probabilities and the lattice simulations for reactions R1 and R2 (Fig. 3) are implemented using a probabilistic scheme (Grassberger, 1982; Provata et al., 1993). A phenomenological algorithm for simulating the enzyme-catalyzed reaction R3 (Fig. 3) has been presented in earlier studies (Kettunen and Kaufman, 1996) and is described here briefly. We propose that reaction rate c (Fig. 3, R3) is a function of the concentration [P] in an area within the range of interaction and model the enzymatic reaction as follows. Any isolated P-particle changes its state to E with a constant average probability c(0), a P-particle with one P-neighbor changes its state with probability c(1), and so on.

The reactions are updated in a completely asynchronous manner: only one reaction can take place during one discrete time step. Particle diffusion is implemented as a nearest-neighbor random walk using asynchronous updating. The schematic simulation algorithm is as follows (Provata et al., 1993; Tretyakov et al., 1995; Kettunen and Kaufman, 1996):

(A) A site i is chosen from the lattice at random. If site i is occupied by a P-particle, (B) is performed. If site i contains an E-particle, (C) is performed. If site i is empty, the algorithm continues from (F).

(B) Read the status of the six neighboring sites of site i. Let n be the number of P-particles found from the neighborhood. Draw a random number z from an uniform deviate (0<z<1). The P-particle on site i changes its state to E, if z<c(n) (Fig. 3, R3).

(C) A neighbor site j is chosen randomly from the six nearest neighbors, with an equal probability of 1/6 each. If site j contains a P-particle, the algorithm continues with (D). If the particle on site j is E, the algorithm continues with (E).

(D) A random number z is drawn; if z<b (Fig. 3, R2), the particle on site i changes its state from E to P.

(E) A random number z is drawn; if z<a (Fig. 3, R1), the particle on site i changes its state from E to P.

(F) A site k is chosen from the lattice at random. If site k is occupied, choose one of the six neighbor sites, say site l, with equal probability 1/6 for each. If site l is vacant, the particle from site k makes a one-step random walk to site l.

As we want to investigate systems in which the diffusion and reaction rates are favored differently, we apply step (F) N times after the reactions are updated (steps A-E), before starting a new time step with (A).

Title | Introduction | Materials and Methods | Results | Conclusion | Comments | Sessions


The main objective of the simulations has been given to investigate if the bistability features that are observed in the mean-field model (Figs. 3 and 4) are maintained when the reactions are simulated in a low-dimensional space. All simulations were carried out on 50 x 50 two-dimensional triangular lattices with periodic boundary conditions, and with parameters a=0.01, b=0.8, and c(n)={ 1.000, 0.154, 0.087, 0.059, 0.045, 0.037, 0.031 } (Kettunen and Kaufman, 1996). Short-range interactions for the chemical reactions and low diffusion does not yield bistability, as shown in Fig. 5. However, a phase transition to relatively high P-particle concentration is located at density 0.56<x<0.6. Initial conditions other than po={0,1} were also tested but they did not lead to bistability. Localized islands of highly active P-particle clusters could be observed at various density values close to the phase transition point (Fig. 6).

FIG5 Click to enlarge

Fig. 5: Simulated final concentrations of [P] as a function of x with low diffusion N=10. Initial conditions p0=0 (blue triangles) and p0=1 (red circles).

FIG6 Click to enlarge

Fig. 6: Spatial snapshot (25 x 25 subimage from 50 x 50 lattice) of a lattice in steady state with low diffusion N=10. Blue and red circles indicate E and P particles, respectively. Total concentration x=0.55 and initial condition p0=1.

Further simulations showed that bistability is recovered with increasing diffusion (Fig. 7). The domain of two stable steady-states is located at the density range 0.57<x<0.8. When compared with Fig. 4 it is apparent that the lattice model exhibits the same qualitative properties as the mean-field model.

Spatial snapshots of the lattices visualize the results more clearly. Fig. 8 shows particle distribution in the steady state at three density values x={0.45, 0.65, 085} and for two different initial conditions p0={0,1}. One can see that for density x=0.45 the ratio of E and P particles in Figs. 8A and 8B are the same and that E-particles are dominating in both cases. At intermediate density x=0.65, the initial condition p0=0 (low initial enzymatic activity) leads to a high E concentration in the steady state (Fig. 8C), while the initial condition p0=1 (high initial enzymatic activity) results in a high catalytic activity in the final state (Fig. 8D). At high particle density x=0.85 most of the particles become finally P independently from the initial conditions (Figs. 8E,F).

FIG7 Click to enlarge

Fig. 7: Simulated final concentrations of [P] as a function of x with high diffusion N=1000. Initial conditions p0=0 (blue triangles) and p0=1 (red circles).

FIG8 Click to enlarge

Fig. 8: Spatial snapshots (25 x 25 subimage from 50 x 50 lattice) of lattices at their steady state with high diffusion N=1000. Blue and red circles indicate E and P particles, respectively. Total concentrations and initial conditions: (A) x=0.45, p0=0; (B) x=0.45, p0=1; (C) x=0.65, p0=0; (D) x=0.65, p0=1; (E) x=0.85, p0=0; and (F) x=0.85, p0=1.

Title | Introduction | Materials and Methods | Results | Conclusion | Comments | Sessions


An earlier lattice model of an autocatalytic phosphorylation cascade was implemented as a spin system: an E-particle becomes P, if at least two of its nearest neighbors are P (Kaufman and Stauffer, 1993). That model did not include diffusion, but simulations yielded first-order transitions when random birth and death of particles was added. The lattice approach, which is presented in this study, gives a more realistic description of the phosphorylation-dephosphorylation reactions. From the biological point of view, simulating reactions in a low-dimensional space with short-range interactions is relevant with respect to the real biological situation. Indeed, PTKs such as p59fyn and p56lck are associated with the TCR which is itself embedded in the cell membrane (Veillette et al, 1988; Rudd et al, 1994). These PTKs have thus almost strictly two-dimensional translational mobility.

In conclusion, it is shown that a lattice model with short range interactions confirm the validity of results obtained in a mean-field theory, provided the diffusion rate is sufficiently high. It should be stressed that the interest of the present study goes well beyond the scope of cellular immunology, given that the role of spatial dimensionality on the behavior of chemical systems is actively being investigated in various fields (Imbihl and Ertl, 1995; Provata et al., 1993; Tretyakov et al., 1995).


Cho E.A., Riley M.P., Sillman A.L., and Quill H., J. Immunol. 151:20 (1993).
Cooper J.A. and Howell B., Cell 73:1051 (1993).
Gajewski T.F., Qian D., Fields P., and Fitch F.W., Proc. Natl. Acad. Sci. USA 91:38 (1994).
Grassberger P., Z. Phys. B Condensed Matter 47:365 (1982).
Kaufman M. and Stauffer D., J. Stat. Phys. 73:843 (1993).
Kaufman M., Andris F. and Leo O., Theoretical insight into antigen-induced T-cell unresponsiveness, in Theoretical and experinmental insight into immunology (Eds. Perelson A.S. and Weisbuch G.), NATO ASI series H66, Springer-Verlag, 93-115 (1992).
Kaufman M., Andris F. and Leo O., Int. Immunology 8:613 (1996).
Kettunen P. and Kaufman M., Math. Mod. and Sci. Comp., in press (1996).
Imbihl R. and Ertl G., Chem. Rev. 95:697 (1995).
Metzger H., J. Immun. 149:1477 (1992).
Murray J.D., Mathematical Biology, Springer-Verlag, Berlin (1989).
Nicolis G. and Prigogine I., Exploring Complexity, Freeman, New York (1987).
Provata A., Turner J.W., and Nicolis G., J. Stat. Phys. 70:1195 (1993).
Rudd C.E., Janssen O., Cai Y.-C., da Silva A.J., Raad M., and Prasad K.V.S., Immunol. Today 15:225 (1994).
Stauffer D., Introduction to Percolation Theory, Taylor & Francis, London (1985).
Schwartz R.H., Science 248:1349 (1990).
Tretyakov A., Provata A. and Nicolis G., J. Phys. Chem. 99:2770 (1995).
Quill H., Riley M.P., Cho E.A., Casnellie J.E., Reed J.C., and Torigoe T., J. Immunol. 149:2887 (1992).
Veillette A., Bookman M.A., Horak E.M., and Bolen J.B., Cell 55:301 (1988).
Veillette A. and Fournel M., Oncogene 5:1455 (1990).

Title | Introduction | Materials and Methods | Results | Conclusion | Comments | Sessions