ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/models2010/modelparallel.py
Revision: 41
Committed: Tue Jan 18 00:35:23 2011 UTC (8 years, 8 months ago) by clausted
File size: 4753 byte(s)
Log Message:
Moved old data class "SPRdataclass" and accompanying surface interaction model modules to /models2010 subdirectory.  The plan is to implement these models for use with the "ba_class" and the modules in the parent directory.  

Should all the models be added to mdl_module or should they each go in their own module?  I am undecided.  
Line File contents
1 """
2 Provide the model for curvefitting.
3 This is a PARALLEL REACTION model.
4 The reaction equations:
5 L1 + A ==ka1=> L1A
6 L1 + A <=kd1== L1A
7 L2 + A ==ka2=> L2A
8 L2 + A <=kd2== L2A
9
10 Rui Hou, Yuhang Wan
11 Last modified on 100427 (yymmdd) by YW
12
13 Typical Pipeline:(work together with "modelclass.py" the father class):
14 >import modelparallel as mp
15 >m = mp.parallelmodel()
16 >...
17
18 """
19 __version__ = "100427"
20
21 import numpy as np
22 import pylab as plt
23 import copy
24 import os
25 import pickle
26 from modelclass import *
27
28
29 class parallelmodel(modelclass):
30 """The data attributes of the class:
31 parainfo: the parameter information
32 sim_data: the simulated data
33 """
34
35 def __init__ (self, parainfo=[], sim_data=[] ):
36 modelclass.__init__(self,parainfo,sim_data)
37
38 def wizard(self, ):
39 '''This function helps you to create a parameter information list.
40 This is a competing reactions model, where there are 10 parameters:
41 rmax1: Maximum analyte binding capacity for L1(RU),
42 ka1: Association rate constant for L1+A=L1A(M-1S-1),
43 kd1: Dissociation rate constant for L1A=L1+A(S-1),
44 rmax2: Maximum analyte binding capacity for L2(RU),
45 ka2: Association rate constant for L2+A=L2A(M-1S-1),
46 kd2: Dissociation rate constant for L2A=L2+A(S-1),
47 ca: Analyte concentration(M),
48 ton: Starting time for sample injection(s),
49 toff: Ending time for sample injection(s),
50 tout: Total time(s).'''
51 print ('')
52 pname = ['rmax1','ka1','kd1','rmax2','ka2','kd2','ca','ton','toff','tout']
53 parainfo = []
54 for i in range(10):
55 print i
56 tmp={}
57 tmp['name'] = pname[i]
58 vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
59 vtmp = vstr.strip().split(',')
60 tmp['number'] = len(vtmp)
61 if len(vtmp) == 1:
62 vtmp = float(vtmp[0])
63 else:
64 vtmp = map(float, vtmp)
65 tmp['value'] = vtmp
66 tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
67
68 parainfo.append(tmp)
69 ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
70 self.parainfo = parainfo
71
72
73 def function(self, t, paralist):
74 '''This function calculates the theoretical curve through the
75 parameter list you give.
76 '''
77 # for parallel reactions model
78
79 ## Assign the parameters to calculate the curve
80 for p in paralist:
81 if p['name'] == 'rmax1': rmax1 = p['value']
82 elif p['name'] == 'ka1': ka1 = p['value']
83 elif p['name'] == 'kd1': kd1 = p['value']
84 elif p['name'] == 'rmax2': rmax2 = p['value']
85 elif p['name'] == 'ka2': ka2 = p['value']
86 elif p['name'] == 'kd2': kd2 = p['value']
87 elif p['name'] == 'ca': ca = p['value']
88 elif p['name'] == 'ton': ton = p['value']
89 elif p['name'] == 'toff': toff = p['value']
90 elif p['name'] == 'tout': tout = p['value']
91 else: print p['name'], p['value']
92 if type(ca1) == list or type(ca2) == list:
93 print "Error: This function can only generate data for a single concentration."
94 return
95
96 t1 = ton ##100 ##234 # Initial injection
97 t2 = toff ##178 ##570 # Wash, begin of the dissociation
98 t3 = tout ##400 ##840 # End wash
99
100 ## Must iterate through data, numerical integration.
101 g = np.zeros(len(t), dtype=float)
102 g1 = np.zeros(len(t), dtype=float)
103 g2 = np.zeros(len(t), dtype=float)
104
105 for i in range(1,len(t)):
106 if (t[i] > t3): break #Speed things up.
107 if (t1 < t[i] < t2):
108 ## Association
109 dG1 = ka1*ca*(rmax1-g1[i-1]) - kd1*g1[i-1]
110 dG2 = ka2*ca*(rmax2-g2[i-1]) - kd2*g2[i-1]
111 elif (t2 < t[i] < t3):
112 ## Dissociation
113 dG1 = 0 - kd1*g1[i-1]
114 dG2 = 0 - kd2*g2[i-1]
115 else:
116 dG1 = 0
117 dG2 = 0
118 if (abs(g1[i]) > 999999999): dG1 = 0
119 if (abs(g2[i]) > 999999999): dG2 = 0
120
121 g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
122 g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
123 g[i] = g1[i] + g2[i]
124
125 return g
126
127
128 ##### End of parallel model class definition.