Simulating data

Simulating data#

Tools#

During these exercises, we will learn some general Python and biochemistry tools:

$A$ $\to$ $B$
- applying numerical integration methods to calculate the change in concentration of $A$ and $B$ with time
$A + B$ $⇌$ $A B$
- applying numerical integration methods to calculate the change in concentration of $A$ , $B$ , and $A B$ with time
$E + S$ $⇌$ $E S$ $\to$ $E + P$
- applying numerical integration methods to calculate the change in concentration of $E$ , $S$ , $E S$ , and $P$ with time

Why?#

This taster session applies numerical integration methods to ordinary differential equations to calculate the change in concentrations for simple biochemical reactions. The same mathematical modeling can be used for biochemical pathways and networks to help us understand how e.g. environmental stimuli induce various phenotypes via sequences of biochemical reactions and time-varying concentrations of biomolecular species.

Interested? Have a look at: Aldridge, B., Burke, J., Lauffenburger, D. et al. Physicochemical modelling of cell signalling pathways. Nat Cell Biol 8, 1195–1203 (2006). https://doi.org/10.1038/ncb1497

Convention#

We will use the convention that $[A]$ is the concentration of free $A$ , $[A]_{0}$ is the initial (or total) concentration of $A$ before any binding or conversion occurs, and $[A]_{e q}$ is the concentration of $A$ at equilibrium.

Import libraries#

Before we start, we need to import the correct libraries, i.e. NumPy, Matplotlib, pandas, and the odeint and curve_fit functions from SciPy.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.integrate import odeint 
from scipy.optimize import curve_fit

$A$ $\to$ $B$ #

Background#

In the following, simple mechanism, $A$ is the reactant and $B$ is the product:

A \overset{k_{f}}{\to} B .

This irreversible first-order reaction is a reaction that proceeds at a reaction rate, $v$ $(M s^{- 1})$ , which depends linearly on only one reactant concentration. Here, the rate of disappearance of $A$ is proportional to the concentration of $A$ itself. The forward rate constant, $k_{f}$ $(s^{- 1})$ , is the proportionality constant that relates the observed reaction rate to concentration. This can be written as:

v = - \frac{d [A]}{d t} = k_{f} [A] .

This is an ordinary differential equation (ODE). We want to find how the concentration of $A$ depends on the time $t$ (independent variable). We can solve this ODE analytically or numerically.

Solving this equation analytically shows an exponential decay defined by

\int_{[A]_{0}}^{[A]} \frac{d [A]}{[A]} = \int_{0}^{t} - k_{f} d t .

l n [A] - l n [A]_{0} = - k_{f} t .

[A] = [A]_{0} e^{- k_{f} t} .

and because $[B]$ = $[A]_{0}$ - $[A]$ :

[B] = [A]_{0} (1 - e^{- k_{f} t}) .

The rate constant can be estimated from the relationship

k_{f} = \frac{l n 2}{t_{1 / 2}},

where the half life or half time, $t_{1 / 2}$ , is the time for $[A]_{0}$ to decrease by 50 %. In molecular terms $t_{1 / 2}$ represents the time in which there is a 50 % probability that a molecule $A$ will be converted to $B$ . Those molecules which have not converted in this time will still have a 50 % probability of converting - so the population will fall from 100 % to 50 % to 25 % to 12.5 % to 6.25 % to … for each half time period.

To solve this ODE numerically, we simply take a small but finite time interval $Δ$ $t = t_{2} - t_{1}$ and calculate the finite change in $[A]$ . In Python, we can do this using the odeint function of the SciPy package. Briefly, odeint requires four main parameters: odeint(model, y0, t, args=), with

model: a function that returns derivative values at requested y and t values: dydt = model(y, t)
y0: the initial conditions of the y states
t: the time points for which the solutions need to be reported
args: when the model function depends on variables other than y and t, these can be passed in a tuple. Of note, a tuple consists of a number of values separated by commas. It is different from a list as a tuple is a collection which is ordered and unchangeable. For more information, see here.

See also the odeint documentation, available here.

Simulation#

