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, 8 months ago) by clausted
File size: 9274 byte(s)
Log Message:
Added fit_module and mdl_module to sphinx documentation.
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 17
139 clausted 18 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 clausted 26 dt = time[i] - time[i-1]
145     ## Treat function as having three separate phases.
146 clausted 18 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 clausted 26 y[i] = y[i-1] + dy*dt
155 clausted 18 elif (t2 < time[i] <= t3):
156 clausted 26 ## Dissociation (conc=0)
157 clausted 18 yb = y[i-1]-base
158     dy = 0 - koff*yb
159 clausted 26 y[i] = y[i-1] + dy*dt
160 clausted 18
161     return y
162 clausted 53 """End of simple1to1() function"""
163 clausted 18
164 clausted 19
165 clausted 26 def simple1to1_mtl(time, data, params):
166     """
167     This function simply models a 1:1 interaction with mass transport limitation.
168 clausted 54 The model parameters are described in the table below.
169 clausted 26
170 clausted 53 Model::
171 clausted 29
172 clausted 53 [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 clausted 26 """
201 clausted 53
202 clausted 26 ## 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 clausted 29 if ('cofa' in params.keys()):
212     conc *= float(params['cofa']['value'])
213    
214 clausted 53 ## Initialize variables.
215 clausted 26 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 clausted 53 """End of simple1to1_mtl() function"""
249 clausted 26
250    
251    
252 clausted 53 ################################# End of module ###############################