ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 65
Committed: Fri May 20 16:33:02 2011 UTC (8 years ago) by clausted
File size: 18709 byte(s)
Log Message:
Added much testing to all modules using docstrings and doctest.  To help this, I also added two optional parameters to mclma(), verbose and maxiter.  For example, use "mclma(rois, verbose=False, maxiter=99)" to run the regression up to 99 iterations without echoing progress to the standard output.

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, numpy as np
48 >>> times, data = np.arange(300), np.zeros(300)
49 >>>
50 >>> param1 = {'rate': {'value':1.0, 'min':-10.0, 'max':10.0, 'fixed':True} }
51 >>> data1 = mdl.drift(times, data, param1)
52 >>> print round(data1[299], 1)
53 299.0
54 >>> param1 = dict(rate=dict(value=-1.0, min=-10.0, max=10.0, fixed=True))
55 >>> data1 = mdl.drift(times, data, param1)
56 >>> print round(data1[299], 1)
57 -299.0
58 >>>
59 >>> param2 = mdl.simple1to1_def_params()
60 >>> param2['rmax'] ['value'] = 1000.0
61 >>> param2['conc']['value'] = 1e-6
62 >>> param2['t1']['value'] = 100
63 >>> param2['kon']['value'] = 1e4
64 >>> param2['t2']['value'] = 200.0
65 >>> param2['koff']['value'] = 1e-2
66 >>> param2['t3']['value'] = 300.0
67 >>> data2 = mdl.simple1to1(times, data, param2)
68 >>> print round(data2[299], 1)
69 160.3
70
71 """
72 __version__ = "110519"
73
74
75 ## Import libraries
76 import numpy as np
77 from copy import deepcopy
78
79
80 def noise(time, data, params):
81 """
82 This function simply adds gaussian noise.
83 It requires numpy arrays of times and starting data values,
84 It only requires one parameter in the params list,
85 ``params['rms']['value']``.
86 """
87 noise = float(params['rms']['value'])
88 return (noise * np.random.randn(len(data)) + data)
89
90
91 def drift(time, data, params):
92 """
93 This function simply models a constant signal drift in units/second.
94 It requires numpy arrays of times and starting data values.
95 It only requires one parameter in the params list,
96 ``params['rate']['value']``.
97 """
98 y = np.zeros(len(time), dtype=float)
99 try:
100 rate = params['rate']['value']
101 except KeyError:
102 print "Error: Rate parameter not supplied to drift() model."
103 return y
104
105 ## Numerical integration.
106 y[0] = data[0]
107 for i in range(1,len(time)):
108 dy = rate * (time[i]-time[i-1])
109 y[i] = y[i-1] + dy
110
111 return y
112 ## End of drift() function.
113
114 def drift_def_params():
115 """
116 Return dictionary of default parameters for drift model.
117 For example::
118
119 >>> from numpy import arange, zeros
120 >>> defaultparam0 = drift_def_params()
121 >>> data0 = drift(arange(100), zeros(100), defaultparam0)
122 >>> print round(data0[99], 1)
123 99.0
124 """
125 dpar = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
126 return dpar
127
128
129 def simple1to1(time, data, params):
130 """
131 This function simply models a 1:1 interaction.
132 The model parameters are described in the table below.
133
134 Model::
135
136 [A] + [L] --kon--> [AL]
137 [A] + [L] <-koff-- [AL]
138
139 Derivation::
140
141 d[AL]/dt = kon*[A]*[L] - koff*[AL]
142 y = [AL]
143 (rmax-y) = [L]
144 conc = [A]
145 rmax = [AL] + [L]
146 dy/dt = conc*kon*(rmax-y) - koff*y
147
148 ======================= ============================================
149 Model Parameter Description
150 ======================= ============================================
151 params['t1']['value'] time of injection for binding, (s)
152 params['rmax']['value'] maximum response, (RIU)
153 params['conc']['value'] concentration of analyte [A], (M)
154 params['kon']['value'] on-rate of analyte, (1/Ms)
155 params['t2']['value'] time of end binding & begin washing, (s)
156 params['koff']['value'] off-rate of analyte, (1/s)
157 params['t3']['value'] time end of washing & data fitting, (s)
158 ----------------------- --------------------------------------------
159 *Optional*
160 params['cofa']['value'] concentration factor, (1/dilution factor)
161 ======================= ============================================
162 """
163
164 ## Skip parameter validation steps for now.
165 t1 = float(params['t1']['value'])
166 rmax = float(params['rmax']['value'])
167 conc = float(params['conc']['value'])
168 kon = float(params['kon']['value'])
169 t2 = float(params['t2']['value'])
170 koff = float(params['koff']['value'])
171 t3 = float(params['t3']['value'])
172 if ('cofa' in params.keys()):
173 conc *= float(params['cofa']['value'])
174
175 ## Initialize variables.
176 errlog = {} ## Error status dictionary.
177 y = np.zeros(len(time), dtype=float)
178 base = y[0] = data[0] ## Baseline SPR signal.
179
180 ## Must iterate through data, numerical integration.
181 for i in range(1,len(time)):
182 dt = time[i] - time[i-1]
183 ## Treat function as having three separate phases.
184 if (time[i] <= t1):
185 ## Pre-binding phase.
186 base = y[i] = data[i]
187 elif (t1 < time[i] <= t2):
188 ## Binding phase
189 yb = y[i-1] - base
190 dy = conc*kon*(rmax-yb) - koff*yb
191 y[i] = y[i-1] + dy*dt
192 ## Check if we overshot the maximum response.
193 if (y[i] > (base+rmax)):
194 y[i] = base+rmax
195 errlog['Response overshot rmax.'] = True
196 elif (t2 < time[i] <= t3):
197 ## Dissociation (conc=0)
198 yb = y[i-1]-base
199 dy = 0 - koff*yb
200 y[i] = y[i-1] + dy*dt
201
202 if any(errlog): print "Errors in simple1to1:", errlog.keys()
203
204 return y
205 """End of simple1to1() function"""
206
207 def simple1to1_def_params():
208 """
209 Return a dictionary of default parameters for simple1to1 model.
210 For example::
211
212 >>> from numpy import arange, zeros
213 >>> defaultparam0 = simple1to1_def_params()
214 >>> data0 = simple1to1(arange(300), zeros(300), defaultparam0)
215 >>> print round(data0[299], 1)
216 160.3
217 """
218 dpar = {}
219 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed='fixed')
220 dpar['conc'] = dict(value=1e-6, min=1e-12, max=1e-1, fixed='fixed')
221 dpar['cofa'] = dict(value=1.0, min=1e-6, max=1e6, fixed='fixed')
222 dpar['t1'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
223 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed='fixed')
224 dpar['t2'] = dict(value=200.0, min=200.0, max=200.0, fixed='fixed')
225 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed='fixed')
226 dpar['t3'] = dict(value=300.0, min=300.0, max=300.0, fixed='fixed')
227 return dpar
228
229
230 def simple1to1_mtl(time, data, params):
231 """
232 This function simply models a 1:1 interaction with mass transport limitation.
233 The model parameters are described in the table below.
234
235 Model::
236
237 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
238 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
239
240 Derivation::
241
242 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
243 y = [AL]
244 (rmax-y) = [L]
245 conc = [Abulk]
246 rmax = [AL] + [L]
247 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
248
249 ======================= ============================================
250 Model Parameter Description
251 ======================= ============================================
252 params['t1']['value'] time of injection for binding, (s)
253 params['rmax']['value'] maximum response, (RIU)
254 params['conc']['value'] concentration of analyte [Abulk], (M)
255 params['kon']['value'] on-rate of analyte, (1/Ms)
256 params['t2']['value'] time of end binding & begin washing, (s)
257 params['koff']['value'] off-rate of analyte, (1/s)
258 params['t3']['value'] time end of washing & data fitting, (s)
259 params['kmtl']['value'] rate of diffusion, (RIU/Ms)
260 ----------------------- --------------------------------------------
261 *Optional*
262 params['cofa']['value'] concentration factor, (1/dilution factor)
263 ======================= ============================================
264 """
265
266 ## Skip parameter validation steps for now.
267 t1 = params['t1']['value']
268 rmax = params['rmax']['value']
269 conc = params['conc']['value']
270 kon = params['kon']['value']
271 t2 = params['t2']['value']
272 koff = params['koff']['value']
273 t3 = params['t3']['value']
274 kmtl = params['kmtl']['value']
275 if ('cofa' in params.keys()):
276 conc *= float(params['cofa']['value'])
277
278 ## Initialize variables.
279 errlog = {} ## Error status dictionary.
280 y = np.zeros(len(time), dtype=float)
281 base = y[0] = data[0] ## Baseline SPR signal.
282
283 ## Error checks.
284 if (kmtl<=0):
285 ## Avoid div/0, assume user wants no mass transport limitation.
286 stat['kmtl must be > 0'] = True
287 kmtl = 1e40 ## An arbitrary very large number.
288
289 ## Must iterate through data, numerical integration.
290 for i in range(1,len(time)):
291 dt = time[i] - time[i-1]
292 ## Treat function as having three separate phases.
293 if (time[i] <= t1):
294 ## Pre-binding phase.
295 base = y[i] = data[i]
296 elif (t1 < time[i] <= t2):
297 ## Binding phase
298 yb = y[i-1] - base
299 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
300 y[i] = y[i-1] + dy*dt
301 ## Check if we overshot the maximum response.
302 if (y[i] > (base+rmax)):
303 y[i] = base+rmax
304 errlog['Response overshot rmax.'] = True
305 elif (t2 < time[i] <= t3):
306 ## Dissociation (conc=0)
307 yb = y[i-1] - base
308 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
309 y[i] = y[i-1] + dy*dt
310
311 if any(errlog): print "Errors in simple1to1:", errlog.keys()
312
313 return y
314 """End of simple1to1_mtl() function"""
315
316 def simple1to1_mtl_def_params():
317 """
318 Return dictionary of default parameters for simple1to1_mtl model.
319 For example::
320
321 >>> from numpy import arange, zeros
322 >>> defaultparam0 = simple1to1_mtl_def_params()
323 >>> data0 = simple1to1_mtl(arange(300), zeros(300), defaultparam0)
324 >>> print round(data0[299], 1)
325 161.1
326 """
327 dpar = simple1to1_def_params() ## Based on simple1to1.
328 dpar['kmtl'] = dict(value=1e9, min=1e2, max=1e10, fixed='fixed')
329 return dpar
330
331
332 def simple1to1_series(time, data, params):
333 """
334 Model up to 8 simple 1:1 SPR injections of varying concentration.
335
336 ======================= ==========================================
337 Model Parameter Description
338 ======================= ==========================================
339 params['kon']['value'] on-rate of analyte, (1/Ms)
340 params['koff']['value'] off-rate of analyte, (1/s)
341 params['rmax']['value'] maximum response, (RIU)
342 params['tp']['value'] time prior to binding, (s)
343 params['ta']['value'] time/duration of association, (s)
344 params['td']['value'] time/duration of dissociation, (s)
345 ----------------------- ------------------------------------------
346 params['t0']['value'] timepoint of first (no.0) injection, (s)
347 ... ...
348 params['t7']['value'] timepoint of eighth (no.7) injection, (s)
349 ----------------------- ------------------------------------------
350 params['c0']['value'] concentration of first injection, (s)
351 ... ...
352 params['c7']['value'] concentration of eighth injection, (s)
353 ======================= ==========================================
354
355 Injections where tN=0 are ignored. For example, if you have only
356 six injections, set t6=0 and t7=0.
357 """
358
359 ## Parameters from params dictionary.
360 ## Skip parameter validation steps for now.
361 kon = float(params['kon']['value'])
362 koff = float(params['koff']['value'])
363 rmax = float(params['rmax']['value'])
364 ## Pre-association, Association duration, Dissociation duration.
365 tp = float(params['tp']['value']) ## Typical -100 s
366 ta = float(params['ta']['value']) ## Typical 300 s
367 td = float(params['td']['value']) ## Typical 300 s
368 td = td + ta ## After 300s binding and 300s washing, 600s elapse.
369 ## How many injections are there? Check valid keys and inj times.
370 tkey = ['t0','t1','t2','t3','t4','t5','t6','t7']
371 ckey = ['c0','c1','c2','c3','c4','c5','c6','c7']
372 for i in range(8):
373 if (tkey[i] in params) and (ckey[i] in params):
374 ## Valid inj time is after sufficient pre-association time.
375 if ((float(params[tkey[i]]['value']) + tp) > 0):
376 inj = i+1 ## Increase inj count.
377 ## Up to eight injection start times, concentrations.
378 ti = [float(params[k]['value']) for k in tkey[0:inj]]
379 c = [float(params[k]['value']) for k in ckey[0:inj]]
380
381 ## Begin.
382 y = np.zeros(len(time), dtype=float)
383 base = y[0] = data[0] ## Baseline SPR signal.
384
385 ## Must iterate through data, numerical integration.
386 for i in range(1,len(time)):
387 dt = time[i] - time[i-1]
388 #if (time[i] > (max(ti)+td)): break ## Speed things up.
389 inmodel = False
390 ## For-loop where j is each injection number.
391 ## Treat function as having three separate phases.
392 for j in range(inj):
393 if ((ti[j]+tp) <= time[i] <= (ti[j])):
394 ## Pre-binding phase.
395 if np.isnan(y[i-1]): y[i-1]=0
396 y[i] = (0.6 * y[i-1]) + (0.4 * data[i])
397 base = y[i]
398 inmodel = True
399 if ((ti[j]) <= time[i] <= (ti[j]+ta)):
400 ## Binding phase
401 yb = y[i-1] - base
402 dy = c[j]*kon*(rmax-yb) - koff*yb
403 y[i] = y[i-1] + dy*dt
404 inmodel = True
405 if (y[i] > (base+rmax)): y[i] = (base+rmax) ## Avoid overshoot.
406 if ((ti[j]+ta) <= time[i] <= (ti[j]+td)):
407 ## Dissociation (conc=0)
408 yb = y[i-1]-base
409 dy = 0 - koff*yb
410 y[i] = y[i-1] + dy*dt
411 inmodel = True
412 ## End "for j" loop.
413 if (inmodel==False):
414 ## Not modeled, e.g. regeneration.
415 y[i] = np.NaN
416 y[i] = 0 ## NaN might break mclma()?
417 ## End "for i" loop.
418 return y
419
420 def simple1to1_series_def_params():
421 """
422 Return a dictionary of default parameters for simple1to1_series model.
423 While parameters for eight injection times are provided, only the
424 first injection time is nonzero. For example::
425
426 >>> from numpy import arange, zeros
427 >>> defaultparam0 = simple1to1_series_def_params()
428 >>> data0 = simple1to1_series(arange(300), zeros(300), defaultparam0)
429 >>> print round(data0[298], 1)
430 160.3
431 """
432 dpar = {}
433 dpar['kon'] = dict(value=1e4, min=1e2, max=1e7, fixed='fixed')
434 dpar['koff'] = dict(value=1e-2, min=1e-6, max=1e-1, fixed='fixed')
435 dpar['rmax'] = dict(value=1e3, min=1e1, max=1e4, fixed='fixed')
436 dpar['tp'] = dict(value=-50.0, min=-50.0, max=-50.0, fixed='fixed')
437 dpar['ta'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
438 dpar['td'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
439 ## Set all eight concentrations to 1uM.
440 for key in ['c0','c1','c2','c3','c4','c5','c6','c7']:
441 dpar[key] = dict(value=1e-6, min=1e-12, max=1e-1, fixed='fixed')
442 ## Just describe the first injection at t0=100s.
443 for key in ['t0','t1','t2','t3','t4','t5','t6','t7']:
444 dpar[key] = dict(value=0.0, min=0.0, max=0.0, fixed='fixed')
445 dpar['t0'] = dict(value=100.0, min=100.0, max=100.0, fixed='fixed')
446 return dpar
447
448
449 def _test():
450 """
451 Automatic Code testing with doctest.
452 Enable use of ellipses as wildcard for returned text.
453 Doctest runs the example code in the docstrings as well as the following.
454 We test that diff eq models here match the simulation in BioSensorArray.
455 Here are the tests::
456
457 >>> from numpy import arange, zeros
458 >>> dpar = simple1to1_series_def_params()
459 >>> values = simple1to1_series(arange(300), zeros(300), dpar)
460 >>> print round(values[298], -1) ## Round to 10s place.
461 160.0
462 >>> import ba_class as ba
463 >>> ba1 = ba.BiosensorArray(1, 300)
464 >>> ba1.roi[0].simulate(bulkRI=0.0, kon=dpar['kon']['value'], \
465 koff=dpar['koff']['value'], conc=dpar['c0']['value'])
466 >>> print round(ba1.roi[0].value[299], -1)
467 160.0
468 """
469 import doctest
470 doctest.testmod(optionflags=doctest.ELLIPSIS)
471 return
472
473 if __name__ == "__main__":
474 """Use 'python mdl_module.py' for code testing."""
475 _test()
476
477 ################################# End of module ###############################