We start be defining a function that returns $\frac{d [A]}{d t}$ and $\frac{d [B]}{d t}$ . The input for our function will be a list with initial concentrations ( $C = {[A_{0}], [B_{0}]}$ , all in $m M$ ), time points for which values of the variables should be returned ( $t$ in $s$ ), and the rate constant ( $k_{f}$ in $s^{- 1}$ ).

def funcAtoB(C, t, kf):   #create the function
    """
    Define a function for A -> B 

    Args:
        C, the concentrations of A and B (float) in mol/L
        t, the time (float) in s
        kf, the forward rate constant (float) in 1/s

    Returns:
        dC/dt, derivative values (float) in (mol/L)/s
    """
    CA = C[0]   #the initial concentration of A is the first element from the list C
    CB = C[1]   #the initial concentration of B is the second element from the list C
    
    dAdt = -kf * CA   #rate equation for [A]
    dBdt = +kf * CA   #rate equation for [B]
    
    return [dAdt, dBdt]

We then provide the initial concentrations, time points, and a rate constant. We use the following parameters:

initial concentration of $A$ : $[A]_{0} = 100$ $m M$
initial concentration of $B$ : $[B]_{0} = 0$ $m M$
rate constant: $k_{f} = 7$ $s^{- 1}$
time points: 1000 points between $0$ and $2 s$ , generated using numpy.linspace(start, stop, num).

CAtoB = [100, 0]   #y0 - initial conditions

tAtoB = np.linspace(0, 2, 1000)   #t - time points, an array with 1000 evenly distributed elements between 0 (included) and 2 (included)

kfAtoB = (7,)   #args - rate constant. Of note, a tuple with one item is constructed by following a value with a comma!

We can now solve the ODE for each time point:

yAtoB = odeint(funcAtoB, CAtoB, tAtoB, kfAtoB)   #solve ODE - odeint(model, y0, t, args=)
print(yAtoB)   #print the outcome

[[1.00000000e+02 0.00000000e+00]
 [9.86083725e+01 1.39162750e+00]
 [9.72361113e+01 2.76388871e+00]
 ...
 [8.55163449e-05 9.99999145e+01]
 [8.43263420e-05 9.99999157e+01]
 [8.31528952e-05 9.99999168e+01]]

Let’s use matplotlib.pyplot.plot to create the $[A]$ and $[B]$ versus $t$ plot.

plt.figure(figsize=(7,5))   #to create a figure object
plt.plot(tAtoB,yAtoB[:,0],'g')   #select the first column of yAtoB, i.e. [A], and plot versus t, use green
plt.plot(tAtoB,yAtoB[:,1],'b')   #select the second column of yAtoB, i.e. [B], and plot versus t, use blue
plt.xlabel('Time (s)', fontsize=18)   #label the X-axis
plt.ylabel('[A] (mM)', fontsize=18)   #label the Y-axis
plt.xlim([0,1])   #set X-axis range
plt.ylim([0,110])   #set Y-axis range
plt.legend(['A','B'], loc='best')   #include a legend
plt.show()   #show the figure object

../_images/997a115f6e8fb803b5888e885092d2d7d83b5d98176bf05fcb201f8a505c4c96.png

Nice! The curve shows the signal we can expect (whether we measure the dissappearance of A or the appearance of B) for an irreversible first-order reaction $A$ $\to$ $B$ , with $k_{f} = 7$ $s^{- 1}$ and $[A]_{0} = 100$ $m M$ .

Exercise 80

How would you determine $k_{f}$ in the lab?

Solution to Exercise 80

If our reactant A absorbs light at a given wavelength (and our product B does not), then we can follow the absorbance and use the Beer-Lambert law to determine the concentration of reactant A, $[A]$ , as a function of time.

A plot of $[A]$ as a function of time will decay exponentially from $[A]_{0}$ down to zero. We can fit our experimental data to this function

[A] = [A]_{0} e^{- k_{f} t}

to determine $k_{f}$ .

When we take the natural logarithm of both sides, we obtain:

l n ([A]) = l n ([A]_{0}) - k_{f} t .

For an irreversible first-order reaction, a plot of $l n ([A])$ as a function of time will yield a straight line and the slope will be equal to $k_{f}$ .

This is demonstrated in the following figure:

Determining the forward rate constant experimentally for irreversible first-order reaction

$A + B$ $⇌$ $A B$ #

