ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 61
Committed: Thu May 19 21:39:02 2011 UTC (8 years, 6 months ago) by clausted
File size: 17674 byte(s)
Log Message:
Update docstrings so html looks nice.
Line File contents
1 """
2 mdl_module
3 -------------------------------------------------------------------------------
4
5 *Example model functions module for SPRI data.*
6 Just a few simple, common interaction models will go here.
7 In the future, each new model will go in its own module/file.
8 Each model is a function taking three parameters.
9
10 **Parameters:**
11 * *time* (numpy array) -- Time points, usually in seconds.
12 * *data* (numpy array) -- Initial SPR signal values.
13 * *params* (dictionary of dictionaries) -- Model parameter description.
14
15 **Returns:** numpy array of calculated SPR signal values.
16
17 The dictionary keys are the names of model parameters
18 (e.g. 'rmax' for maximal response, or 'kon' for kinetic on-rate).
19 Each model parameter is described by a subdictionary containing four entries.
20
21 * ``'value'`` The current value of the parameter.
22 * ``'min'`` The minimum allowable value.
23 * ``'max'`` The maximum allowable value.
24 * ``'fixed'`` Either 'float', 'fixed', or a reference to another ROI.
25
26 The ``min``, ``max``, and ``fixed`` keys are used during automatic
27 curve-fitting. A fixed parameter is not allowed to change, while a
28 float parameter is adjusted until the least-squares algorithm has
29 minimized the sum-squared error. The *fixed* parameter may also be an
30 integer, in which case it is fixed to the value of a parameter of the
31 same name in another ROI.
32
33 The model function returns an array of values obtained by numerical
34 integration. The model is represented by differential equations and
35 integrated using the rectangle rule or, preferentially, using the
36 trapezoidal rule.
37
38 Functions are also available that provide parameter dictionaries
39 containing all the necessary default values.
40
41 .. moduleauthor:: Christopher Lausted,
42 Institute for Systems Biology,
43 OSPRAI developers.
44
45 Examples::
46
47 >>> import mdl_module as mdl
48 >>> import numpy as np
49 >>> times, data = np.arange(300), np.zeros(300)
50 >>>
51 >>> param1 = {'rate': {'value':1.0, 'min':-10.0, 'max':10.0, 'fixed':True} }
52 >>> data1 = mdl.drift(times, data, param1)
53 >>> print round(data1[299], 1)
54 299.0
55 >>> param1 = dict(rate=dict(value=-1.0, min=-10.0, max=10.0, fixed=True))
56 >>> data1 = mdl.drift(times, data, param1)
57 >>> print round(data1[299], 1)
58 -299.0
59 >>>
60 >>> param2 = mdl.simple1to1_def_params()
61 >>> param2['rmax'] ['value'] = 1000.0
62 >>> param2['conc']['value'] = 1e-6
63 >>> param2['t1']['value'] = 100
64 >>> param2['kon']['value'] = 1e4
65 >>> param2['t2']['value'] = 200.0
66 >>> param2['koff']['value'] = 1e-2
67 >>> param2['t3']['value'] = 300.0
68 >>> data2 = mdl.simple1to1(times, data, param2)
69 >>> print round(data2[299], 1)
70 160.3
71
72 """
73 __version__ = "110517"
74
75
76 ## Import libraries
77 import numpy as np
78 from copy import deepcopy
79
80
81 def drift(time, data, params):
82 """
83 This function simply models a constant signal drift in units/second.
84 It requires numpy arrays of times and starting data values,
85 It only requires one parameter in the params list.
86
87 ``params['rate']['value']``
88 """
89 y = np.zeros(len(time), dtype=float)
90 try:
91 rate = params['rate']['value']
92 except KeyError:
93 print "Error: Rate parameter not supplied to drift() model."
94 return y
95
96 ## Numerical integration.
97 y[0] = data[0]
98 for i in range(1,len(time)):
99 dy = rate * (time[i]-time[i-1])
100 y[i] = y[i-1] + dy
101
102 return y
103 ## End of drift() function.
104
105 def drift_def_params():
106 """
107 Return dictionary of default parameters for drift model.
108 For example::
109
110 >>> from numpy import arange, zeros
111 >>> defaultparam0 = drift_def_params()
112 >>> data0 = drift(arange(100), zeros(100), defaultparam0)
113 >>> print round(data0[99], 1)
114 99.0
115 """
116 dpar = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
117 return dpar
118
119
120 def simple1to1(time, data, params):
121 """
122 This function simply models a 1:1 interaction.
123 The model parameters are described in the table below.
124
125 Model::
126
127 [A] + [L] --kon--> [AL]
128 [A] + [L] <-koff-- [AL]
129
130 Derivation::
131
132 d[AL]/dt = kon*[A]*[L] - koff*[AL]
133 y = [AL]
134 (rmax-y) = [L]
135 conc = [A]
136 rmax = [AL] + [L]
137 dy/dt = conc*kon*(rmax-y) - koff*y
138
139 ======================= ============================================
140 Model Parameter Description
141 ======================= ============================================
142 params['t1']['value'] time of injection for binding, (s)
143 params['rmax']['value'] maximum response, (RIU)
144 params['conc']['value'] concentration of analyte [A], (M)
145 params['kon']['value'] on-rate of analyte, (1/Ms)
146 params['t2']['value'] time of end binding & begin washing, (s)
147 params['koff']['value'] off-rate of analyte, (1/s)
148 params['t3']['value'] time end of washing & data fitting, (s)
149 ----------------------- --------------------------------------------
150 *Optional*
151 params['cofa']['value'] concentration factor, (1/dilution factor)
152 ======================= ============================================
153 """
154
155 ## Skip parameter validation steps for now.
156 t1 = float(params['t1']['value'])
157 rmax = float(params['rmax']['value'])
158 conc = float(params['conc']['value'])
159 kon = float(params['kon']['value'])
160 t2 = float(params['t2']['value'])
161 koff = float(params['koff']['value'])
162 t3 = float(params['t3']['value'])
163 if ('cofa' in params.keys()):
164 conc *= float(params['cofa']['value'])
165
166 ## Initialize variables.
167 errlog = {} ## Error status dictionary.
168 y = np.zeros(len(time), dtype=float)
169 base = y[0] = data[0] ## Baseline SPR signal.
170
171 ## Must iterate through data, numerical integration.
172 for i in range(1,len(time)):
173 dt = time[i] - time[i-1]
174 ## Treat function as having three separate phases.
175 if (time[i] <= t1):
176 ## Pre-binding phase.
177 base = y[i] = data[i]
178 elif (t1 < time[i] <= t2):
179 ## Binding phase
180 yb = y[i-1] - base
181 dy = conc*kon*(rmax-yb) - koff*yb
182 y[i] = y[i-1] + dy*dt
183 ## Check if we overshot the maximum response.
184 if (y[i] > (base+rmax)):
185 y[i] = base+rmax
186 errlog['Response overshot rmax.'] = True
187 elif (t2 < time[i] <= t3):
188 ## Dissociation (conc=0)
189 yb = y[i-1]-base
190 dy = 0 - koff*yb
191 y[i] = y[i-1] + dy*dt
192
193 if any(errlog): print "Errors in simple1to1:", errlog.keys()
194
195 return y
196 """End of simple1to1() function"""
197
198 def simple1to1_def_params():
199 """
200 Return a dictionary of default parameters for simple1to1 model.
201 For example::
202
203 >>> from numpy import arange, zeros
204 >>> defaultparam0 = simple1to1_def_params()
205 >>> data0 = simple1to1(arange(300), zeros(300), defaultparam0)
206 >>> print round(data0[299], 1)
207 160.3
208 """
209 dpar = {}
210 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed=True)
211 dpar['conc'] = dict(value=1e-6, min=1e-12, max=1e-1, fixed=True)
212 dpar['cofa'] = dict(value=1.0, min=1e-6, max=1e6, fixed=True)
213 dpar['t1'] = dict(value=100.0, min=100.0, max=100.0, fixed=True)
214 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed=True)
215 dpar['t2'] = dict(value=200.0, min=200.0, max=200.0, fixed=True)
216 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed=True)
217 dpar['t3'] = dict(value=300.0, min=300.0, max=300.0, fixed=True)
218 return dpar
219
220
221 def simple1to1_mtl(time, data, params):
222 """
223 This function simply models a 1:1 interaction with mass transport limitation.
224 The model parameters are described in the table below.
225
226 Model::
227
228 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
229 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
230
231 Derivation::
232
233 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
234 y = [AL]
235 (rmax-y) = [L]
236 conc = [Abulk]
237 rmax = [AL] + [L]
238 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
239
240 ======================= ============================================
241 Model Parameter Description
242 ======================= ============================================
243 params['t1']['value'] time of injection for binding, (s)
244 params['rmax']['value'] maximum response, (RIU)
245 params['conc']['value'] concentration of analyte [Abulk], (M)
246 params['kon']['value'] on-rate of analyte, (1/Ms)
247 params['t2']['value'] time of end binding & begin washing, (s)
248 params['koff']['value'] off-rate of analyte, (1/s)
249 params['t3']['value'] time end of washing & data fitting, (s)
250 params['kmtl']['value'] rate of diffusion, (RIU/Ms)
251 ----------------------- --------------------------------------------
252 *Optional*
253 params['cofa']['value'] concentration factor, (1/dilution factor)
254 ======================= ============================================
255 """
256
257 ## Skip parameter validation steps for now.
258 t1 = params['t1']['value']
259 rmax = params['rmax']['value']
260 conc = params['conc']['value']
261 kon = params['kon']['value']
262 t2 = params['t2']['value']
263 koff = params['koff']['value']
264 t3 = params['t3']['value']
265 kmtl = params['kmtl']['value']
266 if ('cofa' in params.keys()):
267 conc *= float(params['cofa']['value'])
268
269 ## Initialize variables.
270 errlog = {} ## Error status dictionary.
271 y = np.zeros(len(time), dtype=float)
272 base = y[0] = data[0] ## Baseline SPR signal.
273
274 ## Error checks.
275 if (kmtl<=0):
276 ## Avoid div/0, assume user wants no mass transport limitation.
277 stat['kmtl must be > 0'] = True
278 kmtl = 1e40 ## An arbitrary very large number.
279
280 ## Must iterate through data, numerical integration.
281 for i in range(1,len(time)):
282 dt = time[i] - time[i-1]
283 ## Treat function as having three separate phases.
284 if (time[i] <= t1):
285 ## Pre-binding phase.
286 base = y[i] = data[i]
287 elif (t1 < time[i] <= t2):
288 ## Binding phase
289 yb = y[i-1] - base
290 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
291 y[i] = y[i-1] + dy*dt
292 ## Check if we overshot the maximum response.
293 if (y[i] > (base+rmax)):
294 y[i] = base+rmax
295 errlog['Response overshot rmax.'] = True
296 elif (t2 < time[i] <= t3):
297 ## Dissociation (conc=0)
298 yb = y[i-1] - base
299 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
300 y[i] = y[i-1] + dy*dt
301
302 if any(errlog): print "Errors in simple1to1:", errlog.keys()
303
304 return y
305 """End of simple1to1_mtl() function"""
306
307 def simple1to1_mtl_def_params():
308 """
309 Return dictionary of default parameters for simple1to1_mtl model.
310 For example::
311
312 >>> from numpy import arange, zeros
313 >>> defaultparam0 = simple1to1_mtl_def_params()
314 >>> data0 = simple1to1_mtl(arange(300), zeros(300), defaultparam0)
315 >>> print round(data0[299], 1)
316 161.1
317 """
318 dpar = simple1to1_def_params() ## Based on simple1to1.
319 dpar['kmtl'] = dict(value=1e9, min=1e2, max=1e10, fixed=True)
320 return dpar
321
322
323 def simple1to1_series(time, data, params):
324 """
325 Model up to 8 simple 1:1 SPR injections of varying concentration.
326
327 ======================= ==========================================
328 Model Parameter Description
329 ======================= ==========================================
330 params['kon']['value'] on-rate of analyte, (1/Ms)
331 params['koff']['value'] off-rate of analyte, (1/s)
332 params['rmax']['value'] maximum response, (RIU)
333 params['tp']['value'] time prior to binding, (s)
334 params['ta']['value'] time/duration of association, (s)
335 params['td']['value'] time/duration of dissociation, (s)
336 ----------------------- ------------------------------------------
337 params['t0']['value'] timepoint of first (no.0) injection, (s)
338 ... ...
339 params['t7']['value'] timepoint of eighth (no.7) injection, (s)
340 ----------------------- ------------------------------------------
341 params['c0']['value'] concentration of first injection, (s)
342 ... ...
343 params['c7']['value'] concentration of eighth injection, (s)
344 ======================= ==========================================
345
346 Injections where tN=0 are ignored. For example, if you have only
347 six injections, set t6=0 and t7=0.
348 """
349
350 ## Parameters from params dictionary.
351 ## Skip parameter validation steps for now.
352 kon = float(params['kon']['value'])
353 koff = float(params['koff']['value'])
354 rmax = float(params['rmax']['value'])
355 ## Pre-association, Association duration, Dissociation duration.
356 tp = float(params['tp']['value']) ## Typical -100 s
357 ta = float(params['ta']['value']) ## Typical 300 s
358 td = float(params['td']['value']) ## Typical 300 s
359 td = td + ta ## After 300s binding and 300s washing, 600s elapse.
360 ## How many injections are there? Check valid keys and inj times.
361 tkey = ['t0','t1','t2','t3','t4','t5','t6','t7']
362 ckey = ['c0','c1','c2','c3','c4','c5','c6','c7']
363 for i in range(8):
364 if (tkey[i] in params) and (ckey[i] in params):
365 ## Valid inj time is after sufficient pre-association time.
366 if ((float(params[tkey[i]]['value']) + tp) > 0):
367 inj = i+1 ## Increase inj count.
368 ## Up to eight injection start times, concentrations.
369 ti = [float(params[k]['value']) for k in tkey[0:inj]]
370 c = [float(params[k]['value']) for k in ckey[0:inj]]
371
372 ## Begin.
373 y = np.zeros(len(time), dtype=float)
374 base = y[0] = data[0] ## Baseline SPR signal.
375
376 ## Must iterate through data, numerical integration.
377 for i in range(1,len(time)):
378 dt = time[i] - time[i-1]
379 #if (time[i] > (max(ti)+td)): break ## Speed things up.
380 inmodel = False
381 ## For-loop where j is each injection number.
382 ## Treat function as having three separate phases.
383 for j in range(inj):
384 if ((ti[j]+tp) <= time[i] <= (ti[j])):
385 ## Pre-binding phase.
386 if np.isnan(y[i-1]): y[i-1]=0
387 y[i] = (0.6 * y[i-1]) + (0.4 * data[i])
388 base = y[i]
389 inmodel = True
390 if ((ti[j]) <= time[i] <= (ti[j]+ta)):
391 ## Binding phase
392 yb = y[i-1] - base
393 dy = c[j]*kon*(rmax-yb) - koff*yb
394 y[i] = y[i-1] + dy*dt
395 inmodel = True
396 if (y[i] > (base+rmax)): y[i] = (base+rmax) ## Avoid overshoot.
397 if ((ti[j]+ta) <= time[i] <= (ti[j]+td)):
398 ## Dissociation (conc=0)
399 yb = y[i-1]-base
400 dy = 0 - koff*yb
401 y[i] = y[i-1] + dy*dt
402 inmodel = True
403 ## End "for j" loop.
404 if (inmodel==False):
405 ## Not modeled, e.g. regeneration.
406 y[i] = np.NaN
407 y[i] = 0 ## NaN might break mclma()?
408 ## End "for i" loop.
409 return y
410
411 def simple1to1_series_def_params():
412 """
413 Return a dictionary of default parameters for simple1to1_series model.
414 While parameters for eight injection times are provided, only the
415 first injection time is nonzero. For example::
416
417 >>> from numpy import arange, zeros
418 >>> defaultparam0 = simple1to1_series_def_params()
419 >>> data0 = simple1to1_series(arange(300), zeros(300), defaultparam0)
420 >>> print round(data0[298], 1)
421 160.3
422 """
423 dpar = {}
424 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed=True)
425 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed=True)
426 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed=True)
427 dpar['tp'] = dict(value=-50.0, min=-50.0, max=-50.0, fixed=True)
428 dpar['ta'] = dict(value=100.0, min=100.0, max=100.0, fixed=True)
429 dpar['td'] = dict(value=100.0, min=100.0, max=100.0, fixed=True)
430 ## Set all eight concentrations to 1uM.
431 for key in ['c0','c1','c2','c3','c4','c5','c6','c7']:
432 dpar[key] = dict(value=1e-6, min=1e-12, max=1e-1, fixed=True)
433 ## Just describe the first injection at t0=100s.
434 for key in ['t0','t1','t2','t3','t4','t5','t6','t7']:
435 dpar[key] = dict(value=0.0, min=0.0, max=0.0, fixed=True)
436 dpar['t0'] = dict(value=100.0, min=100.0, max=100.0, fixed=True)
437 return dpar
438
439
440 def _test():
441 """
442 Automatic Code testing with doctest.
443 Doctest runs the example code in the docstrings.
444 Enable use of ellipses as wildcard for returned text.
445 """
446 import doctest
447 doctest.testmod(optionflags=doctest.ELLIPSIS)
448
449
450 if __name__ == "__main__":
451 """
452 Code testing.
453 Simply execute 'python mdl_module.py'.
454 """
455 _test()
456
457 ################################# End of module ###############################