<jupyter_start><jupyter_text>13.Computer Lab – Chemical Kinetics 1Goal: The purpose of this lab is to become familiar with modeling a simple ligand binding reaction.IntroductionConsider 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.ProcedurePart I.Set initial concentration of ligand to 10.0 and receptor to 5.0Set rate_eqs=ODE_eqs_part_Iset “Forward Rate Parameter” ($k_{f1}$) value to 0.02set “Reverse Rate Parameter” ($k_{b1}$) value to 0.01Run 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 Set rate_eqs=ODE_eqs_part_II_SS   Run the simulation. How does this plot/simulation differ from the last?Repeat simulation for following values:Set rate_eqs=ODE_eqs_part_II  “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<!--Now, add a new cell block under each simulation in Part II.edit the rate equations using the steady state approximation and compare the simulations:    Set rate_eqs=ODE_eqs_part_II_SSWhat differences do you see?--><jupyter_code>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
#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.01 , “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.02 , “k-prod” ( 𝑘𝑓2 ): 0.05

#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.02, “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.01 , “k-prod” ( 𝑘𝑓2 ): 0.5<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1782106782106782, 'fraction_numerical': 0.03939393939393939, 'mean_word_length': 3.4572347266881027, 'pattern_counts': {'":': 7, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '268562', 'license': 'MIT', 'n_tokens_mistral': 2844, 'n_tokens_neox': 2687, 'n_words': 706, 'path': 'CHEM3405_Physical_Chemistry_Bio/Labs/CHEM3405-Chemical_Kinetics_I.ipynb', 'repo_name': 'robraddi/tu_chem', 'size': '6930'}