Background#

We now use the numerical approach to simulate the reaction profile for a reversible second-order binding reaction:

A + B \underset{k_{r}}{\overset{k_{f}}{⇌}} A B .

Where $k_{f}$ and $k_{r}$ are rate constants for the forward and reverse reactions between $A$ and $B$ . Bear in mind that $k_{f}$ is now a second order rate constant with units $M^{- 1} s^{- 1}$ since the rate of binding is dependent on both $[A]$ and $[B]$ . On the other hand, the reverse rate constant $k_{r}$ is first order ( $s^{- 1}$ ) since it is only dependent on one component, $[A B]$ .

The following differential equations, describing the rate of variation of $[A]$ , $[B]$ and $[A B]$ , can be derived from the binding reaction:

\begin{array}{r} \begin{matrix} \frac{d [A]}{d t} = - k_{f} [A] [B] + k_{r} [A B], \\ \frac{d [B]}{d t} = - k_{f} [A] [B] + k_{r} [A B], \\ \frac{d [A B]}{d t} = k_{f} [A] [B] - k_{r} [A B] . \end{matrix} \end{array}

Note that we need to include a term for the reformation of $[A]$ from $[A B]$ as the reaction is now reversible. It will make little difference to the early stages of the reaction when $[A B] = 0$ or is small, but it will determine the position of the equilibrium at longer times. Indeed, when the reaction is “complete”, and there is no change in $[A]$ , $[B]$ , and $[A B]$ , the reaction has reached its constant equilibrium value. Be aware that there is still a flux of molecules going forwards and backwards between the $A$ and $A B$ states! This can be defined by:

k_{f} [A]_{e q} [B]_{e q} = k_{r} [A B]_{e q} .

Rearranging gives the equilibrium dissociation constant $K_{d}$ :

K_{d} = \frac{[A]_{e q} [B]_{e q}}{[A B]_{e q}} = \frac{k_{r}}{k_{f}} .

There are analytical solutions to $\frac{d [A B]}{d t}$ when component $B$ is in large excess compared to $A$ to give pseudo first order conditions ( $[B]_{0} ≫ [A]_{0}$ ). Under these conditions, the reaction is simplified because $[B]$ stays practically constant and $\sim [B]_{0}$ . Therefore, the profile appears as a single exponential ( $[A B] = [A B]_{e q} (1 - e^{- k_{o b s} t})$ ) with an apparent rate constant $k_{o b s} = k_{f} [B]_{0} + k_{r}$ . However, if $[B]_{0}$ is not much larger than $[A]_{0}$ this approximation gradually breaks down.

Exercise 81

Prove that $[A B] = [A B]_{e q} (1 - e^{- k_{o b s} t})$ with an apparent rate constant $k_{o b s} = k_{f} [B]_{0} + k_{r}$ for pseudo first order conditions ( $[B]_{0} ≫ [A]_{0}$ ).

Solution to Exercise 81

One can also solve this set of ODEs numerically using the odeint function of the SciPy package.

Simulation#

We start be defining a function that returns $\frac{d [A]}{d t}$ , $\frac{d [B]}{d t}$ , and $\frac{d [A B]}{d t}$ . The input for our function will be a list with initial concentrations ( $C = {[A], [B], [A B]}$ , all in $μ$ $M$ ), time $t$ in $s$ , and the rate constants $k_{f}$ (in $μ$ $M$ $s^{- 1}$ ), and $k_{r}$ (in $s^{- 1}$ ).

def funcAandBtofromAB(C, t, kf, kr):   #create the function
    """
    Define a function for A + B -> AB 

    Args:
        C, the concentrations of A, B, and AB (float) in mol/L
        t, the time (float) in s
        kf, the forward rate constant (float) in (mol/L)/s
        kr, the reverse rate constant (float) in 1/s

    Returns:
        dC/dt, derivative values (float) in (mol/L)/s
    """
    CA = C[0] #the initial concentration of A is the first element from the list C
    CB = C[1] #the initial concentration of B is the second element from the list C
    CAB = C[2] #the initial concentration of AB is the thirs element from the list C
    
    dAdt = -kf * CA * CB + kr * CAB #rate equation for [A]
    dBdt = -kf * CA * CB + kr * CAB #rate equation for [B]
    dABdt = kf * CA * CB - kr * CAB #rate equation for [AB]
    
    return [dAdt, dBdt, dABdt]

