13.
Computer Lab – Chemical Kinetics 1
Goal: The purpose of this lab is to become familiar with modeling a simple ligand binding reaction.
Introduction
Consider a binding reaction where ligand ($L$) binds to a protein receptor ($R$) to form a loosely bound encounter complex ($L\cdot R$), and finally a tightly bound product ($P$) with the forward and backward rate constants $k_{f}$ and $k_{b}$, respectively.1
$$L + R \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} L\cdot R \stackrel{k_{f2}}{\rightarrow} P \tag{1}$$First, we need to write out the elementary rate equations:
\begin{align*} \\ \frac{d[L]}{dt} &= \frac{d[R]}{dt} = -k_{f1}[L][R] + k_{b1}[L\cdot R] \tag{1.1}\\ \frac{d[L\cdot R]}{dt} &= k_{f1}[L][R] -k_{b1}[L\cdot R] - k_{f2}[L\cdot R] \tag{1.2}\\ \frac{d[P]}{dt} &= k_{f2}[L\cdot R] \tag{1.3} \\ \end{align*}Next, we apply a steady state approximation on $d[L\cdot R]/dt$:
\begin{equation} k_{b1}[L\cdot R] = k_{f1}[L][R] - k_{f2}[L\cdot R] \implies \boxed{[L\cdot R] = \left( \frac{k_{f1}}{k_{b1}+k_{f2}} [L][R]\right)} \end{equation}Now, we take $[L\cdot R]$ and substitute into equations ${1.1}$ and ${1.3}$.
\begin{align} \frac{d[L]}{dt} &= \frac{d[R]}{dt} = -k_{f1}[L][R] + k_{b1}\left(\frac{k_{f1}}{k_{b1}+k_{f2}} [L][R] \right) \tag{1.1a}\\ \frac{d[P]}{dt} &= k_{f2} \left( \frac{k_{f1}}{k_{b1}+k_{f2}} [L][R]\right) \tag{1.3a} \end{align}Note that this last kinetic step is “one-way”, which violates thermodynamics (!). Don’t worry about this yet (we’re getting to kinetics soon in the lecture). For now, suffice it to say that here we invoke the assumption that the back rate is very small compared to forward rate.↩
Let’s walk through the steps of how to build a model of this receptor-ligand reaction using Python.
Procedure
Part I.
Set initial concentration of ligand to 10.0 and receptor to 5.0
Set
rate_eqs=ODE_eqs_part_Iset “Forward Rate Parameter” ($k_{f1}$) value to 0.02
set “Reverse Rate Parameter” ($k_{b1}$) value to 0.01
Run the simulation and save the graph.
Part II.
Repeat simulation for following values:
Add "k-prod" ($k_{f2}$) and set value to 0.05.
Set
rate_eqs=ODE_eqs_part_IIRun the simulation and save the graph.
Does the reaction reach a steady state?
If not, increase simulation time to 100
Repeat simulation for following values:
“Forward Rate Parameter” ($k_{f1}$): 0.01 , “Reverse Rate Parameter” ($k_{b1}$): 0.02 , “k-prod” ($k_{f2}$): 0.05
“Forward Rate Parameter” ($k_{f1}$): 0.02, “Reverse Rate Parameter” ($k_{b1}$): 0.01 , “k-prod” ($k_{f2}$): 0.5
import numpy as np
import pandas as pd
def ODE_eqs_part_I(C, t, k):
r"""
:math:`L + R \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} L\cdot R`
Args:
C(list): collection of concentrations for each species [L, R, LR]
t(float): iterate over some time-course
k(list): collection of rate constants [kf1, kb1]
"""
kf1, kb1 = k[0], k[1]
dC1dt = -kf1*C[0]*C[1]+kb1*C[2]
dC0dt = dC1dt
dC2dt = -dC1dt
return [dC0dt, dC1dt, dC2dt]
def ODE_eqs_part_II(C, t, k):
r"""See eqs 1.1, 1.2 and 1.3
:math:`L + R \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} L\cdot R \stackrel{k_{f2}}{\rightarrow} P`
Args:
C(list): collection of concentrations for each species [L, R, LR, P]
t(float): iterate over some time-course
k(list): collection of rate constants [kf1, kb1, kf2]
"""
kf1, kb1, kf2 = k[0], k[1], k[2]
dC1dt = -kf1*C[0]*C[1]+kb1*C[2]
dC0dt = dC1dt
dC2dt = k[0]*C[0]*C[1]-k[1]*C[2]-k[2]*C[2]
dC3dt = kf2*C[2]
return [dC0dt, dC1dt, dC2dt, dC3dt]
def ODE_eqs_part_II_SS(C, t, k):
r"""See eqs 1.2, 1.1a and 1.3a
Using the steady state approximation on dLR/dt.
:math:`L + R \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} L\cdot R \stackrel{k_{f2}}{\rightarrow} P`
Args:
C(list): collection of concentrations for each species [L, R, LR, P]
t(float): iterate over some time-course
k(list): collection of rate constants [kf1, kb1, kf2]
"""
kf1, kb1, kf2 = k[0], k[1], k[2]
dC1dt = -kf1*C[0]*C[1]+kb1*(kf1/(kb1+kf2)*C[0]*C[1])
dC0dt = dC1dt
dC2dt = k[0]*C[0]*C[1]-k[1]*C[2]-k[2]*C[2]
dC3dt = kf2*(kf1/(kb1+kf2)*C[0]*C[1])
return [dC0dt, dC1dt, dC2dt, dC3dt]
def simulate(steps, simulation_time, parameters, rate_eqs):
"""Method for integrating the ODEs with SciPy odeint
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
Args:
steps(int): number of iterations
simulation_time(int): duration of simulation (units: s)
parameters(dict): {'species': ['L', 'R', 'LR'],\
'k': [kf1, kb1, kf2, kb2, ...],\ # rate constants
'InitConc': [L, R, P, ...]} # initial concentrations (units: M)
rate_eqs(function): method containing the ordinary differential equations
Returns:
Pandas DataFrame, plot
"""
from scipy.integrate import odeint
print(f"Number of species: {len(parameters['InitConc'])}")
time = np.linspace(0, simulation_time, steps)
C = odeint(rate_eqs, parameters["InitConc"], time, (parameters["k"],)).transpose()
data = {"time": time, **{"%s"%s: C[i] for i,s in enumerate(parameters["species"])}}
df = pd.DataFrame(data)
plot = df.plot(x="time", y=parameters["species"], figsize=(10, 4),
xlim=(data["time"][0], data["time"][-1]))
return df,plot
parameters = {"species": ["L","R","LR"],
"InitConc": [10, 5, 0],
# [kf1, kb1, kf2]
"k": [0.02, 0.01, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=10,
parameters=parameters, rate_eqs=ODE_eqs_part_I)
plot = plot.set_title(r"$L + R ⥂ L\cdot R$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
parameters = {"species": ["L","R", "LR", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [0.02, 0.01, 0.05, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=100,
parameters=parameters, rate_eqs=ODE_eqs_part_II)
plot = plot.set_title(r"$L + R ⥂ L\cdot R \rightarrow P$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
parameters = {"species": ["L","R", "LR", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [0.02, 0.01, 0.05, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=100,
parameters=parameters, rate_eqs=ODE_eqs_part_II_SS)
plot = plot.set_title(r"""Steady State
$L + R ⥂ L\cdot R \rightarrow P$""")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.01 , “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.02 , “k-prod” ( 𝑘𝑓2 ): 0.05
parameters = {"species": ["L","R", "LR", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [0.01, 0.02, 0.05, "NaN"]}
df,plot = simulate(steps=1000, simulation_time=100,
parameters=parameters, rate_eqs=ODE_eqs_part_II)
plot = plot.set_title(r"$L + R ⥂ L\cdot R \rightarrow P$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.02, “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.01 , “k-prod” ( 𝑘𝑓2 ): 0.5
parameters = {"species": ["L","R", "LR", "P"],
"InitConc": [10, 5, 0, 0],
"k": [0.02, 0.01, 0.5, "NaN"]}
df,plot = simulate(steps=1000, simulation_time=100,
parameters=parameters, rate_eqs=ODE_eqs_part_II)
plot = plot.set_title(r"$L + R ⥂ L\cdot R \rightarrow P$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
Bonus
This same reaction can be used to describe Michaelis-Menten enzyme kinetics, in which substrate (S) binds a protein enzyme (E) to form an encounter complex (E∙S), and then chemically converts to product (E + P).
$$E + S \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} E\cdot S \stackrel{k_{cat}}{\to} E + P$$Suppose an enzyme has kf1 = 107 s-1, kb1 = 106 s-1 and kcat = 104 s-1. By adjusting parameters in the function above, which of these rates constants would you say most sensitively controls the overall rate of catalysis?
Answer: $k_{cat}$. See plots below...¶
parameters = {"species": ["E","S", "ES", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [1e7, 1e6, 1e4, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=1e-7,
parameters=parameters, rate_eqs=ODE_eqs_part_II)
plot = plot.set_title(r"$L + R ⥂ L\cdot R \rightarrow P$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
parameters = {"species": ["E","S", "ES", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [1e7, 1e6, 1e8, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=1e-7,
parameters=parameters, rate_eqs=ODE_eqs_part_II)
plot = plot.set_title(r"$L + R ⥂ L\cdot R \rightarrow P$")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
parameters = {"species": ["E","S", "ES", "P"],
"InitConc": [10, 5, 0, 0],
# [kf1, kb1, kf2, ...]
"k": [1e7, 1e6, 1e8, "NaN"]}
#print(pd.DataFrame(parameters))
df,plot = simulate(steps=1000, simulation_time=1e-7,
parameters=parameters, rate_eqs=ODE_eqs_part_II_SS)
plot = plot.set_title(r"""Steady State
$L + R ⥂ L\cdot R \rightarrow P$""")
fig = plot.get_figure()
#fig.savefig("figure.png")
df
