first-order binding kinetics with multiple identical binding sites

Question

I have beads in solution at concentration $B$ which have $N$ independent binding sites for a single target molecule at concentration $T$. $T$ is small compared to $B$, so most beads will probably not have anything bound at all. The on and off rates for a single target binding a single site on the bead are $k_{on}$ and $k_{off}$.

I want to determine the equilibrium distribution of the number of bound targets per bead. I expect it to be a Poisson distribution (right?) with an average binding count that depends on the parameters above.

Define $B_n$ to be the concentration of beads with $n$ bound targets. First order kinetics suggest we should have the coupled system of equations:

$\frac{dB_n}{dt}=k_{on}T(B_{n-1}-B_n)+k_{off}(B_{n+1}-B_n)$,

$\frac{dT}{dt}=\sum_{n=0}^{N-1}\left(-k_{on}TB_n+k_{off}B_{n+1}\right)$

with $B_{-1}\equiv 0\equiv B_{N+1}$

First question: is there an analytical solution for this system of equations? I am perfectly happy to simulate it, but if there is a way to avoid that, it would be nice.

Second question: the wikipedia page for Dissociation Constant has some ambiguous wording suggesting that even if the receptors are all identical, the rate constants will be dependent on $n$ as well, via

$K_d(n) =K_d \frac{n}{N-n+1}$

where $K_d=\frac{k_{off}}{k_{on}}$

How does that play here? Are my equations in need of modification, to something like:

$\frac{1}{k_{on}}\frac{dB_n}{dt}=T(B_{n-1}-B_n)+K_d \frac{n}{N-n+1}(B_{n+1}-B_n)$.

If so, could someone explain where that modification comes from?

EDIT: it seems that this is probably due to the fact that there is a higher change of dissociation if a bead has many targets bound since there is an equal probability for each target to unbind. However, it seems that the fraction presented is 1 for $n=N$, which suggests that the dissociation constant is defined in terms of a fully bound bead and not on a per-site basis. Could someone show me the derivation?

Answer 1

The scheme you need seems to be of the form $$\ce{S + M_i <=>[k_2][k_{-1}] M_{i+1} }$$ where S is the molecule in solution and M the bead with i sites filled. If there is no restriction on how many S may be associated with any bead and as the S are indistinguishable molecules the number of ways of arranging the S among a total of M beads is $$ W=\frac{M!}{\prod M_i!}$$ and the constrains are $\Sigma_0^N M_i = M$ and $\Sigma iM_i = N$ where N is the total number of S molecules. Solving with these constraints is done using Lagrange multipliers, in the same manner as the Statistical Mechanics calculation of Boltzmann distribution, and leads to the distribution of finding i molecules of type S in a bead as, $$p_i = \frac{\bar n ^i}{(1+\bar n)^{i+1}}$$ where $\bar n$ is the average occupancy of a bead. This distribution, which is clearly not Poissonian, has the form that it decreases as i increases even for moderate $\bar n$.

If the kinetic scheme is changed slightly to $$\ce{S + M_i <=>[k_2][(i + 1)k_{-1}] M_{i+1} }$$ where the ‘off’ rate constant now depends on how many S are on a bead then the probability becomes Poissonian. This change/trick effectively makes the molecules distinguishable and now $$p_i = \frac{\bar n^i}{i!}e^{-\bar n}$$ with $\bar n = k_2\ce{[S]}/k_{-1}$.

Which of these or other models applies will have to be determined by experiment.

If you want to solve the equations numerically then it is best to try a Master equation approach.

The set of equations for the population of M is written as $$d\vec M/dt = \vec k \vec M$$ and has the formal solution $\vec M=\vec M_0e^{\vec kt}$ where $\vec M$ is a vector of concentrations of species i and $\vec k$ a matrix of rate constants such as

$$ \bar k=\left [\begin {matrix} -Sk_2 & 2k_{-1} & 0 & 0 & 0 & \cdots\\ Sk_2 & -Sk_2 -2k_{-1} & 3k_{-1} & 0 & 0 & \cdots \\ \vdots& \vdots&\vdots & \vdots& \vdots& \\ \end {matrix}\, \right ]$$

where S is the concentration of S and is assumed to be in excess so that its value is constant. (If you want S to vary then you will need to do a numerical integration on all N equations instead of this method).

To solve, the secular determinant is set up and the eigenvalues $\lambda$ obtained. The populations at time t are given by $$M(t)= x[e^{\lambda t}]x^{-1}M_0$$ where x is the matrix of eigenvectors corresponding to eigenvalues $\lambda$ and $[e^{\lambda t}] $ is a diagonal matrix of $e^{\lambda_1t}, e^{\lambda_2 t} \dots $ etc.

Answer 2

I was able to figure out the answer to my second question using the information provided by @porphyrin in the accepted answer, but I am posting my own solution here for the sake of completeness. Specifically, the reasoning behind the effective rate constants changing when there are multiple bound targets per bead discussed cryptically on the wikipedia page here.

Let $k_{on}^{(0)}$ and $k_{off}^{(1)}$ be the rate constants for binding and unbinging of a target molecule to a single site on the bead. Because we are assuming linear kinetics, if there are multiple binding sites free on the bead, the rate of binding to that bead will simply be multiplied by the number of free binding sites. Similarly, if there are multiple bound targets, the unbinding rate for that bead will be scaled accordingly.

Consider the reaction $$\ce{T + B_{n-1}<=>[k_{on}^{(n-1)}][k_{off}^{(n)}]B_n}$$

Because there are $(N-(n-1))$ free sites on the bead on the left hand side, the effective on rate is $$k_{on}^{(n-1)}=k_{on}^{(0)}(N-(n-1))$$

On the right hand side, because there are $n$ bound targets, the effective off rate for this bead is $$k_{off}^{(n)}=nk_{off}^{(1)}$$.

Therefore, the dissociation constant for this particular reaction is $$K_d^{(n)}=\frac{k_{off}^{(n)}}{k_{on}^{(n-1)}}=\frac{k_{off}^{(1)}}{k_{on}^{(0)}}\frac{n}{N-n+1}$$ as required.