We then provide the initial concentrations, time points, and a rate constant. We use the following parameters:

initial concentration of $A$ : $[A]_{0} = 10$ $μ$ $M$
initial concentration of $B$ : $[B]_{0} = 50$ $μ$ $M$
initial concentration of $A B$ : $[A B]_{0} = 0$ $μ$ $M$
rate constants: $k_{f} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{r} = 10$ $s^{- 1}$ .
time points: 1000 points between $0$ and $0.1 s$

CAandBtofromAB = [10, 50, 0]   #y0 - initial conditions

tAandBtofromAB = np.linspace(0, 0.1, 1000)   #t - time points, an array with 1000 evenly distributed elements between 0 (included) and 0.1 (included)

ksAandBtofromAB = (1, 10)   #args - rate constants

We can now solve the ODE for each time point:

yAandBtofromAB = odeint(funcAandBtofromAB, CAandBtofromAB, tAandBtofromAB, ksAandBtofromAB)   #solve ODE - odeint(model, y0, t, args=)
print(yAandBtofromAB)   #print the outcome

[[1.00000000e+01 5.00000000e+01 0.00000000e+00]
 [9.95012480e+00 4.99501248e+01 4.98752049e-02]
 [9.90059739e+00 4.99005974e+01 9.94026072e-02]
 ...
 [1.95837942e+00 4.19583794e+01 8.04162058e+00]
 [1.95820429e+00 4.19582043e+01 8.04179571e+00]
 [1.95803011e+00 4.19580301e+01 8.04196989e+00]]

We can now plot our traces.

plt.figure(figsize=(7,5))  #to create a figure object
plt.plot(tAandBtofromAB,yAandBtofromAB[:,0],'g')   #select the first column of yAandBtofromAB, i.e. [A], and plot versus t, use green
plt.plot(tAandBtofromAB,yAandBtofromAB[:,1],'b')   #select the second column of yAandBtofromAB, i.e. [B], and plot versus t, use blue
plt.plot(tAandBtofromAB,yAandBtofromAB[:,2],'r')   #select the third column of yAandBtofromAB, i.e. [AB], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18)   #label the X-axis
plt.ylabel('[x] ($\mu$M)', fontsize=18)   #label the Y-axis
plt.xlim([0,0.1])   #set X-axis range
plt.ylim([0,55])   #set Y-axis range
plt.legend(['A','B', 'AB'], loc='best')   #include a legend
plt.show()   #show the figure object

<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:6: SyntaxWarning: invalid escape sequence '\m'
C:\Users\ucbtrv2\AppData\Local\Temp\ipykernel_26512\467036247.py:6: SyntaxWarning: invalid escape sequence '\m'
  plt.ylabel('[x] ($\mu$M)', fontsize=18)   #label the Y-axis

../_images/7090be6f9a2a18315c735878b593af3402ec48fbccf3508721b35149fa79dba4.png

Nice! The curve shows the signal we can expect (whether we measure the dissappearances of A or B or the appearance of AB) for a reversible second-order reaction $A + B$ $⇌$ $A B$ , with $k_{f} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{r} = 10$ $s^{- 1}$ ; $[A]_{0} = 10$ $μ$ $M$ ; and $[B]_{0} = 50$ $μ$ $M$ .

Exercise 82

How would you determine $k_{f}$ and $k_{r}$ in the lab?

Solution to Exercise 82

If our product AB absorbs light at a given wavelength (and our reactants A and B do not), then we can follow the absorbance and use the Beer-Lambert law to determine the concentration of product AB, $[A B]$ , as a function of time after mixing A and B.

When we use “pseudo first order conditions”, $[B]_{0} ≫ [A]_{0}$ , the time course to form AB is an exponential with observed rate constant $k_{o b s}$ that depends on $[B]_{0}$ . Therefore, by varying $[B]_{0}$ , one can evaluate $k_{f}$ and $k_{r}$ from the slope and intercept of a plot of $k_{o b s}$ versus $[B]_{0}$ .

More information about pre-steady state kinetics to determine rate constants can be found in

Pollard TD, De La Cruz EM. Take advantage of time in your experiments: a guide to simple, informative kinetics assays. Mol Biol Cell. 2013;24(8):1103-1110. doi:10.1091/mbc.E13-01-0030
Johnson KA. Transient state kinetic analysis of enzyme reaction pathways. Enzymes. 1992;20:1-61.doi:10.1016/S1874-6047(08)60019-0

$E + S$ $⇌$ $E S$ $\to$ $E + P$ #

Background#

In the following, simple mechanism, $E$ is the enzyme, $S$ is the substrate, $E S$ is the enzyme-substrate complex, and $P$ is the product:

E + S \underset{k_{- 1}}{\overset{k_{+ 1}}{⇌}} E S \overset{k_{2}}{\to} E + P .

Where $k_{+ 1}$ and $k_{- 1}$ are rate constants for the forward and reverse reactions between $E$ and $S$ , respectively, and $k_{2}$ , the rate constant for the $E S \to E + P$ reaction.

We will now simulate this enzyme-catalysed reaction using a numerical approach to see how $[E]$ , $[S]$ , $[E S]$ , and $[P]$ (in $μ$ $M$ ) change as a function of time. These type of curves are called progress curves. This is usually what you measure in the lab when you study enzyme reactions!

We use the following parameters:

initial concentration of enzyme, $[E]_{0} = 1$ $μ$ $M$
initial concentration of substrate, $[S]_{0} = 10$ $μ$ $M$
initial concentration of enzyme-substrate complex, $[E S]_{0} = 0$ $μ$ $M$
initial concentration of product, $[P]_{0} = 0$ $μ$ $M$
rate constants: $k_{+ 1} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{- 1} = 0.1$ $s^{- 1}$ ; $k_{2} = 0.5$ $s^{- 1}$
time points: 10000 points between $0$ and $60 s$

Simulation#

Exercise 83

It is up to you to

Derive the differential equations, describing the rate of variation of $[E]$ , $[S]$ , $[E S]$ , and $[P]$ , from the reaction.
Define a function that returns $\frac{d [E]}{d t}$ , $\frac{d [S]}{d t}$ , $\frac{d [E S]}{d t}$ , and $\frac{d [P]}{d t}$ .
Provide rate constants, initial concentrations, and time points, and evaluate the function for each time point.
Plot the resulting time traces from 0 to 60 $s$ .
Discuss the simulated data.

Solution to Exercise 83

Derive the differential equations, describing the rate of variation of $[E]$ , $[S]$ , $[E S]$ , and $[P]$ , from the reaction.

\begin{array}{r} \begin{matrix} \frac{d [E]}{d t} = - k_{+ 1} [E] [S] + (k_{- 1} + k_{2}) [E S], \\ \frac{d [S]}{d t} = - k_{+ 1} [E] [S] + k_{- 1} [E S], \\ \frac{d [E S]}{d t} = k_{+ 1} [E] [S] - (k_{- 1} + k_{2}) [E S], \\ \frac{d [P]}{d t} = k_{2} [E S] . \end{matrix} \end{array}

Define a function that returns $\frac{d [E]}{d t}$ , $\frac{d [S]}{d t}$ , $\frac{d [E S]}{d t}$ , and $\frac{d [P]}{d t}$ . Here, I used $k_{+ 1} = k 1$ ; $k_{- 1} = k 2$ ; $k_{2} = k 3$ .

#model - function that returns dE/dt, dS/dt, dES/dt, and dP/dt
def funcEandStofromEStoP(C, t, k1, k2, k3) :
    # constants
    CE = C[0]
    CS = C[1]
    CES = C[2]
    CP = C[3]
    
    # equations
    dEdt = -k1 * CE * CS + (k2 + k3) * CES
    dSdt = -k1 * CE * CS + k2 * CES
    dESdt = k1 * CE * CS - (k2 + k3) * CES
    dPdt = k3 * CES
    
    return [dEdt, dSdt, dESdt, dPdt]

Provide rate constants, initial concentrations, and time points, and evaluate the function for each time point.

#y0 - initial conditions
CEandStofromEStoP = [1, 10, 0, 0]

#t - time points
tEandStofromEStoP = np.linspace(0, 60, 10000)

#args - rate constants
ksEandStofromEStoP = (1, 0.1, 0.5)

#solve ODE - odeint(model, y0, t, args=)
yEandStofromEStoP = odeint(funcEandStofromEStoP, CEandStofromEStoP, tEandStofromEStoP, ksEandStofromEStoP)
print(yEandStofromEStoP)

Plot the resulting time traces from 0 to 60 $s$ .

#create plot
plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,0],'g:') #select the first column of yEandStofromEStoP, i.e. [E], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,1],'b-') #select the second column of yEandStofromEStoP, i.e. [S], and plot versus t, use blue
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,3],'y-.') #select the fourth column of yEandStofromEStoP, i.e. [P], and plot versus t, use yellow
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[x] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,60]) #set X-axis range
plt.ylim([0,11]) #set Y-axis range
plt.legend(['E','S', 'ES', 'P'], loc='best') #include a legend
plt.show()

Discuss the simulated data.

Note that the reaction velocity or rate (i.e., the slope of the $[P]$ versus $t$ curve ( $= \frac{d [P]}{d t} = - \frac{d [S]}{d t}$ ) decreases in time due to e.g. substrate depletion. Only early in the progress curve, product formation is linear with time. If we determine the initial reaction rate ( $v_{0}$ ) from the slope of the curve during this early phase, we circumvent substrate depletion.

Note the three phases:

pre-steady state, with a burst in $[E S]$
steady state, with almost constant $[E S]$ , and
end, where $[E S]$ drops. Let us analyse these three phases in detail in the following exercise.

Pre-steady state and steady state phases#

Exercise 84

Let’s analyse $[E]$ and $[E S]$ versus $t$ in more detail. When we analyze the data from the previous exercise, we can distinguish three different phases. Change the X- and Y-axis formats of the plot from the previous exercise to zoom in on these three phases:

The pre-steady state (or transient) phase, in which there is a burst in $[E S]$ . In this case between 0 and ~ 0.5 $s$ .
The steady state phase, in which $[E S]$ remains almost constant. Here, between 0.5 and ~ 10 $s$ .
After ~ 10 $s$ , $[E S]$ drops and there is a progressively nonlinear formation of product.

Solution to Exercise 84

The pre-steady state (or transient) phase, in which there is a burst in $[E S]$ . In this case between 0 and ~ 0.5 $s$ .

#create plot for [E] and [S] from 0 to 0.5 s
plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,0],'g:') #select the first column of yEandStofromEStoP, i.e. [E], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[x] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,0.5]) #set X-axis range
plt.ylim([0,1.1]) #set Y-axis range
plt.legend(['E','ES'], loc='best') #include a legend
plt.show()

The steady state phase, in which $[E S]$ remains almost constant. Here, between 0.5 and ~ 10 $s$ .

#create plot for [E] and [S] from 0.5 to 10 s
plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,0],'g:') #select the first column of yEandStofromEStoP, i.e. [E], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[x] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0.5,10]) #set X-axis range
plt.ylim([0,1.1]) #set Y-axis range
plt.legend(['E', 'ES'], loc='best') #include a legend
plt.show()

After ~ 10 $s$ , $[E S]$ drops and there is a progressively nonlinear formation of product.

#create plot for [E] and [S] from 0 to 60 s
plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,0],'g:') #select the first column of yEandStofromEStoP, i.e. [E], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[x] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,60]) #set X-axis range
plt.ylim([0,1.1]) #set Y-axis range
plt.legend(['E','ES'], loc='best') #include a legend
plt.show()

Experimental design: $[S]_{0}$ relative to $[E]_{0}$ #

Exercise 85

Let’s further investigate what effect changing the value of $[S]_{0}$ relative to $[E]_{0}$ has for the same, simple reaction mechanism.

