**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]** ---

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.

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
*p59 ^{fyn}* and

**
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.
**

**
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.

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.

**
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:

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).

**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.
**

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*.

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 *p _{o}={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).

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

**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
p_{0}=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 *p _{0}={0,1}*. One can see that for
density

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

**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, p_{0}=0;
(B) x=0.45, p_{0}=1;
(C) x=0.65, p_{0}=0;
(D) x=0.65, p_{0}=1;
(E) x=0.85, p_{0}=0; and
(F) x=0.85, p_{0}=1.
**

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
*p59 ^{fyn}* and

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).