ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/mdl_module.py
Revision: 53
Committed: Wed Feb 16 01:56:39 2011 UTC (8 years, 4 months ago) by clausted
File size: 9288 byte(s)
Log Message:
Add mdl_module documentation.  Also add sphinx html output to svn repository.
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 1. *time* An array of time values, usually in seconds.
11 2. *data* An array of initial SPR signal values.
12 3. *params* A dictionary of dictionaries.
13
14 The dictionary keys are the names of model parameters
15 (e.g. 'rmax' for maximal response, or 'kon' for kinetic on-rate).
16 Each model parameter is described by a subdictionary containing four entries.
17
18 1. *'value'* The current value of the parameter.
19 2. *'min'* The minimum allowable value.
20 3. *'max'* The maximum allowable value.
21 4. *'fixed'* Either 'float', 'fixed', or a reference to another ROI.
22
23 The *min*, *max*, and *fixed* keys are used during automatic curve-fitting.
24 A fixed parameter is not allowed to change, while a float parameter is adjusted
25 until the least-squares algorithm has minimized the sum-squared error.
26 The *fixed* parameter may also be an integer, in which case it is fixed to
27 the value of a parameter of the same name in another ROI.
28
29 The model function returns an array of values obtained by numerical integration.
30 The model is represented by differential equations and integrated using the
31 rectangle rule or, preferentially, using the trapezoidal rule.
32
33 .. moduleauthor:: Christopher Lausted,
34 Institute for Systems Biology,
35 OSPRAI developers.
36
37 Examples::
38
39 >>> import mdl_module as mdl
40 >>> import numpy as np
41 >>> times = np.arange(100)
42 >>> data = np.zeros(100)
43 >>>
44 >>> params1 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
45 >>> data1 = drift(time, data, param1)
46 >>>
47 >>> param2 = {'t1': {'value':30.0, 'min':30.0, 'max':30.0, 'fixed':True} }
48 >>> param2['rmax'] = {'value': 100.0}
49 >>> param2['conc'] = {'value': 1e-6}
50 >>> param2['kon'] = {'value': 2e4}
51 >>> param2['t2'] = {'value': 150.0}
52 >>> param2['koff'] = {'value': 1e-3}
53 >>> param2['t3'] = {'value': 270.0}
54 >>> data2 = simple1to1(time, data, param2)
55 """
56 __version__ = "110215"
57
58
59 ## Import libraries
60 import numpy as np
61 from copy import deepcopy
62
63
64 def drift(time, data, params):
65 """
66 This function simply models a constant signal drift in units/second.
67 It requires numpy arrays of times and starting data values,
68 It only requires one parameter in the params list.
69
70 params['rate']['value']
71
72 :param time: Time points.
73 :type time: numpy array
74 :param data: SPR signal values.
75 :type data: numpy array
76 :param params: Binding model parameter description.
77 :type params: dictionary of dictionaries
78 :returns: numpy array
79 """
80 y = np.zeros(len(time), dtype=float)
81 try:
82 rate = params['rate']['value']
83 except KeyError:
84 print "Error: Rate parameter not supplied to drift() model."
85 return y
86
87 ## Numerical integration.
88 y[0] = data[0]
89 for i in range(1,len(time)):
90 dy = rate * (time[i]-time[i-1])
91 y[i] = y[i-1] + dy
92
93 return y
94 ## End of drift() function.
95
96
97 def simple1to1(time, data, params):
98 """
99 This function simply models a 1:1 interaction.
100
101 Model::
102
103 [A] + [L] --kon--> [AL]
104 [A] + [L] <-koff-- [AL]
105
106 Derivation::
107
108 d[AL]/dt = kon*[A]*[L] - koff*[AL]
109 y = [AL]
110 (rmax-y) = [L]
111 conc = [A]
112 rmax = [AL] + [L]
113 dy/dt = conc*kon*(rmax-y) - koff*y
114
115 ======================= ============================================
116 Model Parameter Description
117 ======================= ============================================
118 params['t1']['value'] time of injection for binding, (s)
119 params['rmax']['value'] maximum response, (RIU)
120 params['conc']['value'] concentration of analyte [A], (M)
121 params['kon']['value'] on-rate of analyte, (1/Ms)
122 params['t2']['value'] time of end binding & begin washing, (s)
123 params['koff']['value'] off-rate of analyte, (1/s)
124 params['t3']['value'] time end of washing & data fitting, (s)
125 ----------------------- --------------------------------------------
126 *Optional*
127 params['cofa']['value'] concentration factor, (1/dilution factor)
128 ======================= ============================================
129
130 """
131
132 ## Skip parameter validation steps for now.
133 t1 = float(params['t1']['value'])
134 rmax = float(params['rmax']['value'])
135 conc = float(params['conc']['value'])
136 kon = float(params['kon']['value'])
137 t2 = float(params['t2']['value'])
138 koff = float(params['koff']['value'])
139 t3 = float(params['t3']['value'])
140 if ('cofa' in params.keys()):
141 conc *= float(params['cofa']['value'])
142
143 y = np.zeros(len(time), dtype=float)
144 base = y[0] = data[0] ## Baseline SPR signal.
145
146 ## Must iterate through data, numerical integration.
147 for i in range(1,len(time)):
148 dt = time[i] - time[i-1]
149 ## Treat function as having three separate phases.
150 if (time[i] <= t1):
151 ## Pre-binding phase.
152 base = y[i] = data[i]
153 elif (t1 < time[i] <= t2):
154 ## Binding phase
155 yb = y[i-1] - base
156 dy = conc*kon*(rmax-yb) - koff*yb
157 if (abs(y[i-1]) > 999999999): dy = 0 ## Is this useful?
158 y[i] = y[i-1] + dy*dt
159 elif (t2 < time[i] <= t3):
160 ## Dissociation (conc=0)
161 yb = y[i-1]-base
162 dy = 0 - koff*yb
163 y[i] = y[i-1] + dy*dt
164
165 return y
166 """End of simple1to1() function"""
167
168
169 def simple1to1_mtl(time, data, params):
170 """
171 This function simply models a 1:1 interaction with mass transport limitation.
172
173 Model::
174
175 [Abulk] --km-> [Asurf] + [L] --kon--> [AL]
176 [Abulk] <-km-- [Asurf] + [L] <-koff-- [AL]
177
178 Derivation::
179
180 d[AL]/dt = (kon*[A]*[L] - koff*[AL]) / (1 + kon*[L]/kmtl)
181 y = [AL]
182 (rmax-y) = [L]
183 conc = [Abulk]
184 rmax = [AL] + [L]
185 dy/dt = (kon*conc*(rmax-y) - koff*y) / (1 + kon*(rmax-y)/kmtl)
186
187 ======================= ============================================
188 Model Parameter Description
189 ======================= ============================================
190 params['t1']['value'] time of injection for binding, (s)
191 params['rmax']['value'] maximum response, (RIU)
192 params['conc']['value'] concentration of analyte [Abulk], (M)
193 params['kon']['value'] on-rate of analyte, (1/Ms)
194 params['t2']['value'] time of end binding & begin washing, (s)
195 params['koff']['value'] off-rate of analyte, (1/s)
196 params['t3']['value'] time end of washing & data fitting, (s)
197 params['kmtl']['value'] rate of diffusion, (RIU/Ms)
198 ----------------------- --------------------------------------------
199 *Optional*
200 params['cofa']['value'] concentration factor, (1/dilution factor)
201 ======================= ============================================
202
203 """
204
205 ## Skip parameter validation steps for now.
206 t1 = params['t1']['value']
207 rmax = params['rmax']['value']
208 conc = params['conc']['value']
209 kon = params['kon']['value']
210 t2 = params['t2']['value']
211 koff = params['koff']['value']
212 t3 = params['t3']['value']
213 kmtl = params['kmtl']['value']
214 if ('cofa' in params.keys()):
215 conc *= float(params['cofa']['value'])
216
217 ## Initialize variables.
218 stat = {} ## Error status dictionary.
219 y = np.zeros(len(time), dtype=float)
220 base = y[0] = data[0] ## Baseline SPR signal.
221
222 ## Error checks.
223 if (kmtl<=0):
224 ## Avoid div/0, assume user wants no mass transport limitation.
225 stat['kmtl must be > 0'] = True
226 kmtl = 1e40 ## An arbitrary very large number.
227
228 ## Must iterate through data, numerical integration.
229 for i in range(1,len(time)):
230 dt = time[i] - time[i-1]
231 ## Treat function as having three separate phases.
232 if (time[i] <= t1):
233 ## Pre-binding phase.
234 base = y[i] = data[i]
235 elif (t1 < time[i] <= t2):
236 ## Binding phase
237 yb = y[i-1] - base
238 dy = (conc*kon*(rmax-yb) - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
239 ## Check for integration errors producing wild results.
240 if (abs(y[i-1]) > 9e9): dy = 0; stat['y>9e9'] = True
241 y[i] = y[i-1] + dy*dt
242 elif (t2 < time[i] <= t3):
243 ## Dissociation (conc=0)
244 yb = y[i-1] - base
245 dy = (0 - koff*yb) / (1 + kon*(rmax-yb)/kmtl)
246 y[i] = y[i-1] + dy*dt
247
248 if (len(stat.keys()) > 0): print "Errors in simple1to1:", stat.keys()
249
250 return y
251 """End of simple1to1_mtl() function"""
252
253
254
255 ################################# End of module ###############################