Solution for different modelling mechanisms

Reaction: A + B = C

Jens Hahn - 21/05/2015

1. Deterministic Simulation (ODE)

Ordinary differential equations

Kinetic Rate Law: Mass Action

Initial Values $$[A](0) = 3.0 \textrm{ mM}$$ $$[B](0) = 2.0 \textrm{ mM}$$ $$[C](0) = 2.5 \textrm{ mM}$$

$$\textrm{k}_1 = 0.2 \ \frac{1}{\textrm{mM}\times\textrm{s}}$$$$\textrm{k}_2 = 0.1 \ \frac{1}{\textrm{mM}}$$

Reactions $$\textrm{v}_1 = \textrm{k}_1 \times [A](\textrm{t}) \times [B](\textrm{t})$$

$$\textrm{v}_2 = \textrm{k}_2 \times [C](\textrm{t})$$

Differential Equations $$\frac{\textrm{d}[A]}{\textrm{dt}} = - \textrm{v}_1 + \textrm{v}_2$$

$$\frac{\textrm{d}[B]}{\textrm{dt}} = - \textrm{v}_1 + \textrm{v}_2$$$$\frac{\textrm{d}[C]}{\textrm{dt}} = + \textrm{v}_1 - \textrm{v}_2$$



In [1]:

    
# initial values
time_range = range(0,100)
k1 = 0.2
k2 = 0.1
conc_A = 3.0
conc_B = 2.0
conc_C = 2.5

# function to return changes
def eval_change():
    # reactions
    v1 = k1 * conc_A * conc_B
    v2 = k2 * conc_C
    dadt = -v1 + v2
    dbdt = -v1 + v2
    dcdt = v1 - v2
    return dadt, dbdt, dcdt

# solutions for plotting
solution_ODE = {}
solution_ODE['time'] = time_range
solution_ODE['A'] = []
solution_ODE['B'] = []
solution_ODE['C'] = []

# simulation
for time in time_range:
    solution_ODE['A'].append(conc_A)
    solution_ODE['B'].append(conc_B)
    solution_ODE['C'].append(conc_C)
    
    dadt, dbdt, dcdt = eval_change()
    conc_A += dadt
    conc_B += dbdt
    conc_C += dcdt

Create Plot



In [2]:

    
import matplotlib.pyplot as plt
plt.title('My ODE solution')
plt.plot(solution_ODE['time'], solution_ODE['A'], label='[A]')
plt.plot(solution_ODE['time'], solution_ODE['B'], label='[B]')
plt.plot(solution_ODE['time'], solution_ODE['C'], label='[C]')
plt.xlabel('time [sec]')
plt.ylabel('concentration [mM]')
plt.legend(loc='best', prop={'size': 10})
plt.show()

2. Stochastic Simulation (Gillespie SSA)

Gillespie Stochastic Simulation Algorithm

Inital values

$$A = 30$$$$B = 20$$$$C = 25$$

Combinations h (dependent on molecules involved):

Reactions
R1: second order
$\textrm{A} + \textrm{B} \Rightarrow \textrm{h}_1 = \textrm{N}_A \times \textrm{N}_B$

R2: first order
$\textrm{C} \phantom{ + \textrm{B}} \Rightarrow \textrm{h}_2 = \textrm{N}_C$

Propensity a = h k

$$\textrm{k}_1 = 0.2$$$$\textrm{k}_2 = 0.1$$$$\textrm{a}_1 = \textrm{k}_1 * \textrm{h}_1$$$$\textrm{a}_2 = \textrm{k}_2 * \textrm{h}_2$$



In [4]:

    
import random
import math

# initial values
t_start = 0.
t_end = 100.
num_A = 30
num_B = 20
num_C = 25
k1 = 0.2
k2 = 0.1

solution_SSA = {}
solution_SSA['time'] = []
solution_SSA['A'] = []
solution_SSA['B'] = []
solution_SSA['C'] = []

# simulation
time = t_start
while time <= t_end:
    a1 = num_A * num_B * k1
    a2 = num_C * k2
    a = a1 + a2
    norm_a1 = a1/a
    norm_a2 = a2/a
    
    # plotting
    solution_SSA['time'].append(time)
    solution_SSA['A'].append(num_A)
    solution_SSA['B'].append(num_B)
    solution_SSA['C'].append(num_C)
    
    n1 = random.random()
    
    if n1 <= norm_a1:
        num_A -= 1
        num_B -= 1
        num_C += 1
    else:
        num_A += 1
        num_B += 1
        num_C -= 1
        
    n2 = random.random()
    time += -(1/a) * math.log(n2)

Create Plot



In [5]:

    
import matplotlib.pyplot as plt
plt.title('My Gillespie solution')
plt.plot(solution_SSA['time'], solution_SSA['A'], label='[A]')
plt.plot(solution_SSA['time'], solution_SSA['B'], label='[B]')
plt.plot(solution_SSA['time'], solution_SSA['C'], label='[C]')
plt.xlabel('time [sec]')
plt.ylabel('number of molecules')
plt.legend(loc='best', prop={'size': 10})
plt.show()



In [ ]: