Self-consistent theory of reversible ligand binding to a spherical cell.

Ligand binding to receptors is the initial event in many signaling processes, and a quantitative understanding of this interaction is important for modeling cell behavior. In this paper, we study the kinetics of reversible ligand binding to receptors on a spherical cell surface using a self-consistent stochastic theory. Binding, dissociation, diffusion and rebinding of ligands are incorporated into the theory in a systematic manner. We derive explicitly the time evolution of the ligand-bound receptor fraction p(t) in various regimes. Contrary to the commonly accepted view, we find that the well-known Berg-Purcell scaling for the association rate is modified as a function of time. Specifically, the effective on-rate changes non-monotonically as a function of time and equals the intrinsic rate at very early as well as late times, while being approximately equal to the Berg-Purcell value at intermediate times. The effective dissociation rate, as it appears in the binding curve or measured in a dissociation experiment, is strongly modified by rebinding events and assumes the Berg-Purcell value except at very late times, where the decay is algebraic and not exponential. In equilibrium, the ligand concentration everywhere in the solution is the same and equals its spatial mean, thus ensuring that there is no depletion in the vicinity of the cell. Implications of our results for binding experiments and numerical simulations of ligand-receptor systems are also discussed.


Introduction
A large number of cellular processes are triggered by ligand binding to cell surface receptors. In contrast to reactions taking place within a test tube, ligand-receptor cell binding involves reactants that are not `well-mixed', primarily because the receptors are confined to the cell surface. In addition, the ligand concentration close to the surface is subject to fluctuations arising from association and dissociation processes. The kinetics of ligand-receptor binding could, therefore, deviate substantially from the traditional ordinary differential equation (ODE) based descriptions, which are appropriate for wellmixed reactants whose densities remain constant in space [1].
A related issue is the history-dependence of the binding process arising from the reattachment of dissociated ligands to the cell surface at later times. For this reason, any effective equation describing the binding kinetics has to be necessarily non-local in time, an interesting feature which does not seem to have been widely appreciated in the literature. The traditional approach, starting with the classic paper on bacterial chemosensing by Berg and Purcell [2] and later extended by other authors [3,4,5] has been to find the effective rate constants by solving for the steady state ligand flux to the cell surface. Berg and Purcell showed that, for a spherical cell of radius a , with N receptors on its surface, each binding the ligand with an intrinsic on-rate + k , the net ligand flux (without dissociation) at the surface is (1) is the `effective' on-rate for the process and D is the diffusion coefficient for the ligand. DeLisi and Wiegel [3], and later, Shoup and Szabo [4] showed that the off-rate also undergoes a similar scaling because of rebinding of ligands to the surface, following dissociation at early times. Shoup and Szabo called the factor + + + = Nk aD Nk π γ 4 (2) the absorption probability of a ligand upon contact (and, consequently γ − 1 is the probability of non-absorption). It then follows that the effective rate constants may be expressed as γ πaD , which gives a nice heuristic interpretation of these results.
The scaling relations for the binding and dissociation rates mentioned above have been found to agree well with experimental results [6,7], are widely used in biochemical modeling, and have also been extended to cell surface reactions. In addition, Berg and Purcell's well-known result for the relative error in concentration measurement by a chemotactic cell (like the bacterium E. Coli) follows directly from Eq.1. It is, therefore, of interest to re-examine the validity of the quasi-steady state assumptions (where the ligand concentration is assumed to reach a stationary profile much faster than the receptors) under which these relations were derived and look for possible modifications in a more systematic theory. In addition, we note that in the previous theoretical treatments, association and dissociation processes are treated separately, and, in particular, effects of ligand rebinding on the association kinetics are neglected altogether [8,9]. An attempt towards a unified treatment of binding and dissociation was made by Goldstein and Dembo [10]. However, this technique (method of weighted residuals) is essentially perturbative in nature and does not always predict the correct kinetic behavior. For example, in the pure dissociation problem, the asymptotic decay is algebraic, but MWR predicts exponential decay [10].
In this paper, we generalize a recently introduced self-consistent stochastic formalism [11,12] originally developed to study ligand rebinding to uniformly distributed receptors on a planar surface [12], and later extended to include receptor clusters [11], to study association and rebinding of ligands to receptors on a spherical cell in solution. Being based on the statistics of individual ligand trajectories, this formalism eliminates the need to impose ad hoc boundary conditions for the variation of ligand concentration close to the cell surface. We derive the ligand-bound receptor fraction explicitly as a function of time, and thus derive the effective rate constants. We observe deviations from Berg-Purcell scaling close to equilibrium, and show that the association rate equals the intrinsic rate (and the ligand concentration at the surface equals its spatial mean value on solution) at equilibrium. We discuss experimental implications of our findings, and also test the impact of our predictions in a simple computational model based on the epidermal growth factor and receptor (EGF-EGFR) system.

General formalism
Let us consider a spherical cell of radius a in three dimensions, with 0 2 4 R a N π = receptors distributed randomly on its surface, 0 R being the surface density. The cell is exposed to a ligand in the surrounding medium, whose concentration at 0 = t and at infinite distance from the cell is 0 ρ . Let us denote by the fraction of ligand-bound receptors at time t . The basic kinetic equation for where ) , ( t r ρ is the radially symmetric ligand concentration at distance from the center of the cell, at time t . − k and + k are the dissociation and association rates of the ligand to the receptor. For the rest of this section, we will assume that internalization of ligandreceptor complexes by the cell can be neglected in comparison with dissociation.
is the density contribution `from the bulk', i.e., from the ligands which have not been bound to a receptor until time t (although they may have been `reflected' earlier by boundreceptors or the non-receptor part of the surface). The ligands that make up ) , ( t a r ρ are those which were receptor-bound at an earlier time t < τ and released during the . Let us now define a Green's function for ligand diffusion around the semi-absorbing sphere centered at 0 is the probability density (dimension 1/length 3 ) that a ligand that starts at a radial distance r from the origin at 0 = t is located close to the surface at time t , with possibly multiple visits to the semi-absorbing surface in between (but no absorption by a receptor). Using this Green's function, the above two densities may now be expressed as The first relation is obvious.
. Note that we have assumed complete spherical symmetry of ligand diffusion, which is justified for a uniform distribution of receptors as is the case here.
The separation of the surface ligand density into bulk and rebinding-originated parts is the first new feature of our formalism. The second feature is the expression of these densities in terms of Green's function for individual ligand diffusion trajectories (in the presence of a semi-absorbing surface). The Green's function is computed in the next subsection.

Calculation of
) , ( t r G N The formalism in this section is very similar to what was presented in our earlier paper [11]. Let us first define the first passage probability density ) , ( t r q , which is the probability (per unit volume) that a ligand at radial distance r at t=0 arrives close to the surface for the first time at time t, which is related to ) Here, λ is a microscopic length scale and δ is the time interval for which a ligand will reside inside a volume element 3 λ , before moving away by diffusion. The first term gives the probability that the ligand arrives at the surface for the first time at t, the second term gives the sum of all other events. The factor λ δ η gives the probability of non-binding upon contact with the surface (In general, close to equilibrium, we need to correct for the number of occupied receptors, which can be done approximately by replacing 0 R by ) 1 ( 0 s p R − in η . See also the note following Eq.9b).
In order to calculate ) , ( t r q , we next consider the case 0 0 = R so that the surface is completely reflecting as far as ligand diffusion is concerned. In this case, 1 = η and hence Eq.5 reduces to We will now re-express Eq.5 and 6 in terms of Laplace-transforms ) , ( s r q and ) , ( s r G N , and eliminate ) , ( s r q between the equations, which produces the elegant relation After substituting Eq.8 into Eq.7 and using the basic equations Eq.4a and 4b, we arrive at the following expressions for the Laplace transforms of the densities: For the bulk density, we find where the last expression illustrates the approach to the steady state value at large times, . The density contribution from receptor-released ligands (after using the theorem for convolutions in Eq.4b) is given by For kinetics close to the steady state, the factor N in Eq.2 and the denominators of Eq.7, Eq.9a-c should be corrected for the number of ligand-bound receptors, which is done by replacing N by ) 1 ( s p N − in these places (See Eq.18 and Eq.19 later).

Results
The kinetics in various time regimes may now be worked out by inverting Eq.9a and 9c and substituting in Eq.3. We now study three important cases of interest.
, which is to be expected. In this limit, ) (t r ρ is negligible since release and rebinding events are significant only after a considerable fraction of receptors are bound to ligands. The association curve in this regime is, therefore,

Intermediate regime (
In this case, we consider the regime beyond the (microscopic) time scale D a / 2 , but the receptor occupancy by ligands is still small, i.e., 1 ) ( << t p . In this case, using the i.e., the bulkcontributed density is now only a fraction of the (constant) density at infinity. Eq.9b can also be inverted trivially in this regime to give Let us now substitute for For consistency, we neglect the quadratic term in p(t), and to linear order we have We note from Eq.13 that in this regime, time is `inflated' by a factor 1 ) 1 ( − − γ , which may also be effectively absorbed into the reaction rates to define effective forward and backward rates Note that these effective rates are the same as predicted by the quasi-steady state arguments discussed in the introduction [2,4].

Late time regime and steady state
: We now proceed to analyze the kinetics in the asymptotic regime, when the system is close to equilibrium. Let us first assume that at very late times, the quantities ) ( ), close to the steady state. In terms of Laplace-transformed variables, the equivalent relations are We now substitute these expressions in Eq.9a and Eq.9cand the following relations follow easily. The steady state bound fraction is d s where the equilibrium surface concentration is Note that the scaling factor s γ − 1 appearing in Eq.18 is slightly different from the factor γ − 1 in Eq.9a, and is defined as As mentioned earlier, this (approximate) modification reflects the reduction in the number of available sites for ligand association as binding progresses. Upon solving Eq.17 and Eq.18 together, it is easily seen that 0 ρ γ ρ s s r = and hence i.e., the ligand density close to the surface at equilibrium is the same as the bulk density! (This is, in fact, true everywhere in the solution, see Appendix C). Note that this result is crucially dependent on the inclusion of the rebinding density r ρ and the correction factor in Eq.18 and Eq.19. From Eq.17 it also follows then that ) the equilibrium bound fraction is insensitive to kinetic modifications of rate constants (ie., extrapolating Eq.10 and Eq.13 to ∞ → t gives the same equilibrium value).
Having found the steady state in a precise manner, we now proceed to calculate the leading kinetic terms close to the steady state. After substituting Eq.15 into Eq.3, we find the following equation for ) (t p δ , up to linear terms in the deviations.
Eq.20 is Laplace-transformed as (21) Note that we have used a formal initial condition s p p − = ) 0 ( δ however we recognize that our present analysis is valid only close to the steady state. In order to find It then follows that In the next step, we expand the functions 1 f and 2 In the final step of the procedure, the expansions Eq.24a and Eq.24b are used in Eq.23, which is then substituted into Eq.21. After a few elementary simplifications, we arrive at the final result: We see that this regime is characterized by a time-scale Eq.25 may be explicitly inverted for two limiting cases.
, and At very late times, therefore, the binding curve approaches its equilibrium value algebraically. The power-law approach to the final steady state is in agreement with a general result in the theory of bimolecular reactions, which predicts that, in d dimensions, equilibrium is approached asymptotically as been assumed here. In any realistic situation, however, the number of ligands is finite, and therefore, the absorption by receptors will result in a reduction in their number and hence, the mean concentration. Nevertheless, our main conclusion that there is no depletion region for the ligand concentration close to the cell, should still be valid.
In Fig.1 is also shown for the same parameter values, which clearly displays the non-monotonic variation in time as argued above.

Time-dependent association rate
The time-dependence of the surface ligand concentration may also be framed in terms of a time-dependent association rate, which may be defined as The advantage of this approach, especially for numerical simulations, is that all the effects of concentration fluctuations due to ligand release and rebinding are taken care of in this new effective association rate, and no separate treatment of rebinding (eg. by using a reduced dissociation rate as in Eq.13) is needed. For the dissociation rate, we may therefore simply use the intrinsic value. ) (t r ρ is given by Eq.11. Upon combining the two, and making use of Eq.17, we find that The late-time form is a little different, and involves the binding-corrected factor s γ (Eq.19).
A compact expression that smoothly interpolates between Eq.30 and Eq.31 in the appropriate time-regimes may be obtained by defining the factor )) ( 1 ( 4 using which we find an approximate expression for the effective on-rate for all times. Eq.33 provides a compact, although not completely rigorous (since the introduction of the factor 1-p(t) in Eq.32 is of an ad-hoc nature, and this additional time-dependence has not been included in the Laplace transform method used earlier), expression which smoothly interpolates between the correct limiting forms given by Eq.30 and Eq.31. This form may be suitable for use in computational studies of ligand-receptor interactions, where rebinding and consequent fluctuations in ligand concentration close to the surface are deemed important and may be captured at least approximately. A specific example is discussed in Section 5.

Experimental implications
In experimental studies of ligand-receptor binding, the quantity that is typically measured is the time evolution of ) (t p , for which the effects of ligand release and rebinding on the binding kinetics are more subtle and less direct (as opposed to the ligand concentration itself). Our analysis predicts that as long as the bound fraction is sufficiently small, the binding curve may be very closely approximated to an exponential, with both the on and off-rates scaled by the same factor γ − 1 (Eq.13 and Eq.14). However, closer to the steady state, the time scale of the exponential undergoes a change and switches to s T , in which only the off-rate is scaled (by s γ − 1 ) (Eq.26). Although, there is a significant change in the surface ligand concentration, the qualitative nature of the binding curve is likely to be only subtly altered since both time scales are involved in determining the kinetics of evolution. In particular, the cross-over from Eq.18 to Eq.27 is likely to be difficult to observe in practice. This is particularly true if the equilibrium bound fraction is small, i.e., 1 << s p , for which the transient regime might well extend almost up to equilibrium. For example, in the experiments of Erickson et. al. that verified Berg-Purcell scaling [6] by effectively tuning the number of receptors, this fraction was below 0.5 in most cases. The authors also state that corrections due to the reduction in free receptor number due to binding were never observed, which lends further support to our assertion.
In order to see the differences in time evolution as predicted by Eq.13 and Eq.27 in two different time regimes, it is helpful to define the dimensionless parameter aD Nk . In these cases, the distinction between the intermediate and asymptotic regimes of evolution is unlikely to be sharp.
A precise estimation of the regime where the intermediate regime of evolution (Eq.34) crosses over to the asymptotic regime (Eq.35) may be obtained by demanding that the quadratic term in p(t) in Eq.12 is small compared to the linear term. The condition for applicability of Eq.13a (or Eq.34) then turns out to be To summarize, it is not surprising that the experiments of Erickson et. al. did not notice any difference between the effective on-rate constants measured from (i) the slope of initial rise of the data curve and (ii) by fitting the entire data curve to a theoretical function. In order to see the effect predicted in the present paper, it would be necessary to do an exponential curve-fitting at both ends of the binding curve, i.e., near 0 = t and near s p p ≈ .

Application in a simple computational model
To investigate how the altered association rate expression might impact ligand binding, we simulated a simple single receptor-single ligand binding system, based on epidermal growth factor (EGF) binding parameters for concreteness. The experimental parameter values [17 ] can be found in Table 1, and the complete set of equations used are given in Appendix D. The primary question was whether application of Eq. 33 would lead to differences in EGF binding curves when compared to the Berg-Purcell-Shoup-Szabo effective rate constants (Eq.14) or simply using the intrinsic rate constants Note that the dissociation rate constant was maintained at the intrinsic binding rate when coupled with Eq. 33 since effects of rebinding were pooled into that expression.
A slight difference in the fraction of receptor-ligand complexes formed as a function of time was evident when using Eq. 33, compared to the intrinsic association rate, and negligible difference when compared to using the Berg and Purcell expressions (data not shown). The maximal value of α for this system was ~0.2 , which suggests that the impact would be small here. The differences resulting from the time-dependent association rate expression are better illustrated when the diffusion coefficient is lowered by a factor of ten and the k + is increased by a factor of five, which resulted in 10 ≈ α ( Fig  2A, Table 1) At early times, the new expression (Eq 33) leads to similar kinetics to those of Berg and Purcell, despite differences in the dissociation rate terms, but quickly leads to kinetics mirroring that of the intrinsic rates. However, in all cases the same steady-state level of complexes was obtained. Note that although the intrinsic rate of dissociation was used in the simulations for the dissociation rate with Eq.33, the "effective" dissociation rate found from experiments (i.e. rate found if data was analyzed with a simple ligand-receptor binding model) would be altered from the intrinsic value due to rebinding which, in our analysis, is dealt with through the time-dependent association rate given by Eq.33 We then asked how the rate expressions would impact a more complex system where there were two ligands competing for the single receptor. Using a second "mutant" EGF ligand where the intrinsic association rate was set equal to EGF but the off-rate differed by an order of magnitude, we simulated a competition binding study at 4ºC with equal concentration of EGF and the higher affinity mutant (1 ng/ml) (Fig 2b). The effects are present using baseline EGF values (Table 1, data not shown) but, as with the single ligand system, are more prominent when a lower diffusion coefficient and higher association rate are used. As expected, the level of mutant EGF-EGFR complexes exceeds those of EGF-EGFR with all three rate expression types (intrinsic values, Berg and Purcell expressions, and Eq.33). The new expression matches the Berg and Purcell expression for early times for both mutant EGF complexes and EGF complexes but begins to track with the intrinsic rate in less than 10 min. All three expressions resulted in similar late time (> 60 min) complex levels in agreement with the single ligand system. There was some "overshoot" of the steady-state EGF-EGFR levels evident with all three expressions although none were large. We next simulated conditions where there was an excess of mutant EGF (5 ng/ml mutant to 1 ng/ml EGF) and asked how the unequal concentration impacts complex formation. Trends seen with the EGF mutant with regard to complex formation at equal concentrations added were found with the only difference being the significantly higher binding levels. Complex formation with EGF however was somewhat different in that the "overshoot" was much more pronounced (Fig 2c). In all cases the levels of complexes "overshot" the steady-state level with the new expression resulting in a larger "overshot" compared to either the intrinsic rate or the Berg and Purcell expression. This behavior was evident whether the mutant was of higher or lower affinity (data not shown). Clearly, the differences in the kinetic behavior between Eq.33 and the standard Berg-Purcell rates can be appreciable, and could be important in situations where kinetic measurements play a crucial role.

Discussion
Ligand binding to receptors on a spherical cell surface is a classic problem, and has been studied by several authors, starting with the monumental work of Berg and Purcell [2]. Although differing in details, much previous work has essentially dealt with calculating effective rate constants for the binding and dissociation of ligands [2][3][4][5]. The underlying idea was that the effects of diffusion limited transport of ligands and the nonhomogeneous distribution of receptors in space (i.e., their confinement to the cell surface) could be absorbed into these effective rates. The question of the detailed nature of the kinetic binding curve was not explicitly addressed, although it seems to have been generally assumed that once the intrinsic rates were replaced by these effective rates, the standard mean-field (space-averaged) rate equation would be sufficient to describe the kinetics. This assumption was questioned by Goldstein and Dembo [10] who pointed out that this assumption is of an ad-hoc nature. However, the technique they used, although more systematic than the ad-hoc approach, is perturbative in nature and fails to predict the correct behavior in all the time regimes.
In this paper, we have developed a self-consistent stochastic approach to the binding problem, which has been built upon a similar formalism presented by us earlier to study ligand dissociation and rebinding in the presence of receptor clusters [11]. The principal idea behind this approach is to look at the problem in terms of individual ligand trajectories and derive an expression for the Green's function of ligand diffusion in the presence of the semi-absorbing cell surface. The final expressions for the cell surface ligand concentration (Eq.7 and 8) involve only macroscopic, experimentally measurable parameters like the ligand association rate and the receptor surface density.
Our results confirm the validity of the commonly used quasi-steady state approximation and the associated Berg-Purcell scaling relations over a certain well-defined time regime. However, we also predict that deviations from this relationship will be evident at larger times, close to equilibrium. In particular, we predict that the effective association rate is asymptotically the same as its intrinsic value, but undergoes a non-monotonic change between the Berg-Purcell scaling regime and this late time. The dissociation rate, however, remains a fraction of its true value at all times (beyond a short initial transient period) due to rebinding effects. Our principal prediction is that, if ligands are lost from the surface only through dissociation (i.e., provided other mechanisms like receptormediated internalization can be neglected), the ligand concentration at the cell surface (in fact, anywhere in the solution) at equilibrium is the same as the ligand density at infinity. However, significant variations with time leading to this point do occur, and have been fully characterized in the theory (Eq.18).
Our results are of possible implications in many areas. The quasi-steady state scaling assumption for association and dissociation rates are widely used in computational modeling of biochemical reactions describing ligand-receptor interactions, and it is important to determine if the non-monotonic change in the on-rate and its eventual saturation at the intrinsic rate as predicted in this paper has a significant effect on these types of results. It would also be of interest to extend the ideas in this paper to cellsurface (two-dimensional) ligand-receptor reactions (where the ligand interacts with the receptor after absorption on the cell surface and reaching it by surface diffusion) and is an area we are currently pursuing. This could have direct relevance in the context of receptor clustering, as in, for example, lipid rafts [18]. The results might also be applicable to the question of gradient detection in chemosensory cells. Investigations in this direction are currently in progress and will be reported in the near future. In general, the selfconsistent formalism that we have presented might be found to be of use in many other problems where external ligands/molecules/monomers are absorbed and released simultaneously (eg. motor proteins attaching to and detaching from microtubule filaments).
which leads to Eq.8.

Appendix B
In this appendix, we will solve the pure dissociation problem, where there are no ligands in the bulk solution to start with, but a certain fraction p(0) of receptors are ligand-bound at The decay is now exponential, but the rate has been modified on account of the frequent rebinding events.  In the small s-limit, it is easy to see from Eq.C2 and Eq.C3 that the two density contributions approach the steady state values Ligand-receptor complexes for the single ligand-single receptor system are plotted as a function of time using variable association/dissociation rates: intrinsic k + and k -(triangle), Berg and Purcell analysis (square), and Eq.33 (circle), as described in Appendix D. B. Ligand-receptor complexes in the two ligand-one receptor system with equal concentrations of ligand added (EGF-EGF (filled) mutant EGF-EGFR (open)) are plotted as a function of time using variable association/dissociation ratesas described in A\. C. Ligand-receptor complexes in the two ligand-one receptor system when unequal concentrations (1 ng/ml EGF and 5 ng/ml EGF mutant) are added. Only EGF-EGFR complexes are shown and symbols are the same as listed in A. Equations and Parameters are listed in Appendix D and Table 1.