ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 54
Committed: Wed Feb 16 21:29:02 2011 UTC (8 years, 7 months ago) by clausted
File size: 9274 byte(s)
Log Message:
Added fit_module and mdl_module to sphinx documentation.
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 y = np.zeros(len(time), dtype=float)
140 base = y[0] = data[0] ## Baseline SPR signal.
141
142 ## Must iterate through data, numerical integration.
143 for i in range(1,len(time)):
144 dt = time[i] - time[i-1]
145 ## Treat function as having three separate phases.
146 if (time[i] <= t1):
147 ## Pre-binding phase.
148 base = y[i] = data[i]
149 elif (t1 < time[i] <= t2):
150 ## Binding phase
151 yb = y[i-1] - base
152 dy = conc*kon*(rmax-yb) - koff*yb
153 if (abs(y[i-1]) > 999999999): dy = 0 ## Is this useful?
154 y[i] = y[i-1] + dy*dt
155 elif (t2 < time[i] <= t3):
156 ## Dissociation (conc=0)
157 yb = y[i-1]-base
158 dy = 0 - koff*yb
159 y[i] = y[i-1] + dy*dt
160
161 return y
162 """End of simple1to1() function"""
163
164
165 def simple1to1_mtl(time, data, params):
166 """
167 This function simply models a 1:1 interaction with mass transport limitation.
168 The model parameters are described in the table below.
169
170 Model::
171
172 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
173 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
174
175 Derivation::
176
177 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
178 y = [AL]
179 (rmax-y) = [L]
180 conc = [Abulk]
181 rmax = [AL] + [L]
182 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
183
184 ======================= ============================================
185 Model Parameter Description
186 ======================= ============================================
187 params['t1']['value'] time of injection for binding, (s)
188 params['rmax']['value'] maximum response, (RIU)
189 params['conc']['value'] concentration of analyte [Abulk], (M)
190 params['kon']['value'] on-rate of analyte, (1/Ms)
191 params['t2']['value'] time of end binding & begin washing, (s)
192 params['koff']['value'] off-rate of analyte, (1/s)
193 params['t3']['value'] time end of washing & data fitting, (s)
194 params['kmtl']['value'] rate of diffusion, (RIU/Ms)
195 ----------------------- --------------------------------------------
196 *Optional*
197 params['cofa']['value'] concentration factor, (1/dilution factor)
198 ======================= ============================================
199
200 """
201
202 ## Skip parameter validation steps for now.
203 t1 = params['t1']['value']
204 rmax = params['rmax']['value']
205 conc = params['conc']['value']
206 kon = params['kon']['value']
207 t2 = params['t2']['value']
208 koff = params['koff']['value']
209 t3 = params['t3']['value']
210 kmtl = params['kmtl']['value']
211 if ('cofa' in params.keys()):
212 conc *= float(params['cofa']['value'])
213
214 ## Initialize variables.
215 stat = {} ## Error status dictionary.
216 y = np.zeros(len(time), dtype=float)
217 base = y[0] = data[0] ## Baseline SPR signal.
218
219 ## Error checks.
220 if (kmtl<=0):
221 ## Avoid div/0, assume user wants no mass transport limitation.
222 stat['kmtl must be > 0'] = True
223 kmtl = 1e40 ## An arbitrary very large number.
224
225 ## Must iterate through data, numerical integration.
226 for i in range(1,len(time)):
227 dt = time[i] - time[i-1]
228 ## Treat function as having three separate phases.
229 if (time[i] <= t1):
230 ## Pre-binding phase.
231 base = y[i] = data[i]
232 elif (t1 < time[i] <= t2):
233 ## Binding phase
234 yb = y[i-1] - base
235 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
236 ## Check for integration errors producing wild results.
237 if (abs(y[i-1]) > 9e9): dy = 0; stat['y>9e9'] = True
238 y[i] = y[i-1] + dy*dt
239 elif (t2 < time[i] <= t3):
240 ## Dissociation (conc=0)
241 yb = y[i-1] - base
242 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
243 y[i] = y[i-1] + dy*dt
244
245 if (len(stat.keys()) > 0): print "Errors in simple1to1:", stat.keys()
246
247 return y
248 """End of simple1to1_mtl() function"""
249
250
251
252 ################################# End of module ###############################