ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 29
Committed: Thu May 20 01:14:55 2010 UTC (9 years ago) by clausted
File size: 6940 byte(s)
Log Message:
Add a parameter ('cofa') to the 1to1 models that represents a concentration factor.  This allows for easier analysis of a dilution series.  For example, if two ROIs contain data for the same sample but at a dilutions of 1:2 and 1:4, the set 'conc' value to the molar stock concentration and set the two 'cofa' values to 0.5 and 0.25.
Line File contents
1 """
2 mdl: Example model functions module for SPRI data.
3 Christopher Lausted, Institute for Systems Biology,
4 OSPRAI developers
5 Last modified on 100518 (yymmdd)
6
7 Examples:
8 #import mdl_module as mdl
9 #import numpy as np
10 #times = np.arange(100)
11 #data = np.zeros(100)
12 #
13 #params1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
14 #data1 = drift(time, data, param1)
15 #
16 #param2 = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
17 #param2['rmax'] = {'value': 100.0}
18 #param2['conc'] = {'value': 1e-6}
19 #param2['kon'] = {'value': 2e4}
20 #param2['t2'] = {'value': 150.0}
21 #param2['koff'] = {'value': 1e-3}
22 #param2['t3'] = {'value': 270.0}
23 #data2 = simple1to1(time, data, param2)
24 """
25 __version__ = "100518"
26
27
28 ## Import libraries
29 import numpy as np
30 from copy import deepcopy
31
32
33 def drift(time, data, params):
34 """
35 This function simply models a constant signal drift in units/second.
36 It requires numpy arrays of times and starting data values,
37 It only requires one parameter in the params list.
38 params['rate']['value']
39 """
40 y = np.zeros(len(time), dtype=float)
41 try:
42 rate = params['rate']['value']
43 except KeyError:
44 print "Error: Rate parameter not supplied to drift() model."
45 return y
46
47 ## Numerical integration.
48 y[0] = data[0]
49 for i in range(1,len(time)):
50 dy = rate * (time[i]-time[i-1])
51 y[i] = y[i-1] + dy
52
53 return y
54 ## End of drift() function.
55
56
57 def simple1to1(time, data, params):
58 """
59 This function simply models a 1:1 interaction.
60 [A] + [L] --kon--> [AL]
61 [A] + [L] <-koff-- [AL]
62
63 It requires numpy arrays of times and starting data values.
64 params['t1']['value'] is time of injection for binding, (s)
65 params['rmax']['value'] is maximum response, (RIU)
66 params['conc']['value'] is stock concentration of analyte [A], (M)
67 params['kon']['value'] is on-rate of analyte, (1/Ms)
68 params['t2']['value'] is time of end binding & begin washing, (s)
69 params['koff']['value'] is off-rate of analyte, (1/s)
70 params['t3']['value'] is time end of washing & data fitting, (s)
71
72 It also can take an optional parameter useful for dilution series.
73 params['cofa']['value'] is concentration factor, (1/dilution factor)
74 """
75 ## Skip parameter validation steps for now.
76 t1 = float(params['t1']['value'])
77 rmax = float(params['rmax']['value'])
78 conc = float(params['conc']['value'])
79 kon = float(params['kon']['value'])
80 t2 = float(params['t2']['value'])
81 koff = float(params['koff']['value'])
82 t3 = float(params['t3']['value'])
83 if ('cofa' in params.keys()):
84 conc *= float(params['cofa']['value'])
85
86 """
87 Derivation:
88 d[AL]/dt = kon*[A]*[L] - koff*[AL]
89 y = [AL]
90 (rmax-y) = [L]
91 conc = [A]
92 rmax = [AL] + [L]
93 dy/dt = conc*kon*(rmax-y) - koff*y
94 """
95
96 y = np.zeros(len(time), dtype=float)
97 base = y[0] = data[0] ## Baseline SPR signal.
98
99 ## Must iterate through data, numerical integration.
100 for i in range(1,len(time)):
101 dt = time[i] - time[i-1]
102 ## Treat function as having three separate phases.
103 if (time[i] <= t1):
104 ## Pre-binding phase.
105 base = y[i] = data[i]
106 elif (t1 < time[i] <= t2):
107 ## Binding phase
108 yb = y[i-1] - base
109 dy = conc*kon*(rmax-yb) - koff*yb
110 if (abs(y[i-1]) > 999999999): dy = 0 ## Is this useful?
111 y[i] = y[i-1] + dy*dt
112 elif (t2 < time[i] <= t3):
113 ## Dissociation (conc=0)
114 yb = y[i-1]-base
115 dy = 0 - koff*yb
116 y[i] = y[i-1] + dy*dt
117
118 return y
119 ## End of simple1to1() function
120
121
122 def simple1to1_mtl(time, data, params):
123 """
124 This function simply models a 1:1 interaction with mass transport limitation.
125 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
126 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
127
128 It requires numpy arrays of times and starting data values,
129 params['t1']['value'] is time of injection for binding, (s)
130 params['rmax']['value'] is maximum response, (RIU)
131 params['conc']['value'] is concentration of analyte [Abulk], (M)
132 params['kon']['value'] is on-rate of analyte, (1/Ms)
133 params['t2']['value'] is time of end binding & begin washing, (s)
134 params['koff']['value'] is off-rate of analyte, (1/s)
135 params['t3']['value'] is time end of washing & data fitting, (s)
136 params['kmtl']['value'] is rate of diffusion, (RIU/Ms)
137
138 It also can take an optional parameter useful for dilution series.
139 params['cofa']['value'] is concentration factor, (1/dilution factor)
140 """
141 ## Skip parameter validation steps for now.
142 t1 = params['t1']['value']
143 rmax = params['rmax']['value']
144 conc = params['conc']['value']
145 kon = params['kon']['value']
146 t2 = params['t2']['value']
147 koff = params['koff']['value']
148 t3 = params['t3']['value']
149 kmtl = params['kmtl']['value']
150 if ('cofa' in params.keys()):
151 conc *= float(params['cofa']['value'])
152
153 """
154 Derivation:
155 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
156 y = [AL]
157 (rmax-y) = [L]
158 conc = [Abulk]
159 rmax = [AL] + [L]
160 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
161 """
162
163 stat = {} ## Error status dictionary.
164 y = np.zeros(len(time), dtype=float)
165 base = y[0] = data[0] ## Baseline SPR signal.
166
167 ## Error checks.
168 if (kmtl<=0):
169 ## Avoid div/0, assume user wants no mass transport limitation.
170 stat['kmtl must be > 0'] = True
171 kmtl = 1e40 ## An arbitrary very large number.
172
173 ## Must iterate through data, numerical integration.
174 for i in range(1,len(time)):
175 dt = time[i] - time[i-1]
176 ## Treat function as having three separate phases.
177 if (time[i] <= t1):
178 ## Pre-binding phase.
179 base = y[i] = data[i]
180 elif (t1 < time[i] <= t2):
181 ## Binding phase
182 yb = y[i-1] - base
183 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
184 ## Check for integration errors producing wild results.
185 if (abs(y[i-1]) > 9e9): dy = 0; stat['y>9e9'] = True
186 y[i] = y[i-1] + dy*dt
187 elif (t2 < time[i] <= t3):
188 ## Dissociation (conc=0)
189 yb = y[i-1] - base
190 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
191 y[i] = y[i-1] + dy*dt
192
193 if (len(stat.keys()) > 0): print "Errors in simple1to1:", stat.keys()
194
195 return y
196 ## End of simple1to1_mtl() function
197
198
199
200 ################################# End of module #################################