## Abstract

In this paper, we develop stochastic models of receptor binding by a bivalent ligand. A detailed kinetic study allows us to analyse the role of cross-linking in cell activation by receptor oligomerization. We show how oligomer formation could act to buffer intracellular signalling against stochastic fluctuations. In addition, we put forward the hypothesis that formation of long linear oligomers increases the range of ligand concentration to which the cell is responsive, whereas formation of closed oligomers increases ligand concentration specificity. Thus, different physiological functions requiring different degrees of specificity to ligand concentration would favour formation of oligomers with different lengths and geometries. Furthermore, provided that ligand concentration specificity is taken as a design principle, our model enables us to estimate parameters, such as the minimum proportion of receptors, that must engage in oligomer formation in order to trigger a cellular response.

## 1. Introduction

Cellular responses are frequently triggered by detection of extracellular signalling molecules (ligands) by specialized cell-surface proteins called receptors. Upon receptor binding, an intracellular signalling cascade takes place, eliciting a cellular response.

There are several types of membrane-bound receptor. Broadly speaking, they can be divided in those that depend on oligomer formation for triggering a cellular response and those which do not (e.g. G protein-linked receptors; Helmreich 2001). Our focus in this study will be on the former class of receptors.

Receptor tyrosine kinases (RTKs) bind to hormones, growth factors and cytokines. B cell and T receptors, which are antigen(epitope)-specific and play a fundamental role in triggering the immune response, are closely related to RTKs (Tamir & Cambier 1998). The basis of RTK binding and activation has been a major problem in cellular signalling. Unlike G protein-linked receptors, which are heptamers that cross the cell membrane several times, RTKs are single polypeptide chains that cross the membrane only once. The problem is how a receptor with such structure is activated upon ligand binding (Helmreich 2001). The answer is by oligomerization also commonly referred to as *cross-linking* in the biophysical literature. Hormones, growth factors, cytokines and antigens are usually multivalent ligands, i.e. they have more than one binding site and can engage several receptors (Lauffenburger & Linderman 1993). Following receptor oligomerization, the RTK domains within their cytoplasmic tails cross-activate each other. Upon cross-activation, the phosphate groups attached to selected tyroresidues provide high-affinity docking sites for several tyrosine kinases carrying the SH2 domain, most notably the *src* family of tyrosine kinases (Alberts *et al*. 2002). This process constitutes the early steps in cellular response to an extracellular signal.

The literature on models of multivalent receptor/multivalent ligand binding is extensive and we do not intend to produce a complete review only noting a few works focusing on receptor oligomerization. For an extensive review of early work on receptors, in general, we refer the reader to the monograph by Lauffenburger & Linderman (1993). Goldstein *et al*. (2004) review most recent developments in the field, specially related to immune system receptors.

Early theoretical work on formation of bivalent receptor aggregates by bivalent ligands was conducted by Perelson & DeLisi (1980). This work was later on extended by Posner *et al*. (1995). Perelson & DeLisi (1980) study the formation of linear aggregates of arbitrary length. Introducing the so-called equivalent site approximation (i.e. all binding sites have the same property regardless of the size of the aggregate to which they belong), they reduce an infinite number of ordinary differential equations (ODEs) describing the formation of aggregates of arbitrary size to two equations (for the free and bound ligand concentration) and a conservation law. Moreover, Perelson & DeLisi (1980) show that the concentration of chains of any size can be expressed in terms of the three ligand state concentrations (free, one bound site and two bound sites).

However, these results have been proved by Posner *et al*. (1995) to hold only if linear aggregates are allowed. If rings can form, the equivalent site approximation fails and such general results cannot be obtained. Posner *et al*. (1995) show that if rings are allowed to form only up to a maximum size, *k*, the system of infinite ODEs can be reduced to a system of equations. Posner *et al*. (1995) also show that in this case it is not possible to find an expression for the aggregate size distribution.