Compare $[S]_{0, 1} = 10$ $μ$ $M$ and $[S]_{0, 2} = 1$ $μ$ $M$ for $k_{+ 1, 1} = k_{+ 1, 2} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{- 1, 1} = k_{- 1, 2} = 0.1$ $s^{- 1}$ ; $k_{2, 1} = k_{2, 2} = 0.5$ $s^{- 1}$ , $[E S]_{0, 1} = [E S]_{0, 2} = 0$ $μ$ $M$ , $[P]_{0, 1} = [P]_{0, 2} = 0$ $μ$ $M$ , and $[E]_{0, 1} = [E]_{0, 2} = 1$ $μ$ $M$ . We can see that the initial reaction rate decreases and that $[E S]$ decreases rapidly for $[S]_{0, 2}$ . There is only enough substrate for a single turnover. Where does each curve reach a plateau?
Next, compare $[E]_{0, 1} = 1$ $μ$ $M$ and $[E]_{0, 2} = 0.1$ $μ$ $M$ for $k_{+ 1, 1} = k_{+ 1, 2} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{- 1, 1} = k_{- 1, 2} = 0.1$ $s^{- 1}$ ; $k_{2, 1} = k_{2, 2} = 0.5$ $s^{- 1}$ , $[E S]_{0, 1} = [E S]_{0, 2} = 0$ $μ$ $M$ , $[P]_{0, 1} = [P]_{0, 2} = 0$ $μ$ $M$ , and $[S]_{0, 1} = [S]_{0, 2} = 10$ $μ$ $M$ . We can see that the duration of the steady-state phase increases for $[E]_{0, 2}$ . The enzyme undergoes multiple rounds of catalysis without depleting the substrate. When we study steady-state kinetics, conditions are chosen so that $[S]_{0} ≫ [E]_{0}$ .

Solution to Exercise 85

Compare $[S]_{0, 1} = 10$ $μ$ $M$ and $[S]_{0, 2} = 1$ $μ$ $M$ for $k_{+ 1, 1} = k_{+ 1, 2} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{- 1, 1} = k_{- 1, 2} = 0.1$ $s^{- 1}$ ; $k_{2, 1} = k_{2, 2} = 0.5$ $s^{- 1}$ , $[E S]_{0, 1} = [E S]_{0, 2} = 0$ $μ$ $M$ , $[P]_{0, 1} = [P]_{0, 2} = 0$ $μ$ $M$ , and $[E]_{0, 1} = [E]_{0, 2} = 1$ $μ$ $M$ . We can see that the initial reaction rate decreases and that $[E S]$ decreases rapidly for $[S]_{0, 2}$ . There is only enough substrate for a single turnover. Where does each curve reach a plateau?

CEandStofromEStoP_1 = [1, 10, 0, 0]
CEandStofromEStoP_2 = [1, 1, 0, 0]

yEandStofromEStoP_1 = odeint(funcEandStofromEStoP, CEandStofromEStoP_1, tEandStofromEStoP, ksEandStofromEStoP)
yEandStofromEStoP_2 = odeint(funcEandStofromEStoP, CEandStofromEStoP_2, tEandStofromEStoP, ksEandStofromEStoP)

plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP_1[:,2],'g:') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.plot(tEandStofromEStoP,yEandStofromEStoP_2[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[ES] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,30]) #set X-axis range
plt.ylim([0,1.1]) #set Y-axis range
plt.legend(['[S] = 10 $\mu$M','[S] = 1 $\mu$M'], loc='best') #include a legend
plt.show()

plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP_1[:,3],'g:') #select the fourth column of yAandBtofromAB, i.e. [P], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP_2[:,3],'r--') #select the fourth column of yAandBtofromAB, i.e. [P], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[P] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,30]) #set X-axis range
plt.ylim([0,11]) #set Y-axis range
plt.legend(['[S] = 10 $\mu$M','[S] = 1 $\mu$M'], loc='best') #include a legend
plt.show()

In the first graph ( $[E S]$ versus $t$ ), we see that $[E S]$ decreases rapidly when $[E]_{0} \approx [S]_{0} = 1$ $μ$ $M$ . This is called a single turnover experiment.

In the second graph, ( $[P]$ versus $t$ ), we see that $[P]_{\infty} = [S]_{0}$ .

