1 |
"""
|
2 |
Provide the model for curvefitting.
|
3 |
This is a TWO STATE REACTION model with time variable concentration..
|
4 |
The reaction equations are:
|
5 |
L + A ==ka1=> LA
|
6 |
L + A <=kd1== LA
|
7 |
LA ==ka2=> LA*
|
8 |
LA <=kd2== LA*
|
9 |
|
10 |
Rui Hou, Yuhang Wan
|
11 |
Last modified on 100427 (yymmdd) by YW
|
12 |
|
13 |
Typical Pipeline:(work together with "modelclass.py" the father class):
|
14 |
>import modeltwostate_varyC as mtvc
|
15 |
>m = mtvc.twostatemodel_varyC()
|
16 |
>...
|
17 |
|
18 |
"""
|
19 |
__version__ = "100427"
|
20 |
|
21 |
import numpy as np
|
22 |
import pylab as plt
|
23 |
import copy
|
24 |
import os
|
25 |
import pickle
|
26 |
from modelclass import *
|
27 |
|
28 |
|
29 |
class twostatemodel_varyC(modelclass):
|
30 |
"""The data attributes of the class:
|
31 |
parainfo: the parameter information
|
32 |
sim_data: the simulated data
|
33 |
conc: the time variable concentration data
|
34 |
"""
|
35 |
def __init__ (self, parainfo=[], sim_data=[], conc=[] ):
|
36 |
modelclass.__init__(self,parainfo,sim_data)
|
37 |
self.conc = conc
|
38 |
|
39 |
##------------ not finished yet---------------
|
40 |
def load_conc(self, fname):
|
41 |
##fp = file(fname, 'r+')
|
42 |
conc = np.loadtxt(fname)
|
43 |
self.conc = conc
|
44 |
return
|
45 |
|
46 |
def create_conc(self, t):
|
47 |
# create a step-function as concentration curve for debugging now
|
48 |
c = np.zeros(len(t), dtype=float)
|
49 |
conc = np.zeros((2,len(t)), dtype=float)
|
50 |
t = np.array(t)
|
51 |
t1 = input("start time:")
|
52 |
t2 = input("end time:")
|
53 |
ind = plt.find((t>t1-(t[1]-t[0]))&(t<t2-(t[1]-t[0])))
|
54 |
c[ind] = 1
|
55 |
conc[0] = t
|
56 |
conc[1] = c
|
57 |
self.conc = conc
|
58 |
plt.plot(t,c,'.')
|
59 |
plt.show()
|
60 |
return
|
61 |
|
62 |
##----------------for debugging now------------
|
63 |
|
64 |
def wizard(self, ):
|
65 |
'''This function helps you to create a parameter information list.
|
66 |
This is a competing reactions model, where there are 6 parameters:
|
67 |
rmax: Maximum analyte binding capacity(RU),
|
68 |
ka1: Association rate constant for L+A=LA(M-1S-1),
|
69 |
kd1: Dissociation rate constant for LA=L+A(S-1),
|
70 |
ka2: Forward rate constant for LA=LA*(S-1),
|
71 |
kd2: Backward rate constant for LA*=LA(S-1),
|
72 |
ca: Analyte concentration(M).'''
|
73 |
print ('')
|
74 |
pname = ['rmax','ka1','kd1','ka2','kd2','ca']
|
75 |
parainfo = []
|
76 |
for i in range(6):
|
77 |
print i
|
78 |
tmp={}
|
79 |
tmp['name'] = pname[i]
|
80 |
vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
|
81 |
vtmp = vstr.strip().split(',')
|
82 |
tmp['number'] = len(vtmp)
|
83 |
if len(vtmp) == 1:
|
84 |
vtmp = float(vtmp[0])
|
85 |
else:
|
86 |
vtmp = map(float, vtmp)
|
87 |
tmp['value'] = vtmp
|
88 |
tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
|
89 |
|
90 |
parainfo.append(tmp)
|
91 |
##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
|
92 |
self.parainfo = parainfo
|
93 |
|
94 |
|
95 |
def function(self, t, paralist):
|
96 |
'''This function calculates the theoretical curve through the
|
97 |
parameter list you give.
|
98 |
'''
|
99 |
# for two state reaction model
|
100 |
conc = copy.deepcopy(self.conc)
|
101 |
C = np.interp(t,conc[0],conc[1])
|
102 |
## Assign the parameters to calculate the curve
|
103 |
for p in paralist:
|
104 |
if p['name'] == 'rmax': rmax = p['value']
|
105 |
elif p['name'] == 'ka1': ka1 = p['value']
|
106 |
elif p['name'] == 'kd1': kd1 = p['value']
|
107 |
elif p['name'] == 'ka2': ka2 = p['value']
|
108 |
elif p['name'] == 'kd2': kd2 = p['value']
|
109 |
elif p['name'] == 'ca': ca = p['value']
|
110 |
else: print p['name'], p['value']
|
111 |
if type(ca1) == list or type(ca2) == list:
|
112 |
print "Error: This function can only generate data for a single concentration."
|
113 |
return
|
114 |
|
115 |
## Must iterate through data, numerical integration.
|
116 |
g = np.zeros(len(C), dtype=float)
|
117 |
g1 = np.zeros(len(C), dtype=float)
|
118 |
g2 = np.zeros(len(C), dtype=float)
|
119 |
|
120 |
for i in range(1,len(C)):
|
121 |
dG1 = (ka1*ca*C[i-1]*(rmax-g1[i-1]-g2[i-1]) - kd1*g1[i-1])-(ka2*g1[i-1]- kd2*g2[i-1])
|
122 |
dG2 = ka2*g1[i-1]- kd2*g2[i-1]
|
123 |
|
124 |
if (abs(g1[i]) > 999999999): dG1 = 0
|
125 |
if (abs(g2[i]) > 999999999): dG2 = 0
|
126 |
|
127 |
g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
|
128 |
g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
|
129 |
g[i] = g1[i] + g2[i]
|
130 |
|
131 |
return g
|
132 |
|
133 |
##### End of time variable concentrated two state model class definition. |