help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Guo, C.
Right arrow Articles by Levine, H.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Guo, C.
Right arrow Articles by Levine, H.

Biophys J, November 1999, p. 2358-2365, Vol. 77, No. 5

A Thermodynamic Model for Receptor Clustering

Chinlin Guo and Herbert Levine

Department of Physics, University of California, San Diego, La Jolla, California 92093-0319 USA

    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

Intracellular signaling often arises from ligand-induced oligomerization of cell surface receptors. This oligomerization or clustering process is fundamentally a cooperative behavior between near-neighbor receptor molecules; the properties of this cooperative process clearly affect the signal transduction. Recent investigations have revealed the molecular basis of receptor-receptor interactions, but a simple theoretical framework for using these data to predict cluster formation has been lacking. Here, we propose a simple, coarse-grained, phenomenological model for ligand-modulated receptor interactions and discuss its equilibrium properties via mean-field theory. The existence of a first-order transition for this model has immediate implications for the robustness of the cellular signaling response.

    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

Cell growth, differentiation, migration, and apoptosis are regulated in part by extracellular polypeptide growth factors or cytokines (Heldin, 1995; Stuart and Jones, 1995). As these molecules are unable to pass through the hydrophobic cell membrane, they have to bind to the extracellular domains of specific surface receptors to exert their effects. Much effort has gone into investigating the fundamental question of how the ligand-receptor interaction can trigger the proper intracellular signals. One popular hypothesis is that ligand-induced "clustering" of ligand-receptor complexes can be a key element in the proper activation of downstream signals. (Ashkenazi and Dixit, 1998; Bray et al., 1998; Heldin, 1995; Germain, 1997; Lemmon and Schlessinger, 1994, 1998; Reich et al., 1997; Sakihama et al., 1995).

As an example of this line of reasoning, we consider the signaling cascade mediated by the binding of tumor necrosis factor (TNF) to the receptor TNF-R1. Internally, the cytoplasmic domain of TNF-R1 is "sensed" by a variety of adaptor proteins, namely TRADD, FADD, TRAF2, and RIP; this sensing leads eventually to NF-kappa B/JNK/SAPK activation and apoptosis. To accomplish the downstream signaling, an oligomerization of these adaptor proteins is required (Ashkenazi and Dixit, 1998). One way to facilitate oligomerization is via construction of a molecular scaffolding by TNF-induced TNF-R1 clustering. It is known that TNF-R1 will not aggregate in the absence of TNF; this is due to the association of an inhibitor, "silencer of death domain" (SODD), which normally attaches to TNF-R1 cytoplasmic domains and prevents receptor aggregation (Jiang et al., 1999), or, alternatively, is due to the receptor extracellular domains, inasmuch as spontaneous association of TNF-R1 has been observed in cells that express truncated receptors (Boldin et al., 1995; Vandevoorde et al., 1997). TNF treatment, however, can bring two or more receptors into proximity via its multiple binding capacity (Jones et al., 1990, 1992). This "proximity" might "squeeze" out SODD (Jiang et al., 1999), expose the cytoplasmic "death" domains to adaptor proteins, and thereby stabilize receptor clusters. Thus, a molecular scaffold/nuclei is generated to initiate signaling.

Over a longer time scale, the signaling messages can provide feedback to modify the capability of surface receptor clustering (Humphries, 1996; Wyszynski et al., 1997). This leads to a complex dynamical process involving both the intracellular signaling cascades as well as the surface receptor clustering. The self-organization made possible by these feedback processes has been intensively discussed for signaling cascades (see, e.g., Jafri and Keizer, 1995; Barkai and Leibler, 1997). Much less is understood, however, regarding the role of receptor clustering. It is clear, though, that given the hypothesis that cellular signaling relies on the formation of receptor clusters, the temporal and spatial characteristics of clustering would certainly affect the process of signaling transduction. Thus, modeling the physical properties of receptor clustering is as important as modeling signaling cascades.

Because clustering is due to an interaction between nearest-neighbor receptors, it is obviously a cooperative process. From a physics perspective a system with this type of cooperativity can exhibit a first-order phase transition, corresponding to a jump in the surface density of ligand-receptor complexes. In the coexistence region of this transition, the surface will spontaneously segregate into two phases, dilute and dense. This first-order phase transition endows the signal transduction process with the ability to produce a digital signal in an analog world; this is independent of the details of intracellular cascades, arising instead from the intrinsic cooperativity in ligand-receptor interaction. This has not been adequately addressed in the few models studied to date (Goldstein and Wiegel, 1983; Goldstein and Perelson, 1984; Riley et al., 1995; Coutsias et al., 1997; Shea et al., 1997).

The purpose of this work is to introduce a phenomenological model for the TNF-TNFR1 system to describe the onset of receptor clustering (phase separation). Specifically, we assume that clustering can be described by the statistical mechanics of a simple lattice Hamiltonian, incorporating the fundamental mechanism of a multimeric binding capacity for the ligand. We will calculate (via mean-field theory) a phase diagram and show that clustering will be thermodynamically favored for some range of ligand and receptor densities. Finally, we will do a simple Monte Carlo simulation of this system, showing that receptor diffusion will lead rapidly to cluster formation in the relevant parameter range. We neglect the possibility that there exist long-time feedback processes that modify the clustering capacity, and we ignore some inessential details of the receptor-ligand interaction. More detailed models including these effects, as well as applications to other signaling systems, will be presented in the future.

    THE LATTICE HAMILTONIAN
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

In our model, we treat the cell surface as a lattice with a spacing on the order of a few nm; this is the closest that neighboring receptors can get to each other. Each lattice site i has either one or zero receptor molecules, denoted as ni = 1 or 0. Our receptor has only two states: liganded or unliganded, and the interaction between receptor molecules is determined by their states. This "two-state" model is oversimplified, yet we will see that it gives reasonable predictions for the phase diagram. A "state" label, ti = 1 or 2, to represent unliganded or liganded, then, can be assigned to each occupied receptor. We will further assume that the only ligands on the surface are those bound to receptors. If we let the chemical potential of the ligand be µL and that of the receptor be µR, we then get a contribution to the effective Hamiltonian of the system
H<SUB>1</SUB>({n, t})=<UP>−</UP><LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> &mgr;(t<SUB><UP>i</UP></SUB>)n<SUB><UP>i</UP></SUB> (1)
where µ(1) = µR and µ(2) = µR + µL + gL and -gL is the binding energy between ligand and receptor.

We should clarify the relationship between the parameters used here and those in real experiments. Using standard ideas (Changeux et al., 1967), we notice that with only this term, the partition function can be factorized and reduced to a single site problem,
Z=<LIM><OP>∏</OP><LL><UP>i</UP></LL></LIM> Z<SUB><UP>i</UP></SUB>=<LIM><OP>∏</OP><LL><UP>i</UP></LL></LIM> <FENCE>1+<LIM><OP>∑</OP><LL><UP>k=1</UP></LL><UL><UP>2</UP></UL></LIM> e<SUP><UP>&bgr;&mgr;</UP>(<UP>k</UP>)</SUP></FENCE> (2)
From this, we can immediately obtain the expectation values of the TNF-R1 concentration in the liganded and unliganded states. These are assumed to correspond to the equilibrium condition of the following reaction (Corti et al., 1994; Grell et al., 1998): TNF-R1(m) + TNF right-left-harpoons  TNF · TNF-R1(m), with a corresponding equilibrium dissociation constant, ([TNF-R1(m)]eq[TNF]/[TNF · TNF-R1(m)]eq) = Kdtnf approx  0.59 nM, where the notation TNF-R1(m) means a TNF-R1 molecule distributed on the artificial membrane, and where the brackets [...]eq indicate the equilibrium concentration of the respective molecule. From this, we have ebeta L+gL) = [TNF]/Kdtnf. To obtain the parameters individually, we might employ an "ideal gas law" for the ligand. This yields ebeta µL = [TNF](h2/2pi mtnfkBT)3/2, and gL = kBT ln[(2pi mtnfkBT/h2)3/2/Kdtnfapprox  60kBT, where h is the Planck constant and mtnf is the mass of TNF.

We next add a receptor-receptor interaction term. This takes the general form
H<SUB>2</SUB>({n, t})=<UP>−</UP>½ <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>a(t<SUB><UP>i</UP></SUB>, t<SUB><UP>j</UP></SUB>)n<SUB><UP>i</UP></SUB>n<SUB><UP>j</UP></SUB> (3)
Here, Jij = 1 only when ij are nearest neighbors and is 0 otherwise (Fig. 1). The function a(ti, tj) indicates a "state"-dependent interaction energy between nearest-neighbor receptors, namely, a(1, 1) is the energy between two unliganded receptors, a(1, 2) a(2, 1) is the energy between one liganded and one unliganded receptor, and a(2, 2) is the energy between two liganded receptors. We note that in general, higher order terms might exist, especially considering the "trimeric" nature of the TNF ligand in our model problem. We have similarly neglected the details of the interactions of the cytoplasmic domains, per our earlier discussion. Our goal is to elucidate the basic idea regarding clustering in the simplest possible model, assured that adding more details will not change the basic notion that there exists a first-order transition due to the cooperativity.



View larger version (17K):
[in this window]
[in a new window]
 
FIGURE 1   Interleaved sublattices labeled as filled/open circles on a one-dimensional and a two-dimensional (square/honeycomb) space. On a square lattice, j1,2,3,4 is the nearest neighbor to site i.

a(1, 2) and a(2, 1) are the interaction energies, for which we will use an effective binding strength gE on the order of gL/10, arising via one or two hydrogen bonds between receptors. It is important to realize that our simplified model does not treat explicitly the formation of multimers via multimeric binding. Instead, it arbitrarily assigns the one ligand (e.g., binding two receptors into a dimer) to one of the receptors and describes the dimeric binding as an attraction between a bound and an unbound receptor. Because of this, the model cannot distinguish between this relatively strong interaction and the subsequent much weaker interaction between the dimers. In future work, we will show that this complication does not alter the basic picture presented here.

As discussed above, in the TNF system there is probably a short-range and nonspecific "excluding" interaction between two unliganded or two liganded (with different ligand molecules) receptors. For the sake of simplicity, we will assume that the repulsive energy is on the same order of magnitude as the associative one, i.e., a(1, 1) approx  a(2, 2) approx  -gE. This assumption is not necessary, yet it greatly simplifies the mathematical task for analysis.

The symmetry of a(ti, tj) allows us to introduce a simple matrix notation for the total Hamiltonian H1 + H2. If we use two-component vectors for the state labeling: tau i = [01] for ti = 1, and tau i = [10] for ti = 2, then the Hamiltonian can be rewritten as
H({n<SUB><UP>i</UP></SUB>, &tgr;<SUB><UP>i</UP></SUB>})=<UP>−</UP><LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> n<SUB><UP>i</UP></SUB>[&mgr;(1), &mgr;(2)]&tgr;<SUB><UP>i</UP></SUB>−½g<SUB><UP>E</UP></SUB> <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>n<SUB><UP>i</UP></SUB>n<SUB><UP>j</UP></SUB>&tgr;<SUP><UP>+</UP></SUP><SUB><UP>i</UP></SUB><FENCE><AR><R><C><UP>−</UP>1</C><C>1</C></R><R><C>1</C><C><UP>−</UP>1</C></R></AR></FENCE>&tgr;<SUB><UP>j</UP></SUB> (4)
Here [µ(1), µ(2)] is a 1 × 2 matrix. The simplicity of using this form of the matrix a(ti, tj) can immediately be seen if we make a transformation
&tgr;<SUB><UP>i</UP></SUB>=½<FENCE><AR><R><C>1+&sfgr;<SUB><UP>i</UP></SUB></C></R><R><C>1−&sfgr;<SUB><UP>i</UP></SUB></C></R></AR></FENCE>
with sigma i = ±1. Then
H({n<SUB><UP>i</UP></SUB>, &sfgr;<SUB><UP>i</UP></SUB>})=<UP>−</UP>x <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> n<SUB><UP>i</UP></SUB>−y <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> n<SUB><UP>i</UP></SUB>&sfgr;<SUB><UP>i</UP></SUB>+½g<SUB><UP>E</UP></SUB> <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>n<SUB><UP>i</UP></SUB>n<SUB><UP>j</UP></SUB>&sfgr;<SUB><UP>i</UP></SUB>&sfgr;<SUB><UP>j</UP></SUB> (5)
where x = [µ(1) + µ(2)]/2 is the "averaged" receptor chemical potential, and y = [µ(1) - µ(2)]/2 is directly related to the ligand concentration, e-beta y = <RAD><RCD>[TNF]<IT>/K</IT><SUB>d</SUB><SUP>tnf</SUP></RCD></RAD>. The partition function then reads
Z=<LIM><OP>∑</OP><LL>{<UP>n<SUB>i</SUB>, &sfgr;<SUB>i</SUB></UP>}</LL></LIM> <UP>exp</UP>[<UP>−</UP>&bgr;H({n, &sfgr;})] (6)
where Sigma {ni;sigma i} means ensemble summation over the three different configurations {ni = 0; ni = 1, sigma i = ±1} on each lattice site and beta  = 1/kBT, where kB is the Boltzmann factor and T is the temperature.

If we define a new notation ui = nisigma i, our model would be very similar to a spin-1 antiferromagnetic (AFM) BEG model (Blume et al., 1971),
H({u<SUB><UP>i</UP></SUB>})=<UP>−</UP>x <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> u<SUP><UP>2</UP></SUP><SUB><UP>i</UP></SUB>−y <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> u<SUB><UP>i</UP></SUB>+½g<SUB><UP>E</UP></SUB> <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>u<SUB><UP>i</UP></SUB>u<SUB><UP>j</UP></SUB> (7)
The origin of this AFM behavior is the "negative cooperation" between nearest-neighbor receptors, as we have imposed that a "proximity" of two unliganded or two liganded receptors will cost energy. Similar behavior might occur in the erythropoietin receptor (EPO-R) and the human growth hormone receptor (hGH-R) systems (Heldin, 1995). This negative cooperation will give rise to an absence of clustering in extreme high/low ligand concentration (i.e., y right-arrow ±infinity ) and thereby result in a "bell" shape or window-like signaling response (Elliott et al., 1996).

We should point out that this negative cooperation is not universal. In the case of an EGF-R (epidermal growth factor receptor) system, a ferromagnetic (FM) behavior ("positive cooperation") is more likely, because there clustering requires two or more liganded receptors (Lemmon et al., 1997). Thus the higher the ligand concentration, the more the EGF-R cluster can be formed, and the EGF-EGFR signaling response behaves in a sigmoidal rather than a window-like pattern. It is clear that in both EGF-R and TNF-R (and hGH-R, EPO-R) systems, the ligand multiple binding capacity is the essential ingredient for inducing clustering (of course one should consider the effect of the receptor cytoplasmic domain as well). Which kind of cooperation (negative or positive) one should consider depends on the details of the receptor-receptor interaction (also including the chemical modifications on receptor cytoplasmic domains) and needs to be established experimentally. But, the essential feature of a first-order transition-like behavior in receptor clustering is independent of the sign of this additional cooperativity.

    NUMERICAL SIMULATION
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

To see if our model can generate clustering, we perform a Monte Carlo simulation on a square lattice with the standard Metropolis scheme. For simplicity, we fix the number of liganded and unliganded receptors and do not allow these to fluctuate. Given the rather strong binding, this is not an important constraint. Furthermore, we allow motion only for individual receptors and do not explicitly allow a cluster to move as a whole; this might not be the case in reality. The "jumping" probability that a receptor will move to another lattice site is determined by the Hamiltonian and obeys the detailed balance law. In detail, we pick a receptor at random and try to move it in a randomly chosen direction. The move is accepted if it lowers the energy, and the move is accepted with probability e-beta Delta H if the energy increases.

From Fig. 2, we immediately see that for a given receptor density, changing the ligand concentration moves the system from a nonclustering to a clustering phase. In this figure, the open and filled circles indicate liganded and unliganded receptor molecules, respectively. Note that the open and filled circles are arranged in an alternative way to form the cluster (i.e., inside a cluster, the nearest neighbors of the open circles must be filled circles, and vice versa). This implies that the equilibrium state (which must be translationally invariant) can be described by dividing the system into two interleaved sublattice systems: one sublattice is occupied by one species of receptor molecule (liganded or unliganded), and all of its nearest neighbors belong to the alternative sublattice, which is occupied by another species.



View larger version (59K):
[in this window]
[in a new window]
 
FIGURE 2   Monte Carlo simulation with Metropolis scheme. Here we test the model under a fixed receptor density but different ligand concentrations. In both upper and lower panels, the left figure represents initial conditions, and the right figures are results after 108 Monte Carlo steps. The open and filled circles indicate liganded and unliganded receptor molecules, respectively. There is no stable cluster formation in the upper snapshot, whereas the clustering is stable in the lower one. Here we use gE = 6kBT, density of liganded receptor: upper plane, 0.001; lower plane, 0.03; and density of unliganded receptor: upper plane, 0.059; lower plane, 0.03.

To obtain more insight into the conditions where receptor clustering can take place, we next analyze the partition function via the mean-field approximation.

    MEAN FIELD APPROXIMATION
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

To proceed, we decouple the quadratic term in the Hamiltonian by introducing an auxiliary Gaussian field and employing the standard Hubbard-Stratonovich/Gaussian transformation (see, e.g., Amit, 1993; Parisi, 1988) (Eq. A3). The benefit of this transformation is to decouple the quadratic terms into linear terms such that we can sum over the ensemble configuration ({ni, sigma i}) at each lattice site i independently. This yields (see Appendix for details)
Z=C<LIM><OP>∫</OP></LIM>𝒟&phgr;e<SUP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)∑<SUB>ij</SUB>J<SUB>ij</SUB>&phgr;<SUB>i</SUB>&phgr;<SUB>j</SUB>+∑<SUB>i</SUB>S<SUB>i</SUB>({&phgr;})</SUP> (8)
with
S<SUB><UP>i</UP></SUB>({&phgr;})=<UP>ln</UP><FENCE>1+2z <UP>cosh</UP><FENCE>&bgr;<FENCE>g<SUB><UP>E</UP></SUB> <LIM><OP>∑</OP><LL><UP>j</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>&phgr;<SUB><UP>j</UP></SUB>+y</FENCE></FENCE></FENCE>
where Dphi  = product idphi i, and C is a normalization constant that does not affect the thermodynamic properties of the partition function. The new field phi  ranges from -infinity to +infinity , y = -L + gL)/2, and z ebeta x = ebeta R+(µL+gL)/2] is related to the receptor chemical potential, which remains to be determined (in terms of the receptor density). The first term in Eq. 8 is related to the interaction energy between nearest-neighbor lattice points, whereas the second term is related to the entropy arising because of the available configurations on an individual lattice site.

In mean-field theory, we try to determine a "homogeneous" saddle point approximation for the partition function. For our system, the negative cooperation (i.e., the AFM nature) suggests that the system might prefer having neighboring sites in oppositely liganded states. Thus, we separate the lattice into two interleaved sublattice systems: all nearest neighbors of a lattice site belong to the alternate sublattice (Fig. 1). We then assign two "uniform" order parameters, phi ±, to each sublattice. After this assumption, the exponent of the Boltzmann factor in the partition function (Eq. 8) now becomes right-arrow (N/2)[beta gEDphi +phi - + S(phi +, phi -)], where N is the number of total lattice sites, S(phi +, phi -) Sigma k=± ln[1 + 2z cosh(beta [gEDphi k + y])], and D is the number of nearest neighbors, which depends on the structure of the lattice. For instance, a square lattice yields D = 4, whereas a honeycomb lattice yields D = 3.

Next, we minimize the free energy by varying phi ±. The variation yields the "saddle point" equation
<FENCE><FENCE><FR><NU>∂</NU><DE>∂&phgr;<SUB><UP>±</UP></SUB></DE></FR> [&bgr;g<SUB><UP>E</UP></SUB>D&phgr;<SUB><UP>+</UP></SUB>&phgr;<SUB><UP>−</UP></SUB>+S(&phgr;<SUB><UP>+</UP></SUB>, &phgr;<SUB><UP>−</UP></SUB>)]</FENCE><SUB><UP>&phgr;<SUB>±</SUB>=</UP><A><AC><UP>&phgr;</UP></AC><AC>˜</AC></A><SUB><UP>±</UP></SUB></SUB></FENCE>=0
Working this out explicitly, we find a self-consistent equation for phi~ ±,
<A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>±</UP></SUB>=<UP>−</UP><FR><NU>2z <UP>sinh</UP>[&bgr;(g<SUB><UP>E</UP></SUB>D&phgr;<SUB><UP>∓</UP></SUB>+y)]</NU><DE>1+2z <UP>cosh</UP>[&bgr;(g<SUB><UP>E</UP></SUB>D&phgr;<SUB><UP>∓</UP></SUB>+y)]</DE></FR> (9)
with the free energy density
f(<A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>+</UP></SUB>, <A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>−</UP></SUB>, z)=<UP>−</UP><FR><NU>g<SUB><UP>E</UP></SUB>D</NU><DE>2</DE></FR> <A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>+</UP></SUB><A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>−</UP></SUB>−<FR><NU>k<SUB><UP>B</UP></SUB>T</NU><DE>2</DE></FR>  (10)

<LIM><OP>∑</OP><LL><UP>k=±</UP></LL></LIM> <UP>ln</UP>[1+2z <UP>cosh</UP>(&bgr;[g<SUB><UP>E</UP></SUB>D<A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>k</UP></SUB>+y])]
Finally, the mean-field receptor density is given by < n>  = -partial f(phi~ +, phi~ -, z)/partial µR. Explicitly, we have
⟨n⟩=<LIM><OP>∑</OP><LL><UP>k=±</UP></LL></LIM> <FR><NU>z <UP>cosh</UP>(&bgr;[g<SUB><UP>E</UP></SUB>D<A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>k</UP></SUB>+y])</NU><DE>1+2z <UP>cosh</UP>(&bgr;[g<SUB><UP>E</UP></SUB>D<A><AC>&phgr;</AC><AC>˜</AC></A><SUB><UP>k</UP></SUB>+y])</DE></FR> (11)
We can therefore determine the receptor chemical potential, x (or equivalently, z), in terms of < n> . Thereafter, we can rewrite the free energy density in terms of < n> , phi~ +, and phi~ -.

    THE ONSET OF CLUSTERING
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

There is no closed-form solution for Eq. 9. To get some analytical information, we define phi~ ± = m ± epsilon and, with
U(w)=<FR><NU>2z <UP>sinh</UP>[&bgr;(g<SUB><UP>E</UP></SUB>Dw+y)]</NU><DE>1+2z <UP>cosh</UP>[&bgr;(g<SUB><UP>E</UP></SUB>Dw+y)]</DE></FR>
we have
m=<UP>−</UP><LIM><OP>∑</OP><LL><UP>k=0</UP></LL><UL><UP>∞</UP></UL></LIM> <FR><NU>&egr;<SUP><UP>2k</UP></SUP></NU><DE>(2k)!</DE></FR> U<SUP>(<UP>2k</UP>)</SUP>(m) (12)

&egr;=<LIM><OP>∑</OP><LL><UP>k=0</UP></LL><UL><UP>∞</UP></UL></LIM> <FR><NU>&egr;<SUP><UP>2k+1</UP></SUP></NU><DE>(2k+1)!</DE></FR> U<SUP>(<UP>2k+1</UP>)</SUP>(m) (13)
where U(k)(w) = dkU(w)/dwk. The basic idea of separating out the epsilon  dependence is that solutions with nonzero values of epsilon  represent phases in which the proximity of neighboring receptors gives rise to alternating ligand binding. For very small receptor densities, there are few neighboring receptors, and hence we expect to find a unique solution of the mean-field equations with epsilon  = 0. In fact, it is clear from Eq. 13 that there is a solution with epsilon  = 0 for all values of the parameters, but at larger densities, there may be other, more stable phases. The goal of our analysis will be to understand the general structure of the phase diagram and then to obtain more quantitative detail by numerical means.

To proceed, let us assume that epsilon  is small and solve Eqns. 12 and 13 to order epsilon 2:
m=m<SUB>0</SUB>−m<SUB>1</SUB>&egr;<SUP>2</SUP> (14)
with
m<SUB>0</SUB>=<UP>−</UP>U<SUP>(0)</SUP>(m<SUB>0</SUB>)

m<SUB>1</SUB>=<FR><NU>U<SUP>(2)</SUP>(m<SUB>0</SUB>)</NU><DE>2![1+U<SUP>(1)</SUP>(m<SUB>0</SUB>)]</DE></FR> (15)

&egr;<SUP>3</SUP>=<FR><NU>{U<SUP>(1)</SUP>(m<SUB>0</SUB>)−1}&egr;</NU><DE>m<SUB>1</SUB>U<SUP>(2)</SUP>(m<SUB>0</SUB>)−<FR><NU>1</NU><DE>6</DE></FR>U<SUP>(3)</SUP>(m<SUB>0</SUB>)</DE></FR> (16)
Using the relationship given above for < n> , it is easy to verify that U(1)(m0) = beta gED[< n>  - m02], U(2)(m0) = -m0(beta gED)2[1 - 3< n>  + 2m02], and
U<SUP>(3)</SUP>(m<SUB>0</SUB>)=([1−3⟨n⟩+6m<SUP>2</SUP><SUB>0</SUB>][⟨n⟩−m<SUP>2</SUP><SUB>0</SUB>]

−3m<SUP>2</SUP><SUB>0</SUB>

[1−⟨n⟩])(&bgr;g<SUB><UP>E</UP></SUB>D)<SUP>3</SUP>
We must consider separately the cases where the denominator of Eq. 16 is positive or negative. Let us first imagine it is positive. Then the existence of a nontrivial solution of Eq. 16 requires that {beta gED[< n>  - m02- 1} > 0. At small < n> this condition will clearly fail, and we will have only the trivial solution. Furthermore, this condition will fail at < n> close to 1 for large enough |y|. We can see this by comparing the equation for m0 with the expression for < n> . Note that if y is large enough such that the hyperbolic functions can be replaced by exponentials, we have |m0| = < n> , and the above expression can be replaced by {beta gED[< n>  - < n> 2- 1}; this is negative for the stated condition. As we cross a line in parameter space such that this factor changes sign to positive, there will be new solutions at nonzero epsilon 2, and the one at epsilon  = 0 becomes a local maximum of the free energy. This emergence of a double-well structure with a continuous growth of the nonzero epsilon 2 solution, indicates that the system exhibits a second-order phase transition.

We must next take into account the possibility that {m1U(2)(m0- <FR><NU><IT>1</IT></NU><DE><IT>6</IT></DE></FR>U(3)(m0)} < 0. Having the denominator cross zero gives rise in our current approximation to a large value of epsilon , which thus invalidates the neglect of higher-order terms. Typically, the higher-order terms will stabilize the system at some finite value of epsilon , which thus appears "spontaneously" as some parameter threshold is crossed. This is a first-order phase transition, or equivalently, a triple-well structure for the free energy. If the local minima (for zero and nonzero epsilon 2) have equally low free energy densities, the system can exist in a mixture of the two phases. As we will see, the two coexisting phases differ in their receptor densities. Finally, the points where both {beta gED[< n>  - m02 - 1] = 0} and {m1U(2)(m0- <FR><NU><IT>1</IT></NU><DE><IT>6</IT></DE></FR>U(3)(m0)} = 0 are "critical end-points" points, because they correspond to places where a second-order transition line ends at a first-order line. A diagram of this behavior, generated by the numerical solution of the mean-field equations, is given in Fig. 3.



View larger version (15K):
[in this window]
[in a new window]
 
FIGURE 3   Numerically computed phase diagram, showing that there is a pair of second-order lines, each of which ends on the first-order transition curve (C = critical point; CE = critical endpoint). The phase to the lower right has epsilon  not equal  0. Here we used gE = 6kBT and D = 3. To show the symmetry, we plot the ligand concentration in logarithm units, normalized with respect to the dissociation constant Kdtnf = 0.59 nM. Here the chemical potential µR is related to receptor density. In Fig. 4, we convert the receptor chemical potential into the molecular density.

For a given ligand concentration, we can find the phase coexistence lines arising because of the first-order phase transition. This is done by finding two solutions (solved with differing values of epsilon 2) of the mean-field equations and then fixing z (as a function of y) by requiring that they have equal free energy,
f(&phgr;<SUP>(<UP>d</UP>)</SUP>, &phgr;<SUP>(<UP>d</UP>)</SUP>, z)=f(&phgr;<SUP>(<UP>c</UP>)</SUP><SUB><UP>−</UP></SUB>, &phgr;<SUP>(<UP>c</UP>)</SUP><SUB><UP>+</UP></SUB>, z) (17)
where phi ±(c) are the order parameters for the dense condensed phase and phi (d) are the (equal) ones for the dilute phase. For the condensed phase, the receptor density is close to unity for reasonable values of the cooperativity parameter beta gED. The workings of this system as far as signaling is concerned are shown in Fig. 4. Assume there is some fixed value of the receptor density. As the ligand concentration is increased, we will cross the phase transition boundary and the receptors will segregate into a condensed phase and a dilute one, corresponding to the two coexisting mean-field solutions. Under our basic hypothesis that signaling is affected by having dense clusters, the response will exhibit a sharp jump at a specific threshold ligand concentration. Similarly, as the ligand concentration becomes too high we cross back to the uniform receptor density state and signaling ceases. That is, we have a ligand concentration "window" for receptor clustering.



View larger version (16K):
[in this window]
[in a new window]
 
FIGURE 4   The phase diagram, now shown as a function of the receptor density < n> related to the ligand concentration. To show the symmetry, we plot the ligand concentration in logarithm units, normalized with respect to the dissociation constant Kdtnf = 0.59 nM. The region inside the solid lines is the coexistence region where states of high and low density coexist. As the ligand concentration is altered so as to cross one of these lines, the receptors will spontaneously cluster and thereby allow signaling to occur. < n> min(d) is the minimal receptor density for clustering. For this set of parameters, clustering will occur even for very small overall receptor density. Here a0 is the length scale for the lattice spacing. For surface receptor molecules such as TNF-R1, we might take a0 approx  1 nm.

As can be seen from the figure, the "clustering" window will cease to exist below some minimal receptor density, as we never enter the phase coexistence region. By symmetry, this minimal density can be found by solving the mean field equations for y = 0 where m = 0. This leads after some algebra to the self-consistent equations
⟨n⟩<SUP>(<UP>d</UP>)</SUP><SUB><UP>min</UP></SUB>=<FR><NU>e<SUP><UP>&bgr;g<SUB>E</SUB>D&egr;<SUP>2</SUP>/2</UP></SUP>−1</NU><DE><UP>cosh</UP>[&bgr;g<SUB><UP>E</UP></SUB>D&egr;]−1</DE></FR> (18)
with
&egr;=⟨n⟩<SUP>(<UP>d</UP>)</SUP><SUB><UP>min</UP></SUB><UP> tanh</UP>[&bgr;g<SUB><UP>E</UP></SUB>D&egr;]
The numerical solution of these equations is presented in Fig. 5. As the cooperativity parameter is increased, the minimum density that will support a clustering window goes rapidly to zero. For the TNF-TNFR1 cluster, it has been speculated that the structure of the cluster is a honeycomb-like lattice (Bazzoni and Beutler, 1995; Naismith et al., 1995, 1996), which implies the number of nearest neighbors D = 3. If we use our rough estimate gE sime  gL/10 sime  6kT, we find that < n> min(d) ~<  10-6/a02. Here a0 is the length scale of the lattice spacing. If we take a0 approx  1 nm, on a cell with surface area 100 µm2, this estimate yields a requirement for less than 102 TNF-R1 molecules distributed on the cell surface. Given that an average number of expressed TNF-R1 on the cell surface is ~2000, we find that the cell operates within the desired part of the phase diagram and hence should exhibit strong sensitivity to the application of TNF. However, we should point out that this estimate is very rough, as we have made a number of simplifying assumptions, and this issue needs to be revisited with a more precise model of the receptor interactions.



View larger version (14K):
[in this window]
[in a new window]
 
FIGURE 5   The variation of < n> min(d) as a function of the association energy gE. We find that < n> min(d) rapidly approaches zero once gED >=  15kBT. If we assign D = 3-4, this energy scale corresponds to a single hydrogen bond. Here a0 is the length scale for the lattice spacing. For a surface receptor molecule such as TNF-R1, we might take a0 approx  1 nm.

    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

We have presented a simple model for signal transduction via receptor clustering, based loosely on the TNF-TNFR1 system. Our basic idea is simple. The interaction between receptors can lead to a first-order phase transition with a discontinuous jump in the receptor density as a function of the receptor chemical potential and/or the ligand concentration. Turning this around, this implies that the receptor system will spontaneously phase separate for a range of ligand concentrations. This fact about the thermodynamic equilibrium state will lead under reasonable kinetic assumptions to the rapid formation of receptor clusters. Assuming that these clusters are necessary for the signal to proceed downstream has the immediate consequence that the system exhibits a strong robust response independent of any details of the intracellular signaling cascade. This might provide a simple solution to the problem faced by biological evolution of how to get a digital response in an analog world.

From a physics perspective, there is nothing very surprising about our phase diagram findings. The idea of a "lattice" Hamiltonian with intrinsic "cooperativity" has been proposed before (Changeux et al., 1967), and on general grounds models of this sort can be expected to have first-order phase transitions. What is new here is the connection of the transition to signaling via the idea of receptor clustering. This connects nicely with increasing evidence that clustering is "universal" among many types of receptor classes.

In our model, we have ignored more-than-two receptor interaction, and relevant internal chemical degrees of freedom (such as the dissociation of SODD in the TNF-R1 system). We do not expect these detailed considerations to change the overall picture, but a more sophisticated model will be needed to make more quantitative estimates of ligand thresholds, cluster structures, and, most interestingly, clustering dynamics. We hope to report on these issues in the future, as well as on the extension of our models to other ligand-receptor systems.

Finally, it would be important to extend our work to later-stage dynamics, as that would allow the consideration of processes such as adaptor protein-mediated receptor internalization, cytoskeleton-assisted cluster stabilization, receptor affinity regulation, receptor cross-talk, and adaptation (Barkai and Leibler, 1997; Hahn et al., 1993; Humphries, 1996; Holsinger et al., 1998; Luo and Lodish, 1997; Stewart et al., 1998; Sundberg and Rubin, 1996; Valitutti et al., 1995; Wyszynski et al., 1997). Other possible extensions might involve the inclusion of spatial fluctuations, the explicit treatment of external perturbations (Shoyab and Todaro, 1981), the local heterogeneity of the microenvironment (Bean et al., 1988; Ward and Hammer, 1992), or fluctuations of ligand concentration; all of these issues have been neglected here.

    APPENDIX: THE GAUSSIAN TRANSFORMATION
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES

The identity
<LIM><OP>∫</OP><LL><UP>−∞</UP></LL><UL><UP>∞</UP></UL></LIM><UP>d</UP>x <UP>exp</UP><FENCE><UP>−</UP><FR><NU>1</NU><DE>2g</DE></FR> x<SUP>2</SUP>+isx</FENCE>=C×<UP>exp</UP><FENCE><UP>−</UP><FR><NU>g</NU><DE>2</DE></FR> s<SUP>2</SUP></FENCE> (A1)
with g > 0, i = <RAD><RCD><IT>−1</IT></RCD></RAD>, and C as some constant, can be generalized to
<LIM><OP>∫</OP><LL><UP>−∞</UP></LL><UL><UP>∞</UP></UL></LIM><LIM><OP>∏</OP><LL><UP>i</UP></LL></LIM> <UP>d</UP>x<SUB><UP>i</UP></SUB><UP>exp</UP><FENCE><UP>−</UP><FR><NU>g</NU><DE>2</DE></FR> <LIM><OP>∑</OP><LL><UP>jk</UP></LL></LIM> x<SUB><UP>j</UP></SUB>J<SUB><UP>jk</UP></SUB>x<SUB><UP>k</UP></SUB>+ig <LIM><OP>∑</OP><LL><UP>j</UP></LL></LIM> s<SUB><UP>j</UP></SUB> <LIM><OP>∑</OP><LL><UP>k</UP></LL></LIM> J<SUB><UP>jk</UP></SUB>x<SUB><UP>k</UP></SUB></FENCE> (A2)

=C×<UP>exp</UP><FENCE><UP>−</UP><FR><NU>g</NU><DE>2</DE></FR> <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> s<SUB><UP>i</UP></SUB>J<SUB><UP>ij</UP></SUB>s<SUB><UP>j</UP></SUB></FENCE>
as long as J is a symmetrical positive definite matrix. Thus, Eq. 8 can be obtained by
Z=<LIM><OP>∑</OP><LL>{<UP>n<SUB>i</SUB>, &sfgr;<SUB>i</SUB></UP>}</LL></LIM> <UP>exp</UP><FENCE>&bgr;<FENCE>x <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> n<SUB><UP>i</UP></SUB>+y <LIM><OP>∑</OP><LL><UP>i</UP></LL></LIM> n<SUB><UP>i</UP></SUB>&sfgr;<SUB><UP>i</UP></SUB>−<FR><NU>g<SUB><UP>E</UP></SUB></NU><DE>2</DE></FR> <LIM><OP>∑</OP><LL><UP>ij</UP></LL></LIM> J<SUB><UP>ij</UP></SUB>n<SUB><UP>i</UP></SUB>&sfgr;<SUB><UP>i</UP></SUB>n<SUB><UP>j</UP></SUB>&sfgr;<SUB><UP>j</UP></SUB></FENCE></FENCE> (A3)

=<LIM><OP>∑</OP><LL>{<UP>n<SUB>i</SUB>, &sfgr;<SUB>i</SUB></UP>}</LL></LIM> C<SUB>&thgr;</SUB><LIM><OP>∫</OP></LIM><LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> &thgr;<SUB><UP>j</UP></SUB>e<SUP><UP>−</UP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)<UP>∑<SUB>jk</SUB>J<SUB>jk</SUB>&thgr;<SUB>j</SUB>&thgr;<SUB>k</SUB>+∑<SUB>j</SUB>&bgr;</UP>[<UP>x+&sfgr;<SUB>j</SUB></UP>(<UP>y+ig<SUB>E</SUB>∑<SUB>k</SUB>J<SUB>jk</SUB>&thgr;<SUB>k</SUB></UP>)]<UP>n</UP><SUB><UP>j</UP></SUB></SUP>

=C<SUB>&phgr;</SUB><LIM><OP>∫</OP></LIM><LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> &phgr;<SUB><UP>j</UP></SUB>e<SUP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)<UP>∑<SUB>jk</SUB>J<SUB>jk</SUB>&phgr;<SUB>j</SUB>&phgr;<SUB>k</SUB></UP></SUP>×<FENCE><LIM><OP>∑</OP><LL>{<UP>n<SUB>j</SUB>, &sfgr;<SUB>j</SUB></UP>}</LL></LIM> e<SUP><UP>∑<SUB>j</SUB>&bgr;</UP>[<UP>x+&sfgr;<SUB>j</SUB></UP>(<UP>y+g<SUB>E</SUB>∑<SUB>k</SUB>J<SUB>jk</SUB>&phgr;<SUB>k</SUB></UP>)]<UP>n<SUB>j</SUB></UP></SUP></FENCE>

=C<SUB>&phgr;</SUB><LIM><OP>∫</OP></LIM><LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> &phgr;<SUB><UP>j</UP></SUB>e<SUP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)<UP>∑<SUB>jk</SUB>J<SUB>jk</SUB>&phgr;<SUB>j</SUB>&phgr;<SUB>k</SUB></UP></SUP>×<LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> <FENCE><LIM><OP>∑</OP><LL><UP>n<SUB>j</SUB>, &sfgr;<SUB>j</SUB></UP></LL></LIM> e<SUP><UP>&bgr;</UP>[<UP>x+&sfgr;<SUB>j</SUB></UP>(<UP>y+g<SUB>E</SUB>∑<SUB>k</SUB>J<SUB>jk</SUB>&phgr;<SUB>k</SUB></UP>)]<UP>n<SUB>j</SUB></UP></SUP></FENCE>

=C<SUB>&phgr;</SUB><LIM><OP>∫</OP></LIM><LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> &phgr;<SUB><UP>j</UP></SUB>e<SUP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)<UP>∑<SUB>jk</SUB>J<SUB>jk</SUB>&phgr;<SUB>j</SUB>&phgr;<SUB>k</SUB></UP></SUP>

×<LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> <FENCE>1+2z <UP>cosh</UP><FENCE>&bgr;<FENCE>g<SUB><UP>E</UP></SUB> <LIM><OP>∑</OP><LL><UP>k</UP></LL></LIM> J<SUB><UP>jk</UP></SUB>&phgr;<SUB><UP>k</UP></SUB>+y</FENCE></FENCE></FENCE>

=C<SUB>&phgr;</SUB><LIM><OP>∫</OP></LIM><LIM><OP>∏</OP><LL><UP>j</UP></LL></LIM> &phgr;<SUB><UP>j</UP></SUB>e<SUP>(<UP>&bgr;g<SUB>E</SUB>/2</UP>)<UP>∑<SUB>jk</SUB>J<SUB>jk</SUB>&phgr;<SUB>j</SUB>&phgr;<SUB>k</SUB>+∑<SUB>j</SUB>S<SUB>j</SUB></UP>({<UP>&phgr;</UP>})</SUP>
where i = <RAD><RCD><IT>−1</IT></RCD></RAD>, phi k = itheta k, Sj({phi }) = ln[1 + 2z cosh(beta [gESigma kJjkphi k + y])], and Ctheta , Cphi are integral constants.

    ACKNOWLEDGMENTS

CG acknowledges the LJIS Interdisciplinary Training Program and the Burroughs Wellcome Fund for fellowship support. He also acknowledges Margaret Cheung for help with the numerical simulation. HL acknowledges the support of the U.S. National Science Foundation under grant DMR98-5735.

    FOOTNOTES

Received for publication 10 November 1998 and in final form 7 June 1999.

Address reprint requests to Dr. Herbert Levine, Department of Physics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0319. Tel.: 858-534-4844; Fax: 858-534-7697; E-mail: levine{at}herbie.ucsd.edu.

    REFERENCES
TOP
ABSTRACT
INTRODUCTION
THE LATTICE HAMILTONIAN
NUMERICAL SIMULATION
MEAN FIELD APPROXIMATION
THE ONSET OF CLUSTERING
DISCUSSION
APPENDIX
REFERENCES