ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/models2010/modelcompeting_varyC.py
Revision: 41
Committed: Tue Jan 18 00:35:23 2011 UTC (8 years, 9 months ago) by clausted
File size: 5205 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 COMPETING REACTION model with time variable concentration.
4 The reaction equations are:
5 L + A1 ==ka1=> LA1
6 L + A1 <=kd1== LA1
7 L + A2 ==ka2=> LA2
8 L + A2 <=kd2== LA2
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 modelcompeting_varyC as mcvc
15 >m = mcvc.competingmodel_varyC()
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 class competingmodel_varyC(modelclass):
29 """The data attributes of the class:
30 parainfo: the parameter information
31 sim_data: the simulated data
32 conc: the time variable concentration data
33 """
34
35 def __init__ (self, parainfo=[], sim_data=[], conc=[] ):
36 modelclass.__init__(self,parainfo,sim_data)
37 self.conc = conc
38
39 ##------------ not finished yet---------------
40 def load_conc(self, fname):
41 ##fp = file(fname, 'r+')
42 conc = np.loadtxt(fname)
43 self.conc = conc
44 return
45
46 def create_conc(self, t):
47 # create a step-function as concentration curve for debugging now
48 c = np.zeros(len(t), dtype=float)
49 conc = np.zeros((2,len(t)), dtype=float)
50 t = np.array(t)
51 t1 = input("start time:")
52 t2 = input("end time:")
53 ind = plt.find((t>t1-(t[1]-t[0]))&(t<t2-(t[1]-t[0])))
54 c[ind] = 1
55 conc[0] = t
56 conc[1] = c
57 self.conc = conc
58 plt.plot(t,c,'.')
59 plt.show()
60 return
61
62 ##----------------for debugging now------------
63
64
65 def wizard(self, ):
66 '''This function helps you to create a parameter information list.
67 This is a competing reactions model, where there are 9 parameters:
68 rmax: Maximum analyte binding capacity(RU),
69 mw1: (Relative) molecular weight of A1,
70 ka1: Association rate constant for L+A1=LA1(M-1S-1),
71 kd1: Dissociation rate constant for LA1=L+A1(S-1),
72 mw2: (Relative) molecular weight of A2,
73 ka2: Association rate constant for L+A2=LA2(M-1S-1),
74 kd2: Dissociation rate constant for LA2=L+A2(S-1),
75 ca1: Analyte concentration for A1(M),
76 ca2: Analyte concentration for A2(M).'''
77 print ('')
78 pname = ['rmax','mw1','ka1','kd1','mw2','ka2','kd2','ca1','ca2']
79 parainfo = []
80 for i in range(9):
81 print i
82 tmp={}
83 tmp['name'] = pname[i]
84 vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
85 vtmp = vstr.strip().split(',')
86 tmp['number'] = len(vtmp)
87 if len(vtmp) == 1:
88 vtmp = float(vtmp[0])
89 else:
90 vtmp = map(float, vtmp)
91 tmp['value'] = vtmp
92 tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
93
94 parainfo.append(tmp)
95 ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
96 self.parainfo = parainfo
97
98
99 def function(self, t, paralist):
100 '''This function calculates the theoretical curve through the
101 parameter list you give.
102 '''
103 # for competing reactions model
104 ## Assign the parameters to calculate the curve
105 conc = copy.deepcopy(self.conc)
106 C = np.interp(t,conc[0],conc[1])
107 for p in paralist:
108 if p['name'] == 'rmax': rmax = p['value']
109 elif p['name'] == 'mw1': mw1 = p['value']
110 elif p['name'] == 'ka1': ka1 = p['value']
111 elif p['name'] == 'kd1': kd1 = p['value']
112 elif p['name'] == 'mw2': mw2 = p['value']
113 elif p['name'] == 'ka2': ka2 = p['value']
114 elif p['name'] == 'kd2': kd2 = p['value']
115 elif p['name'] == 'ca1': ca1 = p['value']
116 elif p['name'] == 'ca2': ca2 = p['value']
117 else: print p['name'], p['value']
118 if type(ca1) == list or type(ca2) == list:
119 print "Error: This function can only generate data for a single concentration."
120 return
121
122 rau = float(mw1/mw2)
123
124 ## Must iterate through data, numerical integration.
125 g = np.zeros(len(C), dtype=float)
126 g1 = np.zeros(len(C), dtype=float)
127 g2 = np.zeros(len(C), dtype=float)
128
129 for i in range(1,len(C)):
130 dG1 = ka1*ca1*C[i-1]*(rmax-g1[i-1]-g2[i-1]*rau) - kd1*g1[i-1]
131 dG2 = ka2*ca2*C[i-1]*(rmax/rau-g1[i-1]/rau-g2[i-1]) - kd2*g2[i-1]
132
133 if (abs(g1[i]) > 999999999): dG1 = 0
134 if (abs(g2[i]) > 999999999): dG2 = 0
135
136 g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
137 g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
138 g[i] = g1[i] + g2[i]
139
140 return g
141
142
143 ##### End of time variable concentrated competing model class definition.