Next, compare $[E]_{0, 1} = 1$ $μ$ $M$ and $[E]_{0, 2} = 0.1$ $μ$ $M$ for $k_{+ 1, 1} = k_{+ 1, 2} = 1$ $μ$ $M^{- 1} s^{- 1}$ ; $k_{- 1, 1} = k_{- 1, 2} = 0.1$ $s^{- 1}$ ; $k_{2, 1} = k_{2, 2} = 0.5$ $s^{- 1}$ , $[E S]_{0, 1} = [E S]_{0, 2} = 0$ $μ$ $M$ , $[P]_{0, 1} = [P]_{0, 2} = 0$ $μ$ $M$ , and $[S]_{0, 1} = [S]_{0, 2} = 10$ $μ$ $M$ . We can see that the duration of the steady-state phase increases for $[E]_{0, 2}$ . The enzyme undergoes multiple rounds of catalysis without depleting the substrate. When we study steady-state kinetics, conditions are chosen so that $[S]_{0} ≫ [E]_{0}$ .

CEandStofromEStoP_3 = [1, 10, 0, 0]
CEandStofromEStoP_4 = [0.1, 10, 0, 0]

yEandStofromEStoP_3 = odeint(funcEandStofromEStoP, CEandStofromEStoP_3, tEandStofromEStoP, ksEandStofromEStoP)
yEandStofromEStoP_4 = odeint(funcEandStofromEStoP, CEandStofromEStoP_4, tEandStofromEStoP, ksEandStofromEStoP)

plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP_3[:,2],'g:') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.plot(tEandStofromEStoP,yEandStofromEStoP_4[:,2],'r--') #select the third column of yEandStofromEStoP, i.e. [ES], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[ES] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,30]) #set X-axis range
plt.ylim([0,1.1]) #set Y-axis range
plt.legend(['[S] = 1 $\mu$M','[S] = 0.1 $\mu$M'], loc='best') #include a legend
plt.show()

plt.figure(figsize=(7,5))
plt.plot(tEandStofromEStoP,yEandStofromEStoP_3[:,3],'g:') #select the fourth column of yAandBtofromAB, i.e. [P], and plot versus t, use green
plt.plot(tEandStofromEStoP,yEandStofromEStoP_4[:,3],'r--') #select the fourth column of yAandBtofromAB, i.e. [P], and plot versus t, use red
plt.xlabel('Time (s)', fontsize=18) #label the X-axis
plt.ylabel('[P] ($\mu$M)', fontsize=18) #label the Y-axis
plt.xlim([0,30]) #set X-axis range
plt.ylim([0,11]) #set Y-axis range
plt.legend(['[S] = 1 $\mu$M','[S] = 0.1 $\mu$M'], loc='best') #include a legend
plt.show()

From the first ( $[E S]$ versus $t$ ) and second ( $[P]$ versus $t$ ) graphs, we see that when $[E]_{0} = 0.1$ $μ$ $M$ << $[S]_{0} = 10$ $μ$ $M$ , we have a much longer steady-state phase.

Cytochrome c#

Think about a possible mechanism for cytochrome c peroxidase activity. Don’t forget that native cytochrome c needs to be activated to become a more active peroxidase and that cytochrome c peroxidase can become deactivated as well! Which parameters do you need to be able to simulate the reaction?

There is no right answer here! For examples of possible mechanisms for cytochrome c peroxidase activity, see

Parakra RD , Kleffmann T , Jameson GNL , Ledgerwood EC . The proportion of Met80-sulfoxide dictates peroxidase activity of human cytochrome c. Dalton Trans. 2018;47(27):9128-9135. doi:10.1039/c8dt02185f
Yin V, Shaw GS, Konermann L. Cytochrome c as a Peroxidase: Activation of the Precatalytic Native State by H2O2-Induced Covalent Modifications. J Am Chem Soc. 2017;139(44):15701-15709. doi:10.1021/jacs.7b07106

Simulating data

Contents

Simulating data#

Tools#

Why?#

Convention#

Import libraries#

A → B#

Background#

Simulation#

A+B ⇌ AB#

Background#

Simulation#

E+S ⇌ ES → E+P#

Background#

Simulation#

Pre-steady state and steady state phases#

Experimental design: [S]0 relative to [E]0#

Cytochrome c#

$A$ $\to$ $B$ #

$A + B$ $⇌$ $A B$ #

$E + S$ $⇌$ $E S$ $\to$ $E + P$ #

Experimental design: $[S]_{0}$ relative to $[E]_{0}$ #