ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 55
Committed: Tue Feb 22 00:27:41 2011 UTC (8 years, 7 months ago) by clausted
File size: 9594 byte(s)
Log Message:
Improve mdl_module to avoid SPR responses that overshoot Rmax.
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     >>> times = np.arange(100)
45     >>> data = np.zeros(100)
46     >>>
47 clausted 54 >>> param1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
48 clausted 53 >>> 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 clausted 17 """
59 clausted 54 __version__ = "110216"
60 clausted 17
61    
62     ## Import libraries
63     import numpy as np
64 clausted 18 from copy import deepcopy
65 clausted 17
66    
67 clausted 18 def drift(time, data, params):
68 clausted 17 """
69     This function simply models a constant signal drift in units/second.
70 clausted 18 It requires numpy arrays of times and starting data values,
71 clausted 17 It only requires one parameter in the params list.
72 clausted 53
73 clausted 54 ``params['rate']['value']``
74 clausted 17 """
75 clausted 18 y = np.zeros(len(time), dtype=float)
76 clausted 17 try:
77 clausted 18 rate = params['rate']['value']
78 clausted 17 except KeyError:
79     print "Error: Rate parameter not supplied to drift() model."
80 clausted 18 return y
81 clausted 17
82 clausted 18 ## Numerical integration.
83     y[0] = data[0]
84     for i in range(1,len(time)):
85     dy = rate * (time[i]-time[i-1])
86     y[i] = y[i-1] + dy
87    
88     return y
89 clausted 19 ## End of drift() function.
90 clausted 18
91    
92     def simple1to1(time, data, params):
93     """
94 clausted 26 This function simply models a 1:1 interaction.
95 clausted 54 The model parameters are described in the table below.
96 clausted 26
97 clausted 53 Model::
98 clausted 29
99 clausted 53 [A] + [L] --kon--> [AL]
100     [A] + [L] <-koff-- [AL]
101    
102     Derivation::
103    
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    
126 clausted 18 """
127 clausted 53
128 clausted 18 ## Skip parameter validation steps for now.
129 clausted 29 t1 = float(params['t1']['value'])
130     rmax = float(params['rmax']['value'])
131     conc = float(params['conc']['value'])
132     kon = float(params['kon']['value'])
133     t2 = float(params['t2']['value'])
134     koff = float(params['koff']['value'])
135     t3 = float(params['t3']['value'])
136     if ('cofa' in params.keys()):
137     conc *= float(params['cofa']['value'])
138 clausted 55
139     ## Initialize variables.
140     errlog = {} ## Error status dictionary.
141 clausted 18 y = np.zeros(len(time), dtype=float)
142     base = y[0] = data[0] ## Baseline SPR signal.
143    
144     ## Must iterate through data, numerical integration.
145     for i in range(1,len(time)):
146 clausted 26 dt = time[i] - time[i-1]
147     ## Treat function as having three separate phases.
148 clausted 18 if (time[i] <= t1):
149     ## Pre-binding phase.
150     base = y[i] = data[i]
151     elif (t1 < time[i] <= t2):
152     ## Binding phase
153     yb = y[i-1] - base
154     dy = conc*kon*(rmax-yb) - koff*yb
155 clausted 26 y[i] = y[i-1] + dy*dt
156 clausted 55 ## 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 clausted 18 elif (t2 < time[i] <= t3):
161 clausted 26 ## Dissociation (conc=0)
162 clausted 18 yb = y[i-1]-base
163     dy = 0 - koff*yb
164 clausted 26 y[i] = y[i-1] + dy*dt
165 clausted 55
166     if any(errlog): print "Errors in simple1to1:", errlog.keys()
167 clausted 18
168     return y
169 clausted 53 """End of simple1to1() function"""
170 clausted 18
171 clausted 19
172 clausted 26 def simple1to1_mtl(time, data, params):
173     """
174     This function simply models a 1:1 interaction with mass transport limitation.
175 clausted 54 The model parameters are described in the table below.
176 clausted 26
177 clausted 53 Model::
178 clausted 29
179 clausted 53 [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    
207 clausted 26 """
208 clausted 53
209 clausted 26 ## Skip parameter validation steps for now.
210     t1 = params['t1']['value']
211     rmax = params['rmax']['value']
212     conc = params['conc']['value']
213     kon = params['kon']['value']
214     t2 = params['t2']['value']
215     koff = params['koff']['value']
216     t3 = params['t3']['value']
217     kmtl = params['kmtl']['value']
218 clausted 29 if ('cofa' in params.keys()):
219     conc *= float(params['cofa']['value'])
220    
221 clausted 53 ## Initialize variables.
222 clausted 55 errlog = {} ## Error status dictionary.
223 clausted 26 y = np.zeros(len(time), dtype=float)
224     base = y[0] = data[0] ## Baseline SPR signal.
225    
226     ## Error checks.
227     if (kmtl<=0):
228     ## Avoid div/0, assume user wants no mass transport limitation.
229     stat['kmtl must be > 0'] = True
230     kmtl = 1e40 ## An arbitrary very large number.
231    
232     ## Must iterate through data, numerical integration.
233     for i in range(1,len(time)):
234     dt = time[i] - time[i-1]
235     ## Treat function as having three separate phases.
236     if (time[i] <= t1):
237     ## Pre-binding phase.
238     base = y[i] = data[i]
239     elif (t1 < time[i] <= t2):
240     ## Binding phase
241     yb = y[i-1] - base
242     dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
243     y[i] = y[i-1] + dy*dt
244 clausted 55 ## 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 clausted 26 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 clausted 55 if any(errlog): print "Errors in simple1to1:", errlog.keys()
255 clausted 26
256     return y
257 clausted 53 """End of simple1to1_mtl() function"""
258 clausted 26
259    
260    
261 clausted 53 ################################# End of module ###############################