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 User Rev File contents
1 clausted 17 """
2 clausted 53 mdl_module
3     -------------------------------------------------------------------------------
4 clausted 17
5 clausted 53 **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 clausted 17 """
56 clausted 53 __version__ = "110215"
57 clausted 17
58    
59     ## Import libraries
60     import numpy as np
61 clausted 18 from copy import deepcopy
62 clausted 17
63    
64 clausted 18 def drift(time, data, params):
65 clausted 17 """
66     This function simply models a constant signal drift in units/second.
67 clausted 18 It requires numpy arrays of times and starting data values,
68 clausted 17 It only requires one parameter in the params list.
69 clausted 53
70 clausted 18 params['rate']['value']
71 clausted 53
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 clausted 17 """
80 clausted 18 y = np.zeros(len(time), dtype=float)
81 clausted 17 try:
82 clausted 18 rate = params['rate']['value']
83 clausted 17 except KeyError:
84     print "Error: Rate parameter not supplied to drift() model."
85 clausted 18 return y
86 clausted 17
87 clausted 18 ## 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 clausted 19 ## End of drift() function.
95 clausted 18
96    
97     def simple1to1(time, data, params):
98     """
99 clausted 26 This function simply models a 1:1 interaction.
100    
101 clausted 53 Model::
102 clausted 29
103 clausted 53 [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 clausted 18 """
131 clausted 53
132 clausted 18 ## Skip parameter validation steps for now.
133 clausted 29 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 clausted 17
143 clausted 18 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 clausted 26 dt = time[i] - time[i-1]
149     ## Treat function as having three separate phases.
150 clausted 18 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 clausted 26 y[i] = y[i-1] + dy*dt
159 clausted 18 elif (t2 < time[i] <= t3):
160 clausted 26 ## Dissociation (conc=0)
161 clausted 18 yb = y[i-1]-base
162     dy = 0 - koff*yb
163 clausted 26 y[i] = y[i-1] + dy*dt
164 clausted 18
165     return y
166 clausted 53 """End of simple1to1() function"""
167 clausted 18
168 clausted 19
169 clausted 26 def simple1to1_mtl(time, data, params):
170     """
171     This function simply models a 1:1 interaction with mass transport limitation.
172    
173 clausted 53 Model::
174 clausted 29
175 clausted 53 [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 clausted 26 """
204 clausted 53
205 clausted 26 ## 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 clausted 29 if ('cofa' in params.keys()):
215     conc *= float(params['cofa']['value'])
216    
217 clausted 53 ## Initialize variables.
218 clausted 26 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 clausted 53 """End of simple1to1_mtl() function"""
252 clausted 26
253    
254    
255 clausted 53 ################################# End of module ###############################