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, 5 months ago) by clausted
File size: 9594 byte(s)
Log Message:
Improve mdl_module to avoid SPR responses that overshoot Rmax.
Line File contents
1 """
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__ = "110216"
60
61
62 ## Import libraries
63 import numpy as np
64 from copy import deepcopy
65
66
67 def drift(time, data, params):
68 """
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
73 ``params['rate']['value']``
74 """
75 y = np.zeros(len(time), dtype=float)
76 try:
77 rate = params['rate']['value']
78 except KeyError:
79 print "Error: Rate parameter not supplied to drift() model."
80 return y
81
82 ## 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 ## End of drift() function.
90
91
92 def simple1to1(time, data, params):
93 """
94 This function simply models a 1:1 interaction.
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 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 """
127
128 ## Skip parameter validation steps for now.
129 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
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
144 ## Must iterate through data, numerical integration.
145 for i in range(1,len(time)):
146 dt = time[i] - time[i-1]
147 ## Treat function as having three separate phases.
148 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 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"""
170
171
172 def simple1to1_mtl(time, data, params):
173 """
174 This function simply models a 1:1 interaction with mass transport limitation.
175 The model parameters are described in the table below.
176
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
207 """
208
209 ## 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 if ('cofa' in params.keys()):
219 conc *= float(params['cofa']['value'])
220
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
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 ## 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 any(errlog): print "Errors in simple1to1:", errlog.keys()
255
256 return y
257 """End of simple1to1_mtl() function"""
258
259
260
261 ################################# End of module ###############################