McKeithan (1995) introduces the concept of kinetic proof-reading in the context of T cell receptor binding to explain how small affinity differences between foreign and self-antigens could elicit T cell activation. Hlavacek *et al*. (2001) gives a more recent exploration of the relationship between cross-linking effects and kinetic proof-reading.

Sulzer *et al*. (1996) reconsider the role of cross-linking in the context of the cellular response to a collection of bivalent ligands with different affinities for the receptor as, for example, would be the case in a polyclonal anti-receptor serum. To analyse this system, they defined a *binding field* and a *cross-linking field*. In terms of these quantities, they study the properties of the cross-linking curve and find conditions for clonal expansion under stimulation with a mixture of different ligands.

Further refinement of the models of multivalent receptor/multivalent ligand binding was achieved by Hlavacek *et al*. (1999) when considering steric hindrance by introducing a *steric hindrance factor*, which accounted for the fraction of unbound ligand sites accessible to receptors and therefore available for binding. A detailed analysis of the relationship between receptor clustering and signalling in the particular case of cell-surface IgE-FcϵRI has been recently carried out by Posner *et al*. (2004).

MacGabhann & Popel (2005) studied the system vascular endothelial growth factor (VEGF)/VEGF receptor 2/neuropilin-1. VEGF receptor 2 (VEGFR2) and neuropilin-1 (NRP1) are both found on the surface of endothelial cells. They do not interact directly but can be cross-linked by a VEGF isoform which has binding sites for both VEGFR2 and NRP1.

In this paper, we show how oligomer formation could act to buffer intracellular signalling against stochastic fluctuations when the number of receptors is smaller than in normal conditions. In addition, we put forward the hypothesis that the use of dimers, trimers and rings of cross-linked receptors could serve different physiological purposes requiring the cells to respond to different ligand concentration ranges. Furthermore, if we accept ligand concentration specificity as a valid design principle for receptor function, we show how our model, given that we know other parameter values (mainly affinities and number of receptors), could be used to estimate quantities, such as the minimal proportion of receptors, engaged in oligomer formation in order to trigger a cellular response.

This article is organized as follows. In §2, we formulate our model for cell receptor cross-linking by bivalent ligands and present simulation results. In §3, we summarize and discuss the main results of this analysis and present our conclusions.

## 2. Models of receptor cross-linking by bivalent ligands

In §1, we have reviewed some mathematical models of immune-receptors and the different issues they address. Many of them are formulated in terms of nonlinear ODEs using the law of mass action (LMA; Murray 2001). This approximation is known to be valid for large, homogeneous systems. When either of these two assumptions fail to hold, new phenomena arise that the LMA fails to predict or explain (Durrett & Levin 1994; Togashi & Kaneko 2001; Louzon *et al*. 2003; Schnell & Turner 2004). In this section, we investigate a spatially homogeneous distribution of receptors on the cell surface. However, in our study of the behaviour of oligomerized receptors, we allow the total number of receptors to be small.

The purpose of this is to determine what the impact would be of reducing the expression of receptors on signalling. This might be a relevant issue, for example, to assess the influence of reducing mIg on the treatment efficacy in lymphoma B cells under certain types of immunotherapy (Vitetta & Uhr 1994; Page & Uhr 2005), in which an anti-idiotypic antibody is used either passively or stimulated by a vaccine.

In particular, we formulate a stochastic model of receptor binding and oligomerization by bivalent ligands. Stochastic models of reacting systems can be formulated as Markov processes in terms of a Master equation (Gardiner 1983). For a general chemical system, the corresponding Master equation reads:(2.1)where is a vector whose components are the number of the different molecules involved in the reaction, the components of the vector are the corresponding increments in the number of molecules after reaction *i* has occurred, are the reaction probabilities per unit time and is the probability density of finding the system in a given configuration at time *t* (Kubo *et al*. 1973). The formulation in terms of a Master equation is mathematically very elegant, but its analytical treatment for multi-dimensional systems becomes extremely difficult1 (Gillespie 1977). In chemical systems, the reaction probabilities per unit time satisfy (Kubo *et al*. 1973; Gillespie 1977)(2.2)where *N* is the total number of molecules in the system and . This allows us to rewrite the Master equation (2.1) as(2.3)If , then a Wentzel–Kramer–Brillouin (WKB) ansatz is used to obtain the solution of equation (2.3) (see Kubo *et al*. 1973). However, since we are interested in studying this system for small values of *N*, such approximate solutions are not too interesting. Instead, we perform numerical simulations of the reaction mechanism shown in figure 1 using Gillespie's exact stochastic simulation algorithm (ESSA; see appendix A and Gillespie (1977)), which is equivalent to the Master equation.

We start by formulating a model in which only dimers of receptors can be formed (as would be the case for many growth factors (Alberts *et al*. 2002), in which bivalent ligand binds monovalent receptors). Hereafter this model will be referred to as Model 1. Afterwards we will formulate models allowing the formation of larger oligomers in the form of linear (to be referred to as Model 2) or closed chains (to be referred to as Model 3).

Let *U* be the number of unbound receptors and . Similarly, we define *B* and *b*, and *X* and *x* as the number and proportion of bound receptors and cross-linked receptors, respectively. Hence, the state vector will be given by and, correspondingly, .

The only information we need to specify our Markov model is , According to figure 1*a*, there are four elementary reaction steps (i.e. receptor binding and unbinding and dimer formation and dissociation) characterized by a vector , where , , and are the increments in the number of unbound receptors, bound receptors and dimers, respectively, corresponding to the elementary transition *i*. Note that, due to conservation of the total number of receptors the components of must satisfy , . The reaction probabilities per unit time (transition rates) are given in table 1. In appendix A, we give full details of how they have been obtained. The transition rate corresponding to reaction needs further clarification. The transition rate for this reaction corresponds to the formation of a dimer. It is obtained as the product of two factors, namely, the rate of binding between an unbound receptor and a ligand–receptor heterodimer and the probability of finding another receptor within a distance (see figure 1 for its physical meaning) of the ligand–receptor heterodimer. The later is given by , with being the surface density of receptors on the cell surface and *R*, the average radius of a B cell.

Our first simulation results, shown in figures 2 and 3, show the effect of reducing the total number of surface receptors. We have done simulations with , , and 50. Figure 2 shows particular realizations of the time evolution of the system. Figure 3 shows the stationary probability density for different values of *N*, averaged over a number of realizations (see caption of figure 3 for details).

The results shown in figure 2 show how decreasing the number of receptors has a strong effect on the behaviour of the system. For (figure 2*a*), fluctuations are hardly noticeable. In fact, the stationary probability density is strongly peaked around the average stationary state, exhibiting a Dirac's delta-like behaviour: the variance around the average stationary state is very small. As the total number of receptors decreases, we obtain a different picture: as expected (Kubo *et al*. 1973), the effects of fluctuations become more important. For (figure 2*a*), the behaviour of the system is Gaussian-like. The effects of fluctuations become more dramatic if the number of receptors is reduced further. As shown in figures 2*c* and 3*b* and 2*d* and 3*c*, decreasing further the number of receptors yields a transition between the Gaussian-like behaviour exhibited for larger receptor numbers and a multi-stable behaviour.2 This transition is a type of *noise-induced transition* (Hormsthemke & Lefever 1984), as the corresponding (LMA) deterministic system does not exhibit such a transition. Its nature is similar to the *discreteness-induced transition* observed by Togashi & Kaneko (2001) in a model for small autocatalytic systems.

These simulations show that reducing the number of receptors yields a *diluted* system, in which the predictions of LMA fail to capture some of the essential characteristics of the system. Hence if we want to investigate the response of a B cell to a bivalent antigen as the number of B-cell receptors (BCRs) decreases, the stochastic formulation seems to be more appropriate.

### 2.1 Dynamics of the binding and unbinding process: importance of kinetic proof-reading

One of the aims of the present study is to assess, in a more general setting, whether kinetic proof-reading (McKeithan 1995) explains the importance of receptor cross-linking for building up an appropriate cellular response. In fact, it has been already suggested that in other immune-receptors there may be ways to escape kinetic proof-reading (Hlavacek *et al*. 2001).

Kinetic proof-reading is wielded as a reason why cross-linking is necessary for an efficient response of B cells (and in fact many other receptor binding-mediated cellular responses): the receptor needs to be bound to its ligand for long enough for the activation process of the proteins (usually kinases) and/or second messengers mediating the response to be completed and the signalling pathways initiated (McKeithan 1995).

Let us assume that is the threshold number of cross-links that need to be formed in order to produce signalling within the cell (Sulzer *et al*. 1996). If we assume that ligation of a single receptor can also produce signalling, the threshold number of bound receptors needed for signalling can be assumed to be .

In order to estimate whether oligomers are, on average, more stable than bound ligands, we proceed as follows. Let be the rate of dissociation of a ‘typical’ Src protein tyrosine kinase (PTK) from a bound receptor. Its inverse is a measure of the average time a Src PTK remains bound to a bound/cross-linked receptor. Since RTKs and immune receptors depend on PTKs to relay the signal (Tamir & Cambier 1998; Alberts *et al*. 2002; Goldsby *et al*. 2003), they cannot trigger any signalling pathways, unless they remain bound/cross-linked for long enough (longer than ). Hence, we use the quantities and , i.e. the (steady-state, which in this context means that the moments of the probability density do not change in time) probabilities that at time *X* remains and at time *B* remains , respectively, with () defined as the time of the first occurrence of the event (). In our case, we take . Hlavacek *et al*. (2001) have estimated that, for FcϵRI receptors, is 2.5 times smaller than . We will use . See appendix A for full details of how and are computed.

The results shown in figure 4 demonstrate that dimers are not intrinsically more stable than monomers.3 They simply appear to be stable at different ranges of ligand concentration. Hence, kinetic proof-reading does not seem to be an important factor in explaining oligomer formation in the kind of system under examination here. A reason for this might be that, in T cell receptors, the differences in affinity between self- and foreign antigens are typically very small, hence the necessity for kinetic proof-reading (McKeithan 1995). The receptors we are looking at have high affinities for their ligands and therefore kinetic proof-reading might not be an issue.

### 2.2 Models with formation of longer oligomers

In this section, we extend our previous model to include the formation of longer oligomers, in particular trimers (Model 2). The number of trimers is denoted by *Y* and their concentration4 . The number of dimers bound to two molecules of ligand is denoted by and their concentration (see figure 1). Hence, in this case the state vector will be given by and, correspondingly, . The reaction probabilities per unit time for this model are given in table 2 and have been obtained following the procedure detailed in appendix A.

It has also been observed that trimers can be present as linear (open) chains or can form rings (closed chains) (Lauffenburger & Linderman 1993). In this case, we will consider only formation of chains of up to three cross-linked receptors that then can form rings (Model 3). The state vector will be given by , where *Y* is the number of linear trimers and *R* is the number of rings-like trimers. Correspondingly, , where and . The reaction probabilities per unit time for the model with linear trimers forming rings are given in table 3 (see appendix A for details how they have been obtained).

These simulation results seem to indicate one possible role of cross-linking. Comparing the results shown in figures 2–6 and 10, we can see that forming long oligomers (longer than dimers) reduces the effects of the fluctuations. Hence, it might be argued that cross-linking acts like a noise-suppressor mechanism. It is also interesting to observe that the size of the fluctuations is smaller in Model 3 than in Model 2, even when rings are formed by three cross-linked receptors. This result seems to indicate that not only the number of receptors within a cluster, but also the topology of the cluster appears to be important for noise-suppression (figures 7–11).

On the other hand, from figures 8 and 12, we can reach the same conclusions regarding kinetic proof-reading as in the previous case.

From the results shown in figure 14, we observe that the width of the bell-shaped curves corresponding is bigger for Model 2 (figure 8) than for Model 1. Interestingly, the width of for Model 3 is smaller than the width for for Model 2. Moreover, comparing figure 14*a–c* shows that for Model 1 is biased towards smaller ligand concentrations than for Model 3. These observations might shed some light on the physiological roles of receptors forming oligomers of different lengths and geometries (linear or closed). According to our results, monovalent receptors binding a bivalent ligand (only dimer formation is allowed in such under these conditions) might be used for response to signal which are present in a specific, rather low, range of concentration (see figure 14, dimer, ). This is the case, for example, for growth factor receptors, in which our prediction seems to be supported by the fact that growth factors are usually present only at very low concentrations ( according to Alberts *et al*. 2002). In this scenario, receptors forming longer linear chains would be used to detect signals present in broader range of concentrations, whereas ring formation would be used as a mechanism to recover ligand concentration specificity in cell response.

Another interesting feature of the quantity is observed when looking at figures 5, 9 and 13. As the total number of receptors increases the width of the curve as a function of log(*AL*) changes. In the case of Model 1 it increases, indicating that the receptor is potentially responsive to a wider range of ligand concentration. On the contrary, for Models 2 and 3 the width of as a function of log(*AL*) decreases as *N* grows.

When the total number of receptors is increased in Model 1, the range at which the cell is responsive shifts to higher values of the parameter *AL* (higher ligand concentrations).

### 2.3 Ligand concentration specificity as a design principle

The models described in this paper might be useful to estimate parameter values. In particular, , i.e. the threshold proportion of receptors that need to be engaged in oligomer formation to produce a cellular response, and , the *cross-linking affinity* may be estimated. A further assumption is needed, namely, that ligand concentration specificity (i.e. the cell responds only to a narrow range of values of the ligand concentration) is a ‘receptor design principle’.

Making this assumption, we can observe from figure 14 that for the models with trimer formation and ring formation, the value predicted by our model and the ligand concentration specificity design principle yields , which means that atleast 90% of the receptors are engaged in oligomer formation.

Using this design principle, we might be able to say something about the value of for Models 2 and 3. Concerning Model 2, the results in figure 15*a* show that values of yield a reduced range of ligand concentration in which the cell is responsive to the signal, as compared to the results shown in figure 14*b* obtained with . As for Model 3, specificity is achieved for roughly the same range of values of .

## 3. Conclusions and discussion

We have presented and analysed three models of ligand binding-induced oligomerization of cell surface receptors, accounting for the formation of oligomers of different lengths (dimers and trimers) and different topologies (linear trimers or rings). This approach has allowed us to dissect the effect on cellular signalling of the formation of each one of these three different structures under stimulation by bivalent ligand.

The stochastic formulation of the model has been used to study the behaviour of these models as a function of the number of receptors. Not surprisingly, small numbers of receptors induce fluctuations. However, the intensity of these fluctuations depends on the size and topology of the oligomers formed by cross-linking. In fact, cross-linking appears to buffer the system against noise, as fluctuations are smaller in Model 3, in which rings are allowed to form, than in Model 2 in which only linear trimer formation is considered. Noise intensity in Model 2, in turn, is smaller than in Model 1, where only dimer formation is allowed.

The curves showing as a function of ligand concentration can be interpreted as describing cell responses to extracellular stimulation. Cross-linking produces cellular response only in a bound interval of ligand concentration more or less centred around the receptor affinity. This property has been previously observed in a number of experimental (Banu & Watanabe 1999; Yamazaki *et al*. 2002) and theoretical studies (Lauffenburger & Linderman 1993; Sulzer *et al*. 1996). However, its possible therapeutic relevance has been overlooked in some pathological situations, as for example in cancer.

Two examples in the context of anti-tumour therapy, in which ligand/receptor binding dynamics appears to play an important role, are anti-angiogenic therapy in solid tumours and immunotherapy in B-cell lymphomas. Regarding the former, many of the current protocols under investigation involve drugs which would neutralize the several types of tumour-angiogenic factors, most notably VEGF (table 4). In other words, the aim is to reduce the concentration of ‘active’ growth factor. However, according to the results obtained with Model 1 (which would be applicable to most growth factor/growth factor receptor systems), an equally effective strategy would be to increase the levels of VEGF, in particular, if this is combined with an increase in the minimum number of receptors engaged in oligomer formation to trigger signalling, (see figure 14*a*).

Lymphoma B cells carry on their surface an antibody with a unique idiotype. An antibody to this idiotype can discriminate lymphoma B cells from normal B cells and can be administered passively or produced by the immune system in response to a vaccine. In either case, the anti-idiotype binds the antibody on the surface of the lymphoma B cells, triggering a signalling cascade, which may decide the fate (i.e. apoptosis or quiescence) of these cells. Models 2 and 3 show how kinetic properties of the anti-idiotype/antibody binding process could help to increase or reduce the range of responsiveness of the lymphoma B cells to therapeutic antibody.

The analysis of the feasibility of such approaches to cancer treatment as well as further investigation of the therapeutic implications of the biophysical properties of ligand/receptor binding and elaboration of these arguments will be the subject of future research.

We have found that oligomers of different lengths and topologies yield different intervals of responsiveness to the external stimulus. Formation of longer (linear) oligomers appears to increase the range of ligand concentration to which the cell is responsive. However, inclusion of ring formation in the models produces a reduction of this responsiveness interval, i.e. formation of closed oligomers makes the cell response more specific to ligand concentration. We could argue that mechanisms in which the cell is responding to a wide range of ligand concentrations will depend upon receptors forming long linear oligomers, whereas response to stimuli in a more specific range of ligand concentrations will depend upon receptors forming shorter oligomers or rings.

Furthermore, under the assumption that ligand concentration specificity is a design principle for receptor function, our model allows us to estimate the values of parameter that are not available in the literature, as for example and , i.e. the threshold proportion of oligomerized receptors to produce signal and the cross-linking affinity, respectively.

## Acknowledgments

The authors would like to thank Robin Callard for useful comments. T.A. and K.M.P. thank the EPSRC for financial support under grant GR/509067 and K.M.P. thanks the Joint Research Councils (EPSRC, BBSRC, MRC) for support under grant GR/R47455.

## Footnotes

↵Other approaches in terms of Markov chains may be formulated (see ch. 4 of Lauffenburger & Linderman (1993) and references therein). Nevertheless, as with the Master equation formulation, this approach becomes very difficult to handle analytically for multi-dimensional systems.

↵This means that the probability distribution goes from exhibiting a single maximum (Gaussian-like behaviour) to presenting several maxima, giving rise to multiple stable states and (random) transitions between them.

↵Here, we understand by ‘stability’ that either the number of cross-linked receptors or the number of bound receptors remain above a certain thershold long enough to trigger receptor binding-induced signalling.

↵Hereafter, when the term ‘concentration’ refers to the variables of our models it is to be understood as the corresponding proportion with respect to the total number of receptors. When referred to the ligand, its meaning must be understood to be the usual one.

↵Other algorithms, some of them more computationally efficient, have been recently reviewed by Burrage

*et al*. (2004).↵All the second-order reactions considered in our model are of this type.

↵The LMA is a mean field approximation.

- Received December 20, 2005.
- Accepted January 26, 2006.

- © 2006 The Royal Society