ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
(Generate patch)
# Line 1 | Line 1
1   """
2 < mdl: Example model functions module for SPRI data.
3 < Christopher Lausted, Institute for Systems Biology,
4 < OSPRAI developers
5 < Last modified on 100421 (yymmdd)
6 <
7 < Example:
8 < #import mdl_module as mdl
9 < #params = {'Rate': {'value':1, 'min':-100.0, 'max':100.0, 'fixed':False} }
10 < #params = dict(Rate = dict(value=1, min=0, max=10, fixed=False))
11 < #times = range(100)
12 < #data1 = drift(time, data0, params)
13 < print data1
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__ = "100421"
67 > __version__ = "110517"
68  
69  
70   ## Import libraries
# Line 25 | Line 77
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 <    params['rate']['value']
80 >    
81 >    ``params['rate']['value']``
82      """
83      y = np.zeros(len(time), dtype=float)
84      try:
# Line 41 | Line 94
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 <    It requires numpy arrays of times and starting data values,
117 <    params['t1']['value']   is time of injection for binding
118 <    params['rmax']['value'] is maximum response.
119 <    params['conc']['value'] is time of concentration of analyte
120 <    params['kon']['value']  is on-rate of analyte
121 <    params['t2']['value']   is time of end binding / begin washing
122 <    params['koff']['value']  is off-rate of analyte
123 <    params['t3']['value']   is time end of washing / data fitting.
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']
# Line 62 | Line 255
255      kon = params['kon']['value']
256      t2 = params['t2']['value']
257      koff = params['koff']['value']
258 <    t3 = params['t3']['value']
259 <    
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
285 <            if (abs(y[i-1]) > 999999999): dy = 0  ## Is this useful?
286 <            y[i] = y[i-1] + dy * (time[i] - time[i-1])
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
292 <            yb = y[i-1]-base
293 <            dy = 0 - koff*yb
294 <            y[i] = y[i-1] + dy * (time[i] - time[i-1])
295 <            
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 ###############################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines