<h4><span style='color:blue'>13.</span></h4>


<p><strong>Computer Lab – Chemical Kinetics 1</strong></p>
<p><strong>Goal</strong>: The purpose of this lab is to become familiar with modeling a simple ligand binding reaction.</p>
<p><strong>Introduction</strong></p>
<p>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.<a href="#fn_note1" class="footnote-ref" id="fn_note1"><sup>1</sup></a></p>

$$L + R \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} L\cdot R \stackrel{k_{f2}}{\rightarrow} P \tag{1}$$


<p>First, we need to write out the elementary rate equations:</p>

\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}


<section class="footnotes">
<hr>
<ol>
<li id="fn_note1"><p>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.<a href="#fn_note1" class="footnote-back">↩</a></p></li>
</ol>
<hr>
</section>




<p>Let’s walk through the steps of how to build a model of this receptor-ligand reaction using Python.</p>

<p><strong>Procedure</strong></p>
<p><strong>Part I.</strong></p>
<ol class="incremental" type="1">
<li><p>Set initial concentration of ligand to 10.0 and receptor to 5.0</p></li>
<li><p>Set <code>rate_eqs=ODE_eqs_part_I</code></p></li>
<ol class="incremental" type="a">
<ol class="incremental" type="i">
<li><p>set “Forward Rate Parameter” ($k_{f1}$) value to 0.02</p></li>
<li><p>set “Reverse Rate Parameter” ($k_{b1}$) value to 0.01</p></li>
</ol></li>
</ol></li>
<li><p>Run the simulation and save the graph.</p></li>
</ol></li>
</ol>
<p><strong>Part II.</strong></p>
<ol class="incremental" type="1">
<li><p>Repeat simulation for following values:</p>
<li><p>Add "k-prod" ($k_{f2}$) and set value to 0.05.</p></li>
<li><p>Set <code>rate_eqs=ODE_eqs_part_II</code></p></li>
<li><p>Run the simulation and save the graph.</p></li>
<ol class="incremental" start="1" type="a">
<li><p>Does the reaction reach a steady state?</p></li>
<li><p>If not, increase simulation time to 100 </p></li>
<li><p>Set <code>rate_eqs=ODE_eqs_part_II_SS</code></p></li>   
<li><p>Run the simulation. How does this plot/simulation differ from the last?</p></li>
</ol></li>
<li><p>Repeat simulation for following values:</p>
<ol class="incremental" start="5" type="a">
<li><p>Set <code>rate_eqs=ODE_eqs_part_II</code></p></li>  
<li><p>“Forward Rate Parameter” ($k_{f1}$): 0.01 , “Reverse Rate Parameter” ($k_{b1}$): 0.02 , “k-prod” ($k_{f2}$): 0.05</p></li>
<li><p>“Forward Rate Parameter” ($k_{f1}$): 0.02, “Reverse Rate Parameter” ($k_{b1}$): 0.01 , “k-prod” ($k_{f2}$): 0.5</p></li>
</ol></li>
<!--
<li><p>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:</p>
<ol class="incremental" start="3" type="a">    
<li><p>Set <code>rate_eqs=ODE_eqs_part_II_SS</code></p></li>
<li><p>What differences do you see?</p></li>
</ol></li>
-->
</ol>
    
    


```python
import numpy as np
import pandas as pd
```


```python
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


```


```python
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
```


```python
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
```


```python
#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.01 , “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.02 , “k-prod” ( 𝑘𝑓2 ): 0.05

```


```python
#“Forward Rate Parameter” ( 𝑘𝑓1 ): 0.02, “Reverse Rate Parameter” ( 𝑘𝑏1 ): 0.01 , “k-prod” ( 𝑘𝑓2 ): 0.5

```

<p><strong>Bonus</strong></p>
<p>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).</p>


$$E + S \underset{k_{b1}}{\stackrel{k_{f1}}{⥂}} E\cdot S \stackrel{k_{cat}}{\to} E + P$$




<p>Suppose an enzyme has <em>k</em><sub>f1</sub> = 10<sup>7</sup> s<sup>-1</sup>, <em>k</em><sub>b1</sub> = 10<sup>6</sup> s<sup>-1</sup> and <em>k</em><sub>cat</sub> = 10<sup>4</sup> s<sup>-1</sup>. By adjusting parameters in the function above, which of these rates constants would you say most sensitively controls the overall rate of catalysis?</p>




```python

```

b8dfc4590d8703fa25aff0e4b24aa8b42209de3a
ipynb
Jupyter Notebook
CHEM3405_Physical_Chemistry_Bio/Labs/CHEM3405-Chemical_Kinetics_I.ipynb
robraddi/tu_chem
18b8247d6c00e33f15f040a57a32b5fc2372137a
MIT
2021-04-29T04:26:42.000Z
2021-04-29T04:26:42.000Z
CHEM3405_Physical_Chemistry_Bio/Labs/CHEM3405-Chemical_Kinetics_I.ipynb
robraddi/tu_chem
18b8247d6c00e33f15f040a57a32b5fc2372137a
MIT
CHEM3405_Physical_Chemistry_Bio/Labs/CHEM3405-Chemical_Kinetics_I.ipynb
robraddi/tu_chem
18b8247d6c00e33f15f040a57a32b5fc2372137a
MIT
2020-12-03T17:47:05.000Z
2020-12-03T17:47:05.000Z
Qwen/Qwen-72B
1. YES
2. YES