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