Before we start...¶
We need to install packages required to calculate. To do so, we need to use Terminal as following instruction.
It should be installed but just in case, if you attempt to execute this script at somewhere else
- (assuming you are in the same env as "open lab") Open Terminal from Launcher
- Click Terminal then in the prompt, type the followings
pip install scipy matplotlib ipywidgets ipympl
Or execute the following line:
In [1]:
!pip install scipy matplotlib
In [2]:
## These are the python libraries you will utilize for the calculation
import matplotlib.pylab as plt
import matplotlib.pyplot as plt; plt.rcdefaults()
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
%matplotlib inline
import numpy as np
import scipy
from scipy.optimize import curve_fit
from scipy.integrate import odeint
Ordinary Differential Equation¶
$\dfrac{d[A]}{dt}=-k_1[A]+k_2[B]$
$\dfrac{d[B]}{dt}=-(k_2+k_3)[B]+k_1[A]+k_4[C]$
$\dfrac{d[C]}{dt}=-(k_4+k_5)[C]+k_3[B]+k_6[D]$
$\dfrac{d[D]}{dt}=-k_6[D]+k_5[C]$
- create a function to return $\dfrac{dy}{dt}$
- set parameters (ex. $k_1 = 1$)
- set length of calculation (this is a time-dependent)
- execute ode solver
- plot them
In [3]:
def func1(y,t):
k1 = 1
k2 = 0.5
k3 = 1.2
k4 = 0.7
k5 = 0.8
k6 = 0.7
dAdt = -k1*y[0] + k2*y[1]
dBdt = -(k2+k3)*y[1] + k1*y[0] + k4*y[2]
dCdt = -(k4+k5)*y[2] + k3*y[1] + k6*y[3]
dDdt = -k6*y[3] + k5*y[2]
dydt = [dAdt, dBdt, dCdt, dDdt]
return dydt
In [4]:
def func2(y,t,k):
k1 = k[0]
k2 = k[1]
k3 = k[2]
k4 = k[3]
k5 = k[4]
k6 = k[5]
dAdt = -k1*y[0] + k2*y[1]
dBdt = -(k2+k3)*y[1] + k1*y[0] + k4*y[2]
dCdt = -(k4+k5)*y[2] + k3*y[1] + k6*y[3]
dDdt = -k6*y[3] + k5*y[2]
dydt = [dAdt, dBdt, dCdt, dDdt]
return dydt
In [5]:
def func3(y,t,k):
A,B,C,D = y
k1,k2,k3,k4,k5,k6 = k
dAdt = -k1*A + k2*B
dBdt = -(k2+k3)*B + k1*A + k4*C
dCdt = -(k4+k5)*C + k3*B + k6*D
dDdt = -k6*D + k5*C
dydt = [dAdt, dBdt, dCdt, dDdt]
return dydt
In [6]:
def func4(y,t,J):
return np.dot(J,y)
In [7]:
# Initial condition of calculation
y0 = [10,0,0,0]
# for function 1
# the set of parameter is predetermined
# for function 2,3
# k (reaction rate or transition rate) is one of input variables you need to specify before calculation
# function 2 and 3 handle the list of reaction rate
k = [1, 0.5, 1.2, 0.7, 0.8, 0.7]
# for function 4
# we are using Jacobian method so we will introduce Jacobian matrix using np.array and the function will perform dot operation
k1 = k[0]
k2 = k[1]
k3 = k[2]
k4 = k[3]
k5 = k[4]
k6 = k[5]
Jacobian = np.array([[-k1, k2, 0, 0],
[k1, -(k2+k3), k4, 0],
[0, k3, -(k4+k5), k6],
[0, 0, k5, -k6]])
print(Jacobian)
In [8]:
# calculation time
time_start = 0
time_end = 10
NoOfSteps = 100
time = np.linspace(time_start,time_end,NoOfSteps)
In [9]:
y1 = odeint(func1,y0,time)
y2 = odeint(func2,y0,time,args=(k,))
y3 = odeint(func3,y0,time,args=(k,))
y4 = odeint(func4,y0,time,args=(Jacobian,))
In [10]:
print('shape of the outcome y1 =', np.shape(y1))
print('shape of the outcome y2 =', np.shape(y2))
print('shape of the outcome y3 =', np.shape(y3))
print('shape of the outcome y4 =', np.shape(y4))
In [11]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='func1')
plt.plot(time,y2[:,0],'b-',label='func2')
plt.plot(time,y3[:,0],'g-',label='func3')
plt.plot(time,y4[:,0],'m-',label='func4')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[11]:
In [12]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='A')
plt.plot(time,y1[:,1],'b-',label='B')
plt.plot(time,y1[:,2],'g-',label='C')
plt.plot(time,y1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[12]:
Use of Dictionary¶
In [13]:
def genJacobian(k):
k1 = k['k1']
k2 = k['k2']
k3 = k['k3']
k4 = k['k4']
k5 = k['k5']
k6 = k['k6']
Jacobian = np.array([[-k1, k2, 0, 0],
[k1, -(k2+k3), k4, 0],
[0, k3, -(k4+k5), k6],
[0, 0, k5, -k6]])
return Jacobian
In [14]:
kdict = {'k1':1,'k2':0.5,'k3':1.2,'k4':0.7,'k5':0.8,'k6':0.7}
J = genJacobian(kdict)
y4_1 = odeint(func4,y0,time,args=(J,))
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y4_1[:,0],'r-',label='A')
plt.plot(time,y4_1[:,1],'b-',label='B')
plt.plot(time,y4_1[:,2],'g-',label='C')
plt.plot(time,y4_1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[14]:
Let's try something simpler¶
Creating your own library/package using .py file
In [15]:
import functions as fn
y1 = odeint(fn.func1,y0,time)
y2 = odeint(fn.func2,y0,time,args=(k,))
y3 = odeint(fn.func3,y0,time,args=(k,))
y4 = odeint(fn.func4,y0,time,args=(Jacobian,))
In [16]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='A')
plt.plot(time,y1[:,1],'b-',label='B')
plt.plot(time,y1[:,2],'g-',label='C')
plt.plot(time,y1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[16]:
In [ ]: