Equilibrium Probability Distribution for Number of Bound Receptor-Ligand Complexes
1 Equilibrium Probability Distribution for Number of Bound Receptor-Ligand Complexes
Tuhin Chakrabortty * , Manoj M Varma * Center for Nano Science and Engineering, Indian Institute of Science, Bangalore, India
Robert Bosch Centre for Cyber-Physical Systems, Indian Institute of Science, Bangalore, India
Abstract
The phenomenon of molecular binding, where two molecules, referred to as a receptor and a ligand, bind together to form a ligand-receptor complex, is ubiquitous in biology and essential for the accurate functioning of all life-sustaining processes. The probability of a single receptor forming a complex with any one of πΏ surrounding ligand molecules at thermal equilibrium can be derived from a partition function obtained from the Gibbs-Boltzmann distribution. We extend this approach to a system consisting of π receptors and πΏ ligands to derive the probability density function π(π; π , πΏ) to find π bound receptor-ligand complexes at thermal equilibrium. This extension allows us to illustrate two aspects of this problem which are not apparent in the single receptor problem, namely, a) a symmetry to be expected in the equilibrium distribution of the number of bound complexes under exchange of π and πΏ and b) the number of bound complexes obtained from chemical kinetic equations has an exact correspondence to the maximum probable value of π from the expression for π(π; π , πΏ) . We derive the number fluctuations of π and present a practically relevant molecular sensing application which benefits from the knowledge of π(π; π , πΏ) . Introduction
The phenomenon of molecular binding, where two molecules, generally referred to as a receptor and a ligand, bind together to form a ligand-receptor complex, occurs ubiquitously in biological systems and is essential for the accurate functioning of many life-sustaining processes ranging from DNA transcription to immune response [1, 2]. The ligand-receptor complex formation is a reversible process, which means that while a ligand molecule can associate with a receptor molecule to form a ligand-receptor complex, the converse process, namely the dissociation or breaking down of a ligand-receptor complex into one each of the ligand and receptor molecules can also happen simultaneously. In general, the rates of these two competing processes would be different leading to one process dominating over the other. As we will see in the next section of this article, thermodynamic equilibrium is reached when the rates of these two processes become equal. The equilibrium state is also a steady state where the number of the bound ligand-receptor complexes remains constant over time due to the two equal and opposite processes cancelling each other. While living organisms are non-equilibrium systems which exchange energy and matter with the surroundings, as illustrated beautifully by Rob Philips, the separation of time-scales allows a large number of biological processes involving molecular binding to be viewed as processes in equilibrium [3]. Due to this ubiquitous and critical role of molecular binding in biological functions, we would like to be able to answer questions such as, how many ligand-receptor complexes do we expect to find at equilibrium after we start with π receptors and πΏ ligands. As we will see in the next paragraph, there are a number of different ways to answer this question and related questions such as the magnitude of fluctuations in the number of ligand-receptor complexes at equilibrium. One of the ways to answer the question on the number of ligand-receptor complexes to be expected at equilibrium, starting from π receptors and πΏ ligands can be answered by considering the molecular binding as a reversible chemical reaction and writing down the rate equations by invoking the law of mass action[4]. This approach is worked out in section two of this article. This approach provides not only the number of ligand-receptor complexes at equilibrium but also how the number of complexes change as a function of time from the initial value of zero. A question then arises, namely, if we perform π independent experiments each with π receptors and πΏ ligands and measure the number of ligand-receptor complexes at equilibrium for each of these experiments, will we always get the same number as predicted by the reaction rate approach? As we might expect, the answer is no. There 2 will be variation in the number of ligand-receptor complexes seen in the different experiments which can be quantitatively characterized by the variance or standard deviation of the measured values. The variation in the measured number of complexes arises not just because of experimental uncertainties which may be present but also due to a more fundamental reason, namely, that the binding of a receptor and a ligand molecule to form the complex as well as its dissociation at a later point in time are stochastic events which follow certain probability distributions. As a result, even if the experimental uncertainties were reduced to zero hypothetically, one would still find a distribution of values for the number of ligand-receptor complexes even if one performed identical experiments repeatedly. Therefore, an estimate for the magnitude of the variation is essential to obtain a fundamental understanding of noise in biological systems [5, 6, 7, 8] or in a more applied context to understand and optimize performance of biosensors. Biosensors are devices which rely on molecular binding to detect or quantify the presence of specific molecules in human samples such as blood. The amount of these specific molecules, referred to as biomarkers, indicate the presence of diseases such as cancer [9]. One may like to know the fluctuations in the biosensor signal caused due to the stochastic nature of molecular binding in diagnostic applications, which are generally performed under equilibrium conditions. The magnitude of these fluctuations set a fundamental limit on the lowest While the reaction rate approach can tell us the number of ligand-receptor complexes at equilibrium, it cannot give us any information regarding the statistical distribution of this quantity and related parameters of interest such as the magnitude of fluctuations (quantified by parameters such as the variance or the standard deviation) around the mean number of complexes. It is in this context that the application of statistical physics to this problem helps us answer these questions in detail. We know that equilibrium statistical mechanics allows us to calculate the distribution of particles into the various states of a given system. These states may be discrete or continuous and are characterized by a specific energy level of the state. By enumerating all the states of a given system appropriately scaled by their thermodynamic weights, one can obtain the partition function of the system which then yields the probability distribution of particles in these states. In the molecular binding problem, one could consider the unbound ligand and receptor molecules as one state and the bound complex as another state forming a two-state system and calculate the partition function. This approach has been exactly followed to derive the equilibrium occupation probability of a single receptor surrounded by πΏ ligand molecules [3, 6, 10, 11]. The equilibrium occupation probability of the receptor is the probability that the receptor will be bound by one of the πΏ ligand molecules in equilibrium. The resultant formula is shown to be equivalent to the one which is obtained from the chemical reaction rate equations [12]. This problem serves as a good illustration of how statistical mechanics links macroscopic parameters such as the equilibrium dissociation constant to the microscopic, molecular level details such as the energy level difference between bound and unbound states. Although the problem with one receptor and πΏ ligands serves to illustrate the application of equilibrium statistical mechanics to solve the binding problem, in many real situations we have more than one receptors. For instance, a cell may have multiple receptors on its surface and biosensor certainly has millions of receptors which are used to bind a ligand target. Therefore, we would like to extend the solution of the one receptor problem into one which has π receptors. The formal statement of the problem is to calculate the probability distribution of the number of receptor-ligand complexes π , under thermal equilibrium conditions, given π receptors and πΏ ligands. If the binding of ligands to receptors can be considered as independent events, then the answer is simply π times the probability of binding a single receptor by πΏ ligands which has been solved previously. The assumption of independence is a little tricky because the occupation of one receptor by a ligand implies that for the next complex formation there is one less receptor and one expects that the probabilities of binding would change. However, are there situations where independence of events can be justified? To answer such questions systematically, one way forward is to write the partition function for the system and obtain the required probability distribution. An alternate method of deriving the 3 probability distribution of receptor-ligand complex is from the master equation describing the binding process which is Markovian, and accounting for transitions into and out of the state with π receptor-ligand complexes [13]. One can obtain the partition function from a recurrence relation extracted from the master equation. The advantage of this method is that in addition to the equilibrium probability distribution of the number of complexes, one can also obtain the temporal evolution of this distribution. Here we illustrate the first method of solution. In this article we derive the partition function by direct enumeration and from it the probability density function π(π; π , πΏ) , which is the probability to find π bound receptor-ligand complexes at thermal equilibrium with the initial receptor and ligand concentrations of π and πΏ respectively. The article is divided into the following sections. In section 2, we derive the number of bound complexes π *+ , at equilibrium using the conventional approach based on balancing of forward and backward reaction rates. Here we show that the expression for the number of bound receptors from the binding kinetic equation (π *+ ) is symmetric with respect to the exchange of π and πΏ and discuss why this symmetry is to be expected in this problem. In section 3, we present the partition function for a system consisting of π receptors and πΏ ligands which can bind to form complexes. We show that this expression and the probability density function π(π; π , πΏ) which follows from it is also symmetric (invariant) under exchange of π and πΏ . In section 3, we provide explicit expressions for π(π; π , πΏ) and show that the expression for π *+ derived from the reaction rate approach corresponds to π - which maximizes π(π; π , πΏ) , and not to the expectation value β©πβͺ . This indicates that π(π; π , πΏ) is a skewed distribution, although the skew is not very significant. In section 4, we show that under certain regimes π(π; π , πΏ) , which is a rather complicated expression, can be approximated by a binomial distribution. We then derive the number fluctuations of the bound complex and finally in section 5, we connect this problem to a practical situation involving the sensing of biomolecules in the context of the extremely important application of early diagnosis of disease conditions.
2. The reaction rate approach
Let us consider a system with π receptors and πΏ ligands. The ligands can bind to the receptors reversibly with binding rate π and unbind with a rate π . The general reaction kinetic equation specifying the rate of formation of complexes π follows from the law of mass action as, ππππ‘ = π (π β π)(πΏ β π) β π π (1) Here, the first term in the RHS represents the instantaneous rate of free ligand and free receptor molecules associating to form a complex while the second term represents the instantaneous rate of dissociation of existing ligand-receptor complexes. At equilibrium, these two opposing reaction rates are equal and the RHS becomes equal to zero from which we see that the number of complexes π *+ at equilibrium satisfies the quadratic equation below, π
8π β π *+
98πΏ β π *+ π *+ (2) Solving Eq. (2) for π *+ , we get π *+ = (π + πΏ + πΎ < ) Β± >(π + πΏ + πΎ < ) ? β 4π πΏ2 , (3) where πΎ < = B CDD B CE .
4 Equation (3) describes the general solution for the number of bound complexes at steady state indicating two possible values of π *+ at equilibrium. However, we know that π *+ cannot exceed the maximum of πΏ or π . Therefore, the only physically acceptable solution of π *+ is, π *+ = (π + πΏ + πΎ < ) β >(π + πΏ + πΎ < ) ? β 4π πΏ2 (4) For special cases e.g. π βͺ πΏ and
π β« πΏ , one can modify Eq (2) by ignoring the depletion in the number of receptor and ligand molecules as shown in Eq. 5(a) and Eq. 5(b) respectively, π
8π β π *+
9πΏ = π π *+ π π 8πΏ β π *+ π *+ (5a) (5b) which gives us the familiar solutions for the reaction kinetic equation known as the Michaelis-Menten equation. [14] π *+(JβͺK) = π πΏπΏ + πΎ < π *+(Jβ«K) = π πΏπ + πΎ < (6a) (6b) One of the noteworthy differences between Eq (4) and Eq (6), besides their different forms, is symmetry with respect to exchange of π and πΏ . The form of the solution (Eq (4)) is symmetric with respect to π and πΏ , which is expected given that the receptors and ligands are arbitrary labels. Therefore, interchanging the labels should not affect the final expression. On the other hand, for the special cases (Eq (6)) the asymmetry comes in because we ignore the depletion of one of the species, thereby distinguishing the labels. A symmetric expression similar to the form of equation (6) can be obtained by imposing the condition < ) ? in equation 4 as π *+(LM-) = π πΏπ + πΏ + πΎ < , (7) which reduces to equation (6) for π β« πΏ or π βͺ πΏ . Figure 1 shows the effects of the approximations in Eq (6) and Eq (7) in different regimes. We define relative error πΏ as πΏ = π *+ β π *+(OPPQ1R) π *+ where π *+ is the general form of the number of bound receptors described in Eq (4) and π *+(OPPQ1R) represents the three approximated forms for the number of bound receptors, namely π *+(JβͺK) , π *+(Jβ«K) and π *+(LM-) . The relative error is plotted with respect to the number of receptor molecules π while the ligand concentration πΏ and the dissociation constant πΎ < are kept constant. Therefore, changing π from π βͺ πΏ to π β« πΏ covers all the possible regimes. As expected, while Eq 6(a) and Eq 5 6(b) are very good approximations for Eq (4) for regimes
π βͺ πΏ and
π β« πΏ respectively, the relative errors for both the approximations are quite significant for other regimes. Interestingly, the symmetric expression π *+(LM-) is the most accurate approximation across all regimes. Figure 1.
Effect of approximations: The relative error for the three approximations namely, *+ JβͺK , 8π *+ Jβ«K and 8π *+ LM- are plotted with respect to the number of receptors for all the regimes. The symmetric approximation
V8π *+ LM- W shows the least error across all regimes. The relative error plots of *+ JβͺK and 8π *+ Jβ«K do not diverge symmetrically because the relative error for of *+ JβͺK becomes independent of R for smaller values. [Parameters:
πΏ = 1000; πΎ < = 100 ]
3. Partition function for a system containing πΉ receptors and π³ ligands The number of bound complexes at equilibrium can also be calculated from equilibrium statistical mechanics using Gibbs-Boltzmann distribution. As before, we consider a system with πΏ ligands uniformly distributed in the solution and π receptors attached to a surface. We conceptualize the π receptor πΏ ligand problem in a similar manner as the one receptor πΏ ligand problem solved previously [10]. We discretize the system containing the receptor and ligand molecules in solution as being composed of Ξ© boxes or sites as shown in figure 2. We imagine that the discretization of the system volume into these boxes is such that any given box or site can be occupied by at most one free ligand or receptor molecule or a single ligand-receptor complex. This means that out of a total of Ξ© boxes, π of them would constitute the ones which contain the π receptor molecules present in the system, These are pictorially indicated by the blue colored boxes in figure 2. The ligand molecules, represented by the red discs in figure 2, can be distributed into any of the total Ξ© boxes. Note that a ligand molecule may occupy a colorless (empty) box or blue colored box. These two options represent the two states that a ligand molecule may be found in, namely, as a free molecule (if occupying a colorless box); and as a ligand-receptor complex (if occupying a blue colored box). These two states, namely, the free and bound states of the ligand molecule, are associated with energies π and π ^ respectively, with the sub-scripts π and π denoting the free and bound states. The system may be found in different configurations each with a different total energy. For example, figure 2 shows a specific configuration for a system with π = 7 and
πΏ = 12 and number of complexes π = 3.
The total energy πΈ of the system with the specific configuration shown in figure 2 would be πΈ = 9π + 3π ^ , because out of the 12 6 ligands, 9 are in the free state (occupying colorless boxes) and 3 are in the bound state (occupying blue boxes where a receptor is already present). However, as we can clearly see, there can be many more configurations which have the same value of total energy πΈ because we could rearrange the red discs into any of the other boxes which are not occupied in figure 2. Statistical mechanics provides us with a recipe to calculate the probability of finding the overall state of the system, for instance with π = 3 as shown in figure 2, from the number of configurations, referred to as microstates, leading to the same overall system state or total energy. Figure 2, A schematic representation of the microstates corresponding to a state of the system. The system here is composed of 7 receptors and 12 ligands. A particular configuration corresponding to the state of the system with 3 ligand-receptor complexes is shown here. Other possible combinations are rearrangements of the objects as described in the main text which are enumerated to obtain the partition function. Let us calculate the number of microstates for the state of the system, similar to the one shown in figure 2, where the number of receptors, ligands and complexes are π , πΏ and π respectively. We can do this by calculating the number of ways πΏ indistinguishable ligands can be distributed among the Ξ© β π free sites (colorless boxes in figure 2) and the π binding sites (blue boxes). The total energy of this state would be πΈ = (πΏ β π)π + ππ ^ as there are πΏ β π ligands which are free and π ligands which are bound, with energy levels π and π ^ respectively. As already explained in solution to the one receptor problem, Ξ© β« π , πΏ . Therefore,
Ξ© β π can be well approximated by Ξ© itself. Thus, the number of arrangements or configurations of the state of the system with total energy πΈ = (πΏ β π)π + ππ ^ is simply the number of ways to distribute πΏ β π ligands into
Ξ© β π β Ξ© free sites and π ligands into the π sites occupied by the receptors, which is simply eKfQ JQ where VππW = . The partition function π is then obtained by multiplying this number of microstates with the so-called Boltzmann weight π fkl , where π½ = nB o p with π q the Boltzmann constant, and summing over all possible values of π . The smallest possible value of π is clearly zero (no ligand-receptor complex) and the maximum possible value of π will be the smaller of π and πΏ . So the partition function becomes π = r s Ξ©πΏ β πt uvw(K,J)Qxy π f(KfQ)kz D sπ πt π fQkz { (8) Equation (8) can be further simplified by invoking Ξ© β« (πΏ β π) and therefore
Ξ©!8Ξ© β (πΏ β π)9! = Ξ©
KfQ . (9) Introducing Eq (9) into Eq (8), the expression for π can be simplified as 7 π = 8Ξ©π fkz D K r 1(πΏ β π)! uvw(K,J)Qxy sπ πt 8Ξ©π k|z fQ , (10) where Ξπ = π ^ β π is the binding energy of the receptor-ligand complex. The probability that there are π bound complexes in the system can be written in terms of the partition function as, π(π; πΏ, π ) = 1π ~8Ξ©π fkz D K (πΏ β π)! sπ πt 8Ξ©π k|z fQ (cid:127) (11) Introducing the expression for π as derived in Eq (10) and eliminating the common terms, Eq(11) gives us π(π; πΏ, π ) = π ! πΏ!(π β π)! (πΏ β π)! π! 8Ξ©π k|z fQ β π ! πΏ!(π β π (cid:129) )! (πΏ β π (cid:129) )! π (cid:129) ! (Ξ©π k|z ) fQ (cid:130) uvw(K,J)Q (cid:130) xy (12) where π (cid:129) is the dummy variable used as the summation index. Equation (12) describes the general form of π(π; πΏ, π ) which is easily seen to be symmetric with respect to exchange of πΏ and π and hence consistent with the reaction kinetics-based results mentioned in section 1. Further, the partition function and the probability distribution shown in Eq. (11) is identical to that obtained from the master equation-based approach as expected [13]. Obtaining a closed form expression for π(π; πΏ, π ) in Eq. (12) is rather difficult as illustrated in Appendix I. However, we can obtain explicit expressions for π(π; πΏ, π ) for special cases where π β« πΏ or π βͺ πΏ . For instance, take the case where
π βͺ πΏ as considered commonly in biosensor literature [15, 16, 17]. As
π βͺ πΏ and maximum of π is the minimum of πΏ and π , it follows that πΏ β« π allowing us to make the approximation,
πΏ!(πΏ β π)! = πΏ Q Therefore, the partition function for this case (π JβͺK ) becomes π JβͺK = 8Ξ©π fkz D K (1 + πΌ K ) J πΏ! where πΌ K = Ke* (cid:132)(cid:133)(cid:134)
The probability mass function π JβͺK (π; πΏ, π ) can then be readily obtained as π JβͺK (π; πΏ, π ) = 1(1 + πΌ K ) J (cid:135)sπ πt πΌ KQ (cid:136) (13) Similarly, the probability mass function for π β« πΏ can be written by interchanging the index πΏ with π in Eq (13) as 8 π Jβ« K (π; πΏ, π ) = 1(1 + πΌ J ) K (cid:135)sπΏπt πΌ JQ (cid:136) (14) where πΌ J = Je* (cid:132)(cid:133)(cid:134)
Unlike Eq. (12), the expressions in Eq (13) and Eq (14) are no longer symmetric with respect to the exchange of π and πΏ because we are ignoring the depletion effect of one of the two species. By re-arranging the probability mass functions described in Eq. (13) as shown in Eq. (15), we see that π JβͺK (π; πΏ, π ) is a Binomial distribution of the form B (1 β π) with π = π ; π = π and π = (cid:137) (cid:138) n(cid:139)(cid:137) (cid:138) π JβͺK (π; πΏ, π ) = sπ π t s πΌ K K t Q s 11 + πΌ K t JfQ (15) The term (πΏ/Ξ©) is the concentration of ligand molecules in the solution and the term π k|z can be mapped to the macroscopic dissociation constant πΎ (cid:141) [10]. For a system of volume π£ , the dissociation constant πΎ < described in section 2 can be written in terms of the macroscopic dissociation constant πΎ (cid:141) as πΎ < = π£πΎ (cid:141) . Introducing these relations in Eq (14), the expression for π JβͺK (π; πΏ, π ) can be obtained as, π JβͺK (π; πΏ, π ) = sπ πt s πΏπΏ + πΎ < t Q s πΎ < πΏ + πΎ < t JfQ (16) The Binomial distribution
π(π) = 8 B (1 β π) for π β β and π β 0 with the product ππ equal to a finite value π , can be approximated by a Poisson distribution given by π(π) = (cid:147) (cid:148)(cid:149) * (cid:148)(cid:150) B! . Here, the quantity π is the expectation value or the mean of the distribution denoted by β¨πβ© . The number of receptor molecules in many situations is fairly large (of the order of 10 or more; see section 5 for example) and we may consider π β β . Further if πΎ < β« πΏ , then π = KK(cid:139)(cid:153) (cid:154) is close to zero. The product of π and π is π = JKK(cid:139)(cid:153) (cid:154) β JK(cid:153) (cid:154) . Similarly, for
π βͺ πΏ , one can also show that the above argument is valid if
πΏ β β and πΎ < β« π . Therefore, the Poisson approximation of the Binomial distribution will be given by π (cid:153) (cid:154) β«J,K (π; πΏ, π ) = Vπ πΏπΎ < W Q π fJK(cid:153) (cid:154) π! (17) Equation 17 is commonly used in the case of multistep biological systems such as DNA transcription to simplify the analysis [18]. The expression is also relevant for describing systems where the number of bound complexes is extremely low compared to both ligand and receptors e.g. interaction of non-specific ligand molecules and cell surface receptors [12]. The steady state distribution described in Eq 17 can also be derived from the chemical master equation (CME) as described in reference [18]. To summarize, in this section we showed that under different assumptions, the number of bound complexes at equilibrium, starting from π receptors and πΏ ligands follows different probability distributions. Table 1 summarizes the analysis of this section. 9 Condition Binding kinetic equation Equilibrium Statistical Mechanics
π βͺ πΏ ππππ‘ = π (π β π)(πΏ) β π π π *+ = π πΏπΏ + πΎ < π JβͺK (π; πΏ, π ) = sπ πt s πΏπΏ + πΎ < t Q s πΎ < πΏ + πΎ < t JfQ β¨πβ© = π πΏπΏ + πΎ < π - = π πΏπΏ + πΎ < (Binomial distribution) π β« πΏ ππππ‘ = π (π )(πΏ β π) β π π π *+ = π πΏπ + πΎ < π Jβ«K (π; πΏ, π ) = sπΏπt s π π + πΎ < t Q s πΎ < π + πΎ < t KfQ β¨πβ© = π πΏπ + πΎ < π - = π πΏπ + πΎ < (Binomial distribution) πΎ < β« π , πΏ ππππ‘ = π π πΏ β π π π *+ = π πΏπΎ < π (cid:153) (cid:154) β«J,K (π; πΏ, π ) = VπΏπ πΎ < W Q π fKJ(cid:153) (cid:154) π! β¨πβ© = π πΏπΎ < π - = π πΏπΎ < (Poisson distribution) π β πΏ ππππ‘ = π (π β π)(πΏ β π) β π π π *+ = (π + πΏ + πΎ < ) β >(π + πΏ + πΎ < ) ? β 4π πΏ2 π(π; πΏ, π ) = π ! πΏ!(π β π)! (πΏ β π)! π! 8Ξ©π k|z fQ β π ! πΏ!(π β π (cid:129) )! (πΏ β π (cid:129) )! π (cid:129) ! (Ξ©π k|z ) fQ (cid:130) uvw(K,J)Q (cid:130) xy β¨πβ© = πΏ + πΎ < πΏ π (1 β πΏ, βπΏ + π + 2, βπΎ < )π (βπΏ, βπΏ + π + 1, βπΎ < ) π - = (π + πΏ + πΎ < ) β >(π + πΏ + πΎ < ) ? β 4π πΏ2 Table 1: Landscape of receptor-ligand interactions: We have compiled the results obtained from the binding kinetic equation and equilibrium statistical mechanics for all the possible conditions. For the statistical mechanics method, we have the equilibrium probability distribution along with the mean (β¨πβ©) and the mode (π - ) of the distribution which is derived in the article. One interesting point to 10 observe here is that for all the conditions, the expression for the equilibrium solution for the binding kinetic equation *+ is equivalent to the mode (π - ) of the distribution.
4. Moments of π(π; π³, πΉ) : How are the probability distributions of the number of complexes summarized in Table 1 related to the number of complexes in Eq. 3 and Eq. 6, derived from the reaction rate approach? To answer this question, we need to calculate the moments of the distributions. Moments are single numbers which characterize an aspect of the distribution, for instance the first moment (see equation below) of π(π; πΏ, π ) gives us the expectation value or the mean number of bound receptors from the probability mass function. Let us calculate the expectation value for the case
π βͺ πΏ , denoted by the symbol β¨πβ©
JβͺK from the probability mass function provided in Table 1. We can write, β¨πβ©
JβͺK = r π
JQxy π JβͺK (π; πΏ, π ) = r π
JQxy sπ πt s πΏπΏ + πΎ < t Q s πΎ < πΏ + πΎ < t JfQ
This sum can be readily evaluated or one can use the standard result that the mean of the Binomial distribution
π(π) = 8 B (1 β π) is ππ to yield, β¨πβ© JβͺK = π πΏπΏ + πΎ < , (18) We see that Eq. (18) is identical to Eq. 6(a) obtained from reaction rate approach. Similarly, one can also calculate β¨πβ© Jβ«K by interchanging πΏ and π and see that the expectation value of π from the probability distribution is identical to the number of complexes at equilibrium obtained from the reaction rate. What about the Poisson approximation mentioned in Eq. (17)? We will only discuss the case π βͺ πΏ as the solution to the
π β« πΏ case can be obtained by simply exchanging π and πΏ . As discussed in the previous section, Poisson approximation is valid when πΎ < β« πΏ . This means that the off-rate is significantly high compared to the depletion effect of the ligands, i.e. the ligand-receptor complexes break down into free ligands and free receptors rapidly and therefore, one may assume that π receptors and πΏ ligands are always available for association. The rate of formation of complexes ignoring the depletion of receptors and ligands reads ππππ‘ = π π πΏ β π π (19) From which, equating the LHS to zero at equilibrium gives the equilibrium number of complexes as π *+ = π πΏπΎ < (20) Eq. 20 is identical to the mean of the Poisson distribution. So for all these cases we see that the mean of the probability distribution of the number of ligand-receptor complexes is identical to the number of ligand-receptor complexes derived from the reaction rate approach. Let us now calculate the mean of the distribution for the general case described in equation 12. For the general case, the expression for the mean number of bound complexes can be written as 11 β¨πβ© = r π uvw(K,J)Qxy π(π; πΏ, π ) = β π π ! πΏ!(π β π)! (πΏ β π)! π! 8Ξ©π k|z fQuvw(K,J)Qxy β π ! πΏ!(π β π)! (πΏ β π)! π! (Ξ©π k|z ) fQuvw(K,J)Qxy (21) From the calculations described in Appendix I, β¨πβ© can be written as β¨πβ© = πΏ + πΎ < πΏ π (1 β πΏ, βπΏ + π + 2, βπΎ < )π (βπΏ, βπΏ + π + 1, βπΎ < ) (22) where π(π, π, π§) is the confluent hypergeometric function of the second kind. Interestingly, the form of Eq. (22) is quite different from the result obtained using the rate equation in Eq (3). In Fig. 3a we have plotted the expression for β¨πβ© from Eq. (22) and the number of equilibrium ligand-receptor complexes from the reaction rate approach from Eq. 3 to numerically compare them. It is clear that while these expressions are very different algebraically, they match each other quite well numerically. The two expressions diverge marginally for
πΏ β π . Figure 3: a.
Comparison of the general forms of the solution obtained from the binding kinetic equation and equilibrium statistical mechanics. For special cases i.e. |πΏ β π | β« 0 , both the solutions match perfectly. However, for
πΏ β π , there is some difference between them as shown in the subfigure, particularly for low values of πΎ < . The difference reduces with increase in πΎ < . b. Probability distribution for different conditions which are shown by coloured circles in Fig a. Irrespective of the shape of the distribution, the mean and mode are very close to each other. [Parameters:
πΏ = 10 ] While the number of ligand-receptor complexes calculated using reaction kinetics approach numerically matches the expectation value obtained from the probability mass function of the number of complexes, the algebraic expressions are very different. Is there any other quantity derivable from the probability mass function of the ligand-receptor complexes which exactly matches the solution form the reaction rate approach? Indeed, we see that the mode of π(π; πΏ, π ) , denoted by π - , exactly matches Eq. (3) derived from the reaction rate approach. To calculate the mode of a distribution, we need to maximize the distribution with respect to π . We maximize log8π(π; πΏ, π )9 taking advantage of the monotonicity of the log function. Introducing log in equation (11), we get a b log π(π; πΏ, π ) = log sπ ! πΏ!π t β log8(π β π)!9 β log8(πΏ β π)!9 β log(π!) β π log(Ξ©π k|z ) (23) Using Stirlingβs approximation log π₯! = π₯ log π₯ β π₯ in equation (23) and maximizing with respect to π gives us log Β₯(πΏ β π - )(π β π - )π - Ξ©π k|z Ζ = 0 (24) where π - is the mode of the probability distribution. As explained before, Ξ©π k|z can be mapped to the dissociation constant πΎ < described in section 2. Therefore, Eq (24) becomes (πΏ β π - )(π β π - ) = π - πΎ < , (25) which can be solved for π - to obtain π - = (π + πΏ + πΎ < ) Β± >(π + πΏ + πΎ < ) ? β 4π πΏ2 (26) Similar to π *+ , the upper bound of π - is max(πΏ, π ). Therefore, we only consider the negative quadratic to be the physically acceptable solution. Therefore, π - becomes π - = (π + πΏ + πΎ < ) β >(π + πΏ + πΎ < ) ? β 4π πΏ2 (27) Therefore, the number of bound receptors at steady state π *+ derived in section 1 is the number of bound receptors one is expected to obtain experimentally with maximum probability and not the average number of bound receptors which is commonly used in the literature [18]. Then, why does the number of bound complexes obtained from reaction rate equations match the expected value β¨πβ© obtained from the probability distribution function π(π; πΏ, π ) numerically even though their analytical expressions are quite different? To answer this, we note that all the distributions in Table 1, including special cases such as π βͺ πΏ have their mode close to the mean irrespective of asymmetries in the distributions. Figure 3b illustrates the point mentioned above.
5. Higher moments of π(π; π³, πΉ) : One advantage of the equilibrium statistical mechanics method over the reaction rate approach is that it allows us to calculate all the moments of the probability distribution function at equilibrium. Having established that π JβͺK (π; πΏ, π ) is the binomial distribution for the special case of
π βͺ πΏ considered here, we can readily find the variance in the number of bound complexes, π JβͺK ? , as π JβͺK ? = π πΏπΎ < (πΏ + πΎ < ) ? (28) 13 Similarly, for π β« πΏ , the variance can be obtained by interchanging πΏ and π in Eq 28. Although, obtaining closed form expressions for the higher moments for the general form of π(π; πΏ, π ) is not possible, we can obtain a general expression for the higher moments. To estimate the higher moments, we define a moment generating function πΊ(π§) = r π§ Q π(π; πΏ, π ) uvw(K,J)Qxy (29) The π Β«βΉ moment π(π) of π(π; πΏ, π ) can be obtained by calculating the π Β«βΉ partial differential of πΊ(π§) at π§ = 1 . The π Β«βΉ partial differential of πΊ(π§) at π§ = 1 gives us π ο¬2 πΊ(π§)| ο¬xn = r (cid:176)β(π β π) π(π; πΏ, π )Β· uvw(K,J)Qxy (30) Expanding Eq (30), we get π ο¬2 πΊ(π§)| ο¬xn = r π(π; πΏ, π ) uvw(K,J)Qxy ~π β r ππ + r r πππ β¦ β¦ (cid:127) (31) In literature, the coefficient of π B in equation (31) is called Stirlingβs number of first kind and is denoted as π (π, π) [19]. Therefore, Eq (32) can be rewritten in terms of Stirlingβs numbers as π ο¬2 πΊ(π§)| ο¬xn = r π(π; πΏ, π ) uvw(K,J)Qxy r π (π, π)π
B2Bxy (32) The π Β«βΉ moment π(π) now can be calculated as
π(π) = π ο¬2 πΊ(π§)| ο¬xn β r π(π β π)π ο¬β‘ πΊ(π§)| ο¬xn2fnβ‘xn (33) where π(π) = π (π, π) β β π(π)π (π β π, π β π)
Bfnβ‘xn with π(0) = 1, and π(1) = π (π, 1)
6. A numerical example: Micro/Nano particle-based molecular sensors
Early, pre-symptomatic, detection of cancer is an unmet demand which can potentially save millions of lives annually [20]. Detection of a class of molecules called cytokines has been pursued as a means to achieve early detection of cancer [21]. In such sensors, a sample volume π L containing the target ligands interacts with a system containing receptors and the number of bound complexes translates 14 into a physically measurable quantity such as a change in electrical conductance or optical absorbance [22]. For a designer of such a molecular sensing system, it is important to understand the effect of design choices such as the area of the sensing region used to capture target ligand molecules from sample, the volume of the sample to be used etc. on the noise due to the stochastic nature of molecular binding. Using what we learnt about the variance of molecular binding in the previous section, we can begin to answer such questions an example of which is provided below. Figure 4, Schematic representation of a molecular sensing device. Receptors immobilized in a sensing area bind target ligands in a sample. Let us consider a sensor system, as schematically shown in figure 4, consisting of a sensor surface of radius π where receptor molecules are attached. Typically, the surface of the particle is saturated with a monolayer of receptor molecules. Therefore, the number of receptor molecules on each bead can be estimated to be π = β¦O β° (cid:141) β° where π is the size of a single receptor molecule which we can take to be of the order of 10 nm [23]. The number of ligand molecules is πΏ = ππ L where π is the concentration of ligand molecules in the sample. Given these conditions, how should one choose the area of the sensing region (size π in figure 4) so that the fluctuation in the number of the bound receptors is less than 1% of the average number of the bound receptors? For a real system, π L is typically of the order of 100 Β΅ L. Taking π to be 10 pico-molar (pM), the number of ligand molecules in the system becomes 6x10 . The mean number of bound receptors, Β΅ , and the variance in the number of bound receptors, Ο ? , on the sensing region can be obtained as, Β΅ = π πΏπΏ + πΎ < Ο ? = π πΏπΎ < (πΏ + πΎ < ) ? (34a) (34b) We define relative fluctuation per particle as Ο (cid:192)`Β΄ = ΛΛ as Ο (cid:192)`Β΄ = Β―πΎ < π πΏ (35) Substituting the expression for π and πΏ from above, we get, π JlK = Β― πΎ < π ? Οπ ? ππ L (36) 15 We take the value of πΎ < to be around 6x10 , corresponding to a dissociation constant of 100 pM. If the desired fluctuation of the bound number of receptors (which will generate the sensor signal) is to be less than 1%, from Eq. (37), we see that π > Β― πΎ < π ? Οπ JlK ππ L Putting the numerical values we see that the size of the sensing region must be greater than about 200 nm to achieve a relative fluctuation of less than 1% in the number of bound receptors. Why does this happen? This is because reducing the size of the sensing region decreases the number of receptors available in the sensing region and increases the relative fluctuations as seen from Eq. (36). This exercise immediately tells us that extreme miniaturization as seen in single molecule level sensors such as those based on nanopores or nano field-effect transistors [24, 25] are likely to suffer from large relative fluctuations due to the small number of receptors available.
7. Conclusions
In this article, we derived the equilibrium probability distribution of the number of ligand-receptor complexes in a system consisting of π receptors and πΏ ligands. This work is an extension of a similar problem with one receptor and πΏ ligands which is now widely available in the references provided in this article. However, a treatment of the general π receptor case and the thorough description of the connections between equilibrium probability distributions obtained and their connections with the number of complexes obtained using the reaction rate approach under various special cases is lacking in current literature. We hope that this article has addressed this gap to some extent. For instance, in Table 1, we summarized different situations relevant to molecular sensing and compared the equilibrium statistical mechanics and the reaction rate approaches and showed how these relate to one another, thus allowing the student to appreciate the equivalence of these methods. Further, this article helped us to illustrate two aspects of this problem which are not directly apparent in the solution to the single receptor problem, namely, a) a symmetry of the equilibrium probability distribution of the number of bound complexes under exchange of π and πΏ and b) the number of bound complexes obtained from chemical kinetic equations has an exact correspondence to the maximum probable value of π . Overall, we hope that this article introduces the application of statistical mechanics to molecular binding problems in a broader setting than previous works. Appendix: Expression for π(π; π³, πΉ) in the general case
From Abramowitz and Stegun [19] Eq 13.5.2, we obtain
π(π, π, π§) = π§ fO Β¨r (π) Q (1 + π β π) Q π! (βπ§) fQ + π(|π§| fp ) pfnQxy Λ (A1) Here π(π, π, π§) is the confluent hypergeometric function of second kind and (π₯) is the Pochhammer symbol defined as (π₯) = Ξ(π₯ + π)Ξ(π₯) , Ξ(π₯) is the Gamma function. 16 The summation in Eq (A1) equivalent to the partition function π defined in Eq (9) for π = min(πΏ, π ) +1, π = βπΏ, π = 1 + π β πΏ, and π§ = βπΎ < . The error term can be neglected for πΎ < β« 0 or min(πΏ, π ) β« 0 . Introducing the parameters in Eq (A1) and using the property (βπ₯) = (β1) (π₯ β π + 1) , we obtain π = r s πΏ!(πΏ β π)! π!t s π !(π β π)!t (πΎ < ) fQuvw(K,J)Qxy = (βπΎ < ) fK π(βπΏ, 1 + π β πΏ, βπΎ < ) (A2) To derive the expression for the mean of the distribution (π) , we use the moment generating function πΊ(π) defined in Eq (28) in the main text, which becomes
πΊ(π) = π K π VβπΏ, 1 + π β πΏ, β πΎ < π Wπ(βπΏ, 1 + π β πΏ, βπΎ < ) (A3) π can now be obtained from Eq (A3) as π = ππΊ(π)ππ | βxn Using the relation [19] πππ§ π(π, π, π§) = βπ π(π + 1, π + 1, π§) we obtain π as π = πΏ + πΎ < πΏ π (1 β πΏ, βπΏ + π + 2, βπΎ < )π (βπΏ, βπΏ + π + 1, βπΎ < ) (A4) References: [1] J. T. Margaret, B. K. Chu, M. Roy, and E. L. Read, βDna-binding kinetics determines the mechanism of noise-induced switching in gene networks,β Biophysical journal , vol. 109, no. 8, pp. 1746β1757, 2015. [2] P. A. GonzΓ‘lez, L. J. CarreΓ±o, D. Coombs, J. E. Mora, E. Palmieri, B. Goldstein, S. G. Nathenson, and A. M. Kalergis, βT cell receptor binding kinetics required for t cell activation depend on the density of cognate ligand on the antigen-presenting cell,β
Proceedings of the National Academy of Sciences , vol. 102, no. 13, pp. 4824β4829, 2005. [3] R. Phillips, βNapoleon is in equilibrium,β
Annu. Rev. Condens. Matter Phys. , vol. 6, no. 1, pp. 85β111, 2015. 17 [4] A. B. Koudriavtsev, R. F. Jameson, and W. Linert,
The law of mass action . Springer Science & Business Media, 2011. [5] K. Ghosh, K. A. Dill, M. M. Inamdar, E. Seitaridou, and R. Phillips, βTeaching the principles of statistical dynamics,β
American journal of physics , vol. 74, no. 2, pp. 123β133, 2006. [6] R. G. Endres,
Physical principles in sensing and signaling: with an introduction to modeling in biology . Oxford University Press, 2013. [7] T. Mora, βPhysical limit to concentration sensing amid spurious ligands,β
Physical review letters , vol. 115, no. 3, p. 038102, 2015. [8] P. R. ten Wolde, N. B. Becker, T. E. Ouldridge, and A. Mugler, βFundamental limits to cellular sensing,β
Journal of Statistical Physics , vol. 162, no. 5, pp. 1395β1424, 2016. [9] J. V. Bonventre and V. S. Vaidya,
Biomarkers: In Medicine, Drug Discovery, and Environmental Health . Wiley, 2010. [10] H. G. Garcia, J. Kondev, N. Orme, J. A. Theriot, and R. Phillips, βA first exposure to statistical mechanics for life scientists,β arXiv preprint arXiv:0708.1899 , 2007. [11] F. Bai, X. Pi, P. Li, P. Zhou, H. Yang, X. Wang, M. Li, Z. Gao, and H. Jiang, βA statistical thermodynamic model for ligands interacting with ion channels: Theoretical model and experimental validation of the kcnq2 channel,β
Frontiers in pharmacology , vol. 9, p. 150, 2018. [12] D. A. Lauffenburger and J. J. Linderman,
Receptors: models for binding, trafficking, and signaling . Oxford University Press on Demand, 1996. [13] K. Ghosh, βStochastic dynamics of complexation reaction in the limit of small numbers,β
The Journal of chemical physics , vol. 134, no. 19, p. 05B606, 2011. [14] A. Cornish-Bowden,
Fundamentals of enzyme kinetics . John Wiley & Sons, 2013. [15] S. De Picciotto, B. Imperiali, L. G. Griffith, and K. D. Wittrup, βEquilibrium and dynamic design principles for binding molecules engineered for reagentless biosensors,β
Analytical biochemistry , vol. 460, pp. 9β15, 2014. [16] T. M. Squires, R. J. Messinger, and S. R. Manalis, βMaking it stick: convection, reaction and diffusion in surface-based biosensors,β
Nature biotechnology , vol. 26, no. 4, pp. 417β426, 2008. [17] B. Saha, T. H. Evers, and M. W. Prins, βHow antibody surface coverage on nanoparticles determines the activity and kinetics of antigen capturing for biosensing,β
Analytical chemistry , vol. 86, no. 16, pp. 8158β8166, 2014. [18] S. Iyer Biswas, βApplications of methods of non-equilibrium statistical physics to models of stochastic gene expression,β Ph.D. dissertation, The Ohio State University, 2009. [19] M. Abramowitz and I. Stegun, βHandbook of mathematical functions with formulas, graphs, and mathematical tables, june 1964,β
National Bureau of Standards Applied Mathematics Series
Advanced healthcare materials , vol. 8, no. 4, p. 1801478, 2019. [22] L. Cohen and D. R. Walt, βHighly sensitive and multiplexed protein measurements,β
Chemical reviews , vol. 119, no. 1, pp. 293β321, 2018. [23] M. A. G. de Castro, C. HΓΆbartner, and F. Opazo, βAptamers provide superior stainings of cellular receptors studied under super-resolution microscopy,β
PLoS One , vol. 12, no. 2, p. e0173050, 2017. [24] P. Waduge, R. Hu, P. Bandarkar, H. Yamazaki, B. Cressiot, Q. Zhao, P. C. Whitford, and M. Wanunu, βNanopore-based measurements of protein size, fluctuations, and conformational changes,β
ACS nano , vol. 11, no. 6, pp. 5706β5716, 2017. [25] M. Kaisti, βDetection principles of biological and chemical FET sensors,β