ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
(Generate patch)
# Line 41 | Line 41
41  
42    >>> import mdl_module as mdl
43    >>> import numpy as np
44 <  >>> times = np.arange(100)
45 <  >>> data = np.zeros(100)
44 >  >>> times, data = np.arange(300), np.zeros(300)
45    >>>
46 <  >>> param1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
47 <  >>> data1 = drift(time, data, param1)
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 = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
56 <  >>> param2['rmax'] =  {'value': 100.0}
57 <  >>> param2['conc'] =  {'value': 1e-6}
58 <  >>> param2['kon'] =   {'value': 2e4}
59 <  >>> param2['t2'] =    {'value': 150.0}
60 <  >>> param2['koff'] =  {'value': 1e-3}
61 <  >>> param2['t3'] =    {'value': 270.0}
62 <  >>> data2 = simple1to1(time, data, param2)
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__ = "110216"
67 > __version__ = "110517"
68  
69  
70   ## Import libraries
# Line 88 | Line 96
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      """
# Line 168 | Line 189
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      """
# Line 255 | Line 297
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 ###############################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines