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 100518 (yymmdd)
6 <
7 < Examples:
8 < #import mdl_module as mdl
9 < #import numpy as np
10 < #times = np.arange(100)
11 < #data = np.zeros(100)
12 < #
13 < #params1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
14 < #data1 = drift(time, data, param1)
15 < #
16 < #param2 = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
17 < #param2['rmax'] =  {'value': 100.0}
18 < #param2['conc'] =  {'value': 1e-6}
19 < #param2['kon'] =   {'value': 2e4}
20 < #param2['t2'] =    {'value': 150.0}
21 < #param2['koff'] =  {'value': 1e-3}
22 < #param2['t3'] =    {'value': 270.0}
23 < #data2 = simple1to1(time, data, param2)
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 = np.arange(100)
45 >  >>> data = np.zeros(100)
46 >  >>>
47 >  >>> param1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
48 >  >>> data1 = drift(time, data, param1)
49 >  >>>
50 >  >>> param2 = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
51 >  >>> param2['rmax'] =  {'value': 100.0}
52 >  >>> param2['conc'] =  {'value': 1e-6}
53 >  >>> param2['kon'] =   {'value': 2e4}
54 >  >>> param2['t2'] =    {'value': 150.0}
55 >  >>> param2['koff'] =  {'value': 1e-3}
56 >  >>> param2['t3'] =    {'value': 270.0}
57 >  >>> data2 = simple1to1(time, data, param2)
58   """
59 < __version__ = "100518"
59 > __version__ = "110216"
60  
61  
62   ## Import libraries
# Line 35 | Line 69
69      This function simply models a constant signal drift in units/second.
70      It requires numpy arrays of times and starting data values,
71      It only requires one parameter in the params list.
72 <    params['rate']['value']
72 >    
73 >    ``params['rate']['value']``
74      """
75      y = np.zeros(len(time), dtype=float)
76      try:
# Line 57 | Line 92
92   def simple1to1(time, data, params):
93      """
94      This function simply models a 1:1 interaction.
95 <    [A] + [L] --kon--> [AL]
96 <    [A] + [L] <-koff-- [AL]
95 >    The model parameters are described in the table below.
96 >    
97 >    Model::
98 >    
99 >      [A] + [L] --kon--> [AL]
100 >      [A] + [L] <-koff-- [AL]
101 >    
102 >    Derivation::
103      
104 <    It requires numpy arrays of times and starting data values.
105 <    params['t1']['value']   is time of injection for binding, (s)
106 <    params['rmax']['value'] is maximum response, (RIU)
107 <    params['conc']['value'] is stock concentration of analyte [A], (M)
108 <    params['kon']['value']  is on-rate of analyte, (1/Ms)
109 <    params['t2']['value']   is time of end binding & begin washing, (s)
110 <    params['koff']['value'] is off-rate of analyte, (1/s)
111 <    params['t3']['value']   is time end of washing & data fitting, (s)
104 >      d[AL]/dt = kon*[A]*[L] - koff*[AL]
105 >      y = [AL]
106 >      (rmax-y) = [L]
107 >      conc = [A]
108 >      rmax = [AL] + [L]
109 >      dy/dt = conc*kon*(rmax-y) - koff*y
110 >    
111 >    ======================= ============================================
112 >    Model Parameter         Description
113 >    ======================= ============================================
114 >    params['t1']['value']   time of injection for binding, (s)
115 >    params['rmax']['value'] maximum response, (RIU)
116 >    params['conc']['value'] concentration of analyte [A], (M)
117 >    params['kon']['value']  on-rate of analyte, (1/Ms)
118 >    params['t2']['value']   time of end binding & begin washing, (s)
119 >    params['koff']['value'] off-rate of analyte, (1/s)
120 >    params['t3']['value']   time end of washing & data fitting, (s)
121 >    ----------------------- --------------------------------------------
122 >    *Optional*            
123 >    params['cofa']['value'] concentration factor, (1/dilution factor)
124 >    ======================= ============================================
125      
72    It also can take an optional parameter useful for dilution series.
73    params['cofa']['value'] is concentration factor, (1/dilution factor)
126      """
127 +    
128      ## Skip parameter validation steps for now.
129      t1 = float(params['t1']['value'])
130      rmax = float(params['rmax']['value'])
# Line 82 | Line 135
135      t3 = float(params['t3']['value'])
136      if ('cofa' in params.keys()):
137          conc *= float(params['cofa']['value'])
138 <    
139 <    """
140 <    Derivation:
88 <    d[AL]/dt = kon*[A]*[L] - koff*[AL]
89 <    y = [AL]
90 <    (rmax-y) = [L]
91 <    conc = [A]
92 <    rmax = [AL] + [L]
93 <    dy/dt = conc*kon*(rmax-y) - koff*y
94 <    """
95 <    
138 >        
139 >    ## Initialize variables.
140 >    errlog = {}  ## Error status dictionary.
141      y = np.zeros(len(time), dtype=float)
142      base = y[0] = data[0]  ## Baseline SPR signal.
143      
# Line 107 | Line 152
152              ## Binding phase
153              yb = y[i-1] - base
154              dy = conc*kon*(rmax-yb) - koff*yb
110            if (abs(y[i-1]) > 999999999): dy = 0  ## Is this useful?
155              y[i] = y[i-1] + dy*dt
156 +            ## Check if we overshot the maximum response.
157 +            if (y[i] > (base+rmax)):
158 +                y[i] = base+rmax
159 +                errlog['Response overshot rmax.'] = True
160          elif (t2 < time[i] <= t3):
161              ## Dissociation (conc=0)
162              yb = y[i-1]-base
163              dy = 0 - koff*yb
164              y[i] = y[i-1] + dy*dt
165 +    
166 +    if any(errlog): print "Errors in simple1to1:", errlog.keys()
167              
168      return y
169 <    ## End of simple1to1() function
169 >    """End of simple1to1() function"""
170  
171  
172   def simple1to1_mtl(time, data, params):
173      """
174      This function simply models a 1:1 interaction with mass transport limitation.
175 <    [Abulk]  --km->  [Asurf] + [L]  --kon-->  [AL]
126 <    [Abulk]  <-km--  [Asurf] + [L]  <-koff--  [AL]
175 >    The model parameters are described in the table below.
176      
177 <    It requires numpy arrays of times and starting data values,
178 <    params['t1']['value']   is time of injection for binding, (s)
179 <    params['rmax']['value'] is maximum response, (RIU)
180 <    params['conc']['value'] is concentration of analyte [Abulk], (M)
181 <    params['kon']['value']  is on-rate of analyte, (1/Ms)
182 <    params['t2']['value']   is time of end binding & begin washing, (s)
183 <    params['koff']['value'] is off-rate of analyte, (1/s)
184 <    params['t3']['value']   is time end of washing & data fitting, (s)
185 <    params['kmtl']['value'] is rate of diffusion,  (RIU/Ms)
177 >    Model::
178 >    
179 >      [Abulk]  --km->  [Asurf] + [L]  --kon-->  [AL]
180 >      [Abulk]  <-km--  [Asurf] + [L]  <-koff--  [AL]
181 >    
182 >    Derivation::
183 >    
184 >      d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
185 >      y = [AL]
186 >      (rmax-y) = [L]
187 >      conc = [Abulk]
188 >      rmax = [AL] + [L]
189 >      dy/dt  = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
190 >    
191 >    ======================= ============================================
192 >    Model Parameter         Description
193 >    ======================= ============================================
194 >    params['t1']['value']   time of injection for binding, (s)
195 >    params['rmax']['value'] maximum response, (RIU)
196 >    params['conc']['value'] concentration of analyte [Abulk], (M)
197 >    params['kon']['value']  on-rate of analyte, (1/Ms)
198 >    params['t2']['value']   time of end binding & begin washing, (s)
199 >    params['koff']['value'] off-rate of analyte, (1/s)
200 >    params['t3']['value']   time end of washing & data fitting, (s)
201 >    params['kmtl']['value'] rate of diffusion,  (RIU/Ms)
202 >    ----------------------- --------------------------------------------
203 >    *Optional*            
204 >    params['cofa']['value'] concentration factor, (1/dilution factor)
205 >    ======================= ============================================
206      
138    It also can take an optional parameter useful for dilution series.
139    params['cofa']['value'] is concentration factor, (1/dilution factor)
207      """
208 +    
209      ## Skip parameter validation steps for now.
210      t1 = params['t1']['value']
211      rmax = params['rmax']['value']
# Line 150 | Line 218
218      if ('cofa' in params.keys()):
219          conc *= float(params['cofa']['value'])
220          
221 <    """
222 <    Derivation:
155 <    d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
156 <    y = [AL]
157 <    (rmax-y) = [L]
158 <    conc = [Abulk]
159 <    rmax = [AL] + [L]
160 <    dy/dt  = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
161 <    """
162 <    
163 <    stat = {}  ## Error status dictionary.
221 >    ## Initialize variables.
222 >    errlog = {}  ## Error status dictionary.
223      y = np.zeros(len(time), dtype=float)
224      base = y[0] = data[0]  ## Baseline SPR signal.
225      
# Line 181 | Line 240
240              ## Binding phase
241              yb = y[i-1] - base
242              dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
184            ## Check for integration errors producing wild results.
185            if (abs(y[i-1]) > 9e9): dy = 0; stat['y>9e9'] = True
243              y[i] = y[i-1] + dy*dt
244 +            ## Check if we overshot the maximum response.
245 +            if (y[i] > (base+rmax)):
246 +                y[i] = base+rmax
247 +                errlog['Response overshot rmax.'] = True
248          elif (t2 < time[i] <= t3):
249              ## Dissociation (conc=0)
250              yb = y[i-1] - base
251              dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
252              y[i] = y[i-1] + dy*dt
253      
254 <    if (len(stat.keys()) > 0): print "Errors in simple1to1:", stat.keys()
254 >    if any(errlog): print "Errors in simple1to1:", errlog.keys()
255      
256      return y
257 <    ## End of simple1to1_mtl() function
257 >    """End of simple1to1_mtl() function"""
258  
259  
260  
200 ################################# End of module #################################
261 + ################################# End of module ###############################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines