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 (9 years, 2 months 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 User Rev File contents
1 clausted 17 """
2 clausted 53 mdl_module
3     -------------------------------------------------------------------------------
4 clausted 17
5 clausted 54 *Example model functions module for SPRI data.*
6 clausted 53 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 clausted 54 Each model is a function taking three parameters.
9 clausted 53
10 clausted 54 **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 clausted 53 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 clausted 54 * ``'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 clausted 53
26 clausted 54 The ``min``, ``max``, and ``fixed`` keys are used during automatic curve-fitting.
27 clausted 53 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 clausted 60 >>> times, data = np.arange(300), np.zeros(300)
45 clausted 53 >>>
46 clausted 60 >>> 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 clausted 53 >>>
55 clausted 60 >>> 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 clausted 17 """
67 clausted 60 __version__ = "110517"
68 clausted 17
69    
70     ## Import libraries
71     import numpy as np
72 clausted 18 from copy import deepcopy
73 clausted 17
74    
75 clausted 18 def drift(time, data, params):
76 clausted 17 """
77     This function simply models a constant signal drift in units/second.
78 clausted 18 It requires numpy arrays of times and starting data values,
79 clausted 17 It only requires one parameter in the params list.
80 clausted 53
81 clausted 54 ``params['rate']['value']``
82 clausted 17 """
83 clausted 18 y = np.zeros(len(time), dtype=float)
84 clausted 17 try:
85 clausted 18 rate = params['rate']['value']
86 clausted 17 except KeyError:
87     print "Error: Rate parameter not supplied to drift() model."
88 clausted 18 return y
89 clausted 17
90 clausted 18 ## 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 clausted 19 ## End of drift() function.
98 clausted 18
99 clausted 60 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 clausted 18
112 clausted 60
113 clausted 18 def simple1to1(time, data, params):
114     """
115 clausted 26 This function simply models a 1:1 interaction.
116 clausted 54 The model parameters are described in the table below.
117 clausted 26
118 clausted 53 Model::
119 clausted 29
120 clausted 53 [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 clausted 18 """
148 clausted 53
149 clausted 18 ## Skip parameter validation steps for now.
150 clausted 29 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 clausted 55
160     ## Initialize variables.
161     errlog = {} ## Error status dictionary.
162 clausted 18 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 clausted 26 dt = time[i] - time[i-1]
168     ## Treat function as having three separate phases.
169 clausted 18 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 clausted 26 y[i] = y[i-1] + dy*dt
177 clausted 55 ## 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 clausted 18 elif (t2 < time[i] <= t3):
182 clausted 26 ## Dissociation (conc=0)
183 clausted 18 yb = y[i-1]-base
184     dy = 0 - koff*yb
185 clausted 26 y[i] = y[i-1] + dy*dt
186 clausted 55
187     if any(errlog): print "Errors in simple1to1:", errlog.keys()
188 clausted 18
189     return y
190 clausted 53 """End of simple1to1() function"""
191 clausted 18
192 clausted 60 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 clausted 19
213 clausted 60
214 clausted 26 def simple1to1_mtl(time, data, params):
215     """
216     This function simply models a 1:1 interaction with mass transport limitation.
217 clausted 54 The model parameters are described in the table below.
218 clausted 26
219 clausted 53 Model::
220 clausted 29
221 clausted 53 [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 clausted 26 """
250 clausted 53
251 clausted 26 ## 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 clausted 29 if ('cofa' in params.keys()):
261     conc *= float(params['cofa']['value'])
262    
263 clausted 53 ## Initialize variables.
264 clausted 55 errlog = {} ## Error status dictionary.
265 clausted 26 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 clausted 55 ## 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 clausted 26 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 clausted 55 if any(errlog): print "Errors in simple1to1:", errlog.keys()
297 clausted 26
298     return y
299 clausted 53 """End of simple1to1_mtl() function"""
300 clausted 60
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 clausted 26
316    
317 clausted 60 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 clausted 26
433 clausted 60
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 clausted 53 ################################# End of module ###############################