Protein binding kinetics

The aim of these posts is to explore protein and ligand binding kinetics and how they can be calculated.

We will use python to derive and visualization of steady-state concentration laws.

Monomer - Dimer kinetics

Lets start with the simple example of dimer formation from two identical protomers:

\[ \ce{P1 + P1 <=>[k_{on}][k_{off}] P2} \]

From the kinetic scheme above, we can write down the differential equations describing this system. These differential equations tell us how the concentrations of the reaction (in this case dimerization) change in time.

\[ \frac{\partial [P_1]}{\partial t} = - 2 k_{on} [P_1][P_1] + 2 k_{off} [P_2] \]

\[ \frac{\partial [P_2]}{\partial t} = k_{on} [P_1][P_1] - k_{off} [P_2] \]

Where square brackets are used to indicate the concentration of the species, ie \([P_1]\) is the concentration of the monomer at a given time point \(t\).

To sanity check our results, we can multiply the second equation with -2, such that both right-hand sides become equal. We then find that:

\[ \frac{\partial [P_1]}{\partial t} = - 2 \frac{\partial [P_2]}{\partial t} \]

When reading the equation above, for the equality to hold true, we find that the dimer \(P_2\) must appear at half the rate at which the monomer \(P_1\) disappears. Cross-checking this with our chemical intuition this makes sense, since two monomers are required to make one dimer.

In order to now find which concentrations of monomer and dimer are formed given how much protomer is in the tube, we need three more pieces of information.

First, we are looking for the concentration of the monomer and dimer at steady-state: the system has reached equillibrium and no more changes in the concentration of either the monomer or the dimer occur. This means we can set both differential equations to zero.

Second, we realize that the total amount of protomer in the tube never goes up or down; its either monomer or as one of the protomers in a dimer. The total amount of protomer is thus:

\[ [P_T] = [P_1] + 2[P_2] \]

Since extinction coefficients are additive (\(\epsilon_{P_2} = 2\epsilon_{P_1}\)) the total protomer concentration can be measured directly by UV-VIS spectroscopy.

Third, we need to have some value for the forward \(k_{on}\) and backward rates \(k_{off}\). These rates describe how fast the reaction takes place, but at the moment we are only interested in the steady-state concentrations. We can define a single new quantity which depends on the ratio of the forward and backward rates:

\[ K_d = \frac{k_{off}}{k_{on}} \]

If we look at our differential equations, the left-hand side describes the change of concentration of the monomer over time. The units are therefore molar per second (\(M s^{-1}\)). To make the units match on the left and right-hand side, we can deduce that the units of \(k_{on}\) must be \(M^{-1} s^{-1}\). Similarly, the units of \(k_{off}\) are \(s^{-1}\). From this unit analysis, we can see that the lifetime or off rate from a complex is independent of concentratation and we can recognize \(K_d\) as the dissociation constant with units \(M\).

We can now use sympy to do some mathematics and figure out what the steady-state concentrations are.

import sympy as sp

P1, P2, PT, k_on, k_off, kD = sp.symbols("P_1 P_2 P_T k_on k_off k_D", positive=True)

Now, we can solve the system of equations as we defined them above. Because of steady-state conditions, the differential equations are zero. We only use one here because the first is just the second multiplied by a constant factor, and thus does now have any additional information. To input equations into sympy, we must set the right hand side to zero and input the resulting left hand side:

sol = sp.solve(
    [
        -2 * k_on * P1 * P1 + 2 * k_off * P2,
        P1 + 2 * P2 - PT,
        (k_off / k_on) - kD,
    ],
    [P1, P2, k_on, k_off],
    dict=True,
)
sol

[{P_1: sqrt(k_D)*sqrt(8*P_T + k_D)/4 - k_D/4,
  P_2: P_T/2 - sqrt(k_D)*sqrt(8*P_T + k_D)/8 + k_D/8,
  k_off: k_D*k_on}]

This returns a single solution, because we have told sympy that all input symbols are positive. We can have a look at the solution for the concentration of the dimer:

sol[0][P2]

\(\displaystyle \frac{P_{T}}{2} - \frac{\sqrt{k_{D}} \sqrt{8 P_{T} + k_{D}}}{8} + \frac{k_{D}}{8}\)

We can see that indeed the solution only depends on the dissociation constant \(K_D\), and not the individual rates and indeed the ‘mass balance’ equation still holds (\(P_1 + 2P_2 = P_T\)), excercise for the reader.

Next, lets use solara to quickly make an interactive component so that we can easily calculate the monomer and dimer concentrations, given a total protomer concentration \(P_T\) and the dissociation constant.

First, we take the symbolic solutions from sympy and lambdify them so that we can input numbers and calculate the output:

solve_for = [P1, P2]
inputs = [PT, kD]

lambdas = {s: sp.lambdify(inputs, sol[0][s]) for s in solve_for}

Next, we make a solara component, where a user can input the values for \(P_T\) and \(K_{D}\) and directly obtain \(P_1\) and \(P_2\) at steady-state:

import solara

@solara.component
def Page():
    PT = solara.use_reactive(10.)
    kD = solara.use_reactive(1.)

    ans = {k: ld(PT.value, kD.value) for k, ld in lambdas.items()}

    solara.InputFloat('PT', value=PT)
    solara.InputFloat('kD', value=kD)

    for k, v in ans.items():
        solara.Text(f'{k.name}: {v}')

Page()

Note that the values are all unitless, so if we enter a dissociation constant of 1 \(\mu M\), the concentration \(P_T\) should also be given in \(\mu M\), and the outputs are then \(\mu M\) as well.

This little widget now tells us that if the monomer concentration if 10 times over the value of \(k_D\), in the steady-state equillibrium one out of three species is still a monomer!

Lets try if we can visualize the ratio of dimer to monomer over a larger range of concentrations. To do so we evalulate the monomer/dimer concentrations over a range spanning from 0.1 \(k_D\) to 1000 \(k_D\).

import numpy as np
import pandas as pd


# create a new function which calculates total molarity (monomer + dimer)
ld_total = sp.lambdify(inputs, sol[0][P1] + sol[0][P2])

# create a vector of PT values ranging from 0.1 times kD to 1000 times kD
PT_values = np.logspace(-1, 3, endpoint=True, num=100)
ans = {k: ld(PT_values, 1) / ld_total(PT_values, 1) for k, ld in lambdas.items()}

# put the results in a dataframe, together with input PT values
df = pd.DataFrame(dict(PT=PT_values) | {k.name: v for k, v in ans.items()})

df.head()

	PT	P_1	P_2
0	0.100000	0.921311	0.078689
1	0.109750	0.915248	0.084752
2	0.120450	0.908825	0.091175
3	0.132194	0.902035	0.097965
4	0.145083	0.894871	0.105129

The resulting values are fractions of monomer and dimer respectively, rather than their concentrations. We can plot the results with altair, adapting the gallery example for a multi-line tooltip graph.

import altair as alt

source = df.melt("PT", var_name="species", value_name="y")

# Create a selection that chooses the nearest point & selects based on x-value
nearest = alt.selection_point(nearest=True, on="pointerover",
                              fields=["PT"], empty=False)

# The basic line
line = alt.Chart(source).mark_line(interpolate="basis").encode(
    x=alt.X("PT:Q", scale=alt.Scale(type="log"), title='Ratio PT/kD'),
    y=alt.Y("y:Q", title='Fraction of total'),
    color="species:N",
)

# Draw points on the line, and highlight based on selection
points = line.mark_point().encode(
    opacity=alt.condition(nearest, alt.value(1), alt.value(0))
)

# Draw a rule at the location of the selection
rules = alt.Chart(source).transform_pivot(
    "species",
    value="y",
    groupby=["PT"]
).mark_rule(color="black").encode(
    x="PT:Q",
    opacity=alt.condition(nearest, alt.value(0.3), alt.value(0)),
    tooltip=[alt.Tooltip(c, type="quantitative", format=".2f") for c in df.columns],
).add_params(nearest)


# Put the five layers into a chart and bind the data
alt.layer(
    line, points, rules
).properties(
    width=600, height=300
)

From the graph, we can see that if we are at \(k_D\), the ratio monomer to dimer is two-to-one, and 50% dimerization is reached at 3 times \(k_D\), and for 95% saturation we need to be at > 500 times \(k_D\)!

So what about the ‘10x above \(k_D\) is full complex formation’ rule? This rule only applies some situations where we study the binding between different partners, such as the formation of an enzyme-ligand complex, and will be the subject of the next post (to be published).

Meanwhile, we can think about why its so hard to reach full dimerization even if concentrations used are far higher than the dissociation constant. If we think about the reaction from an maximum entropy point of view, out intuition might tell us that the entropic equivalent of making all promomers into a dimer is the equivalent of putting all ‘air molecules’ into one corner of the room: it has a vanishingly low probability of happing because its so far away from the maximium entropy state of the system.

Second, a closer look at the differential equation describing the change of \([P_2]\) also tells us the answer. There are two terms in this equation: one is positive and depends on \([P_1]\), the second is negative and depends on \([P_2]\). Therefore, even when there is a lot of total protomer \([P_T]\) compared to the \(k_D\), as more and more dimer is formed \([P_1]\) will go down while \([P_2]\) goes up. Thus the positive term becomes smaller while the negative term becomes bigger, slowing down the formation of dimer, and equillibrium is reached at a point with still a large fraction of \(P_1\) in solution.