ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/modeltwostate.py
Revision: 25
Committed: Wed Apr 28 20:22:24 2010 UTC (9 years, 2 months ago) by rjaynes
File size: 4634 byte(s)
Log Message:
Add py and obj files to allow modeling of more SPR experiments with converter and curvefitting modules.  This is the work of Yuhang Wan and Rui Hou.

1. In "converter.py": 
      Add the saving and reading function for the sprclass data object.
      Also add function "keyfile_read_fake" to provide default information for SPRit and ICM formats in case of the bug when do background_subtract.
      Fix the bugs in "background_subtract".
      Tested by DAM and ICM formats.
2. In model modules:
      "modelclass.py" is the parent class for all the other model classes that performs the theoretical simulating, loading and saving of the parameter or simulated data. Rui and I also add some other model modules like competing model, twostate model, parallel model, and the time variable concentrated models, where the simulated result is compared with Clamp's simulation to make sure the equations are correct. 
       The basicmodel and basicmodel_varyC class are tested. 
3. In "curvefitting.py":
      Add typical pipeline for operation. The examples are packed with the file. 
      Add function to show the Elapsed time for each fitting.
Line File contents
1 """
2 Provide the model for curvefitting.
3 This is a TWO STATE REACTION model.
4 The reaction equations are:
5 L + A ==ka1=> LA
6 L + A <=kd1== LA
7 LA ==ka2=> LA*
8 LA <=kd2== LA*
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 modeltwostate as mt
15 >m = mt.twostatemodel()
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 twostatemodel(modelclass):
30 """The data attributes of the class:
31 parainfo: the parameter information
32 sim_data: the simulated data
33 """
34 def __init__ (self, parainfo=[], sim_data=[] ):
35 modelclass.__init__(self,parainfo,sim_data)
36
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 9 parameters:
41 rmax: Maximum analyte binding capacity(RU),
42 ka1: Association rate constant for L+A=LA(M-1S-1),
43 kd1: Dissociation rate constant for LA=L+A(S-1),
44 ka2: Forward rate constant for LA=LA*(S-1),
45 kd2: Backward rate constant for LA*=LA(S-1),
46 ca: Analyte concentration(M),
47 ton: Starting time for sample injection(s),
48 toff: Ending time for sample injection(s),
49 tout: Total time(s).'''
50 print ('')
51 pname = ['rmax','ka1','kd1','ka2','kd2','ca','ton','toff','tout']
52 parainfo = []
53 for i in range(9):
54 print i
55 tmp={}
56 tmp['name'] = pname[i]
57 vstr = raw_input("input the value of parameter '%s', if a series, seperate with ',': " %pname[i])
58 vtmp = vstr.strip().split(',')
59 tmp['number'] = len(vtmp)
60 if len(vtmp) == 1:
61 vtmp = float(vtmp[0])
62 else:
63 vtmp = map(float, vtmp)
64 tmp['value'] = vtmp
65 tmp['fixed'] = bool(input("Is '%s' fixed? (1/0): " %pname[i]))
66
67 parainfo.append(tmp)
68 ##parainfo.append({'name':'', 'value':0., 'fixed':0, 'number':0})
69 self.parainfo = parainfo
70
71
72 def function(self, t, paralist):
73 '''This function calculates the theoretical curve through the
74 parameter list you give.
75 '''
76 # for two state reaction model
77
78 ## Assign the parameters to calculate the curve
79 for p in paralist:
80 if p['name'] == 'rmax': rmax = p['value']
81 elif p['name'] == 'ka1': ka1 = p['value']
82 elif p['name'] == 'kd1': kd1 = p['value']
83 elif p['name'] == 'ka2': ka2 = p['value']
84 elif p['name'] == 'kd2': kd2 = p['value']
85 elif p['name'] == 'ca': ca = p['value']
86 elif p['name'] == 'ton': ton = p['value']
87 elif p['name'] == 'toff': toff = p['value']
88 elif p['name'] == 'tout': tout = p['value']
89 else: print p['name'], p['value']
90 if type(ca1) == list or type(ca2) == list:
91 print "Error: This function can only generate data for a single concentration."
92 return
93
94 t1 = ton ##100 ##234 # Initial injection
95 t2 = toff ##178 ##570 # Wash, begin of the dissociation
96 t3 = tout ##400 ##840 # End wash
97
98 ## Must iterate through data, numerical integration.
99 g = np.zeros(len(t), dtype=float)
100 g1 = np.zeros(len(t), dtype=float)
101 g2 = np.zeros(len(t), dtype=float)
102
103 for i in range(1,len(t)):
104 if (t[i] > t3): break #Speed things up.
105 if (t1 < t[i] < t2):
106 ## Association
107 dG1 = (ka1*ca*(rmax-g1[i-1]-g2[i-1]) - kd1*g1[i-1])-(ka2*g1[i-1]- kd2*g2[i-1])
108 dG2 = ka2*g1[i-1]- kd2*g2[i-1]
109 elif (t2 < t[i] < t3):
110 ## Dissociation
111 dG1 = 0 - kd1*g1[i-1] - ka2*g1[i-1] + kd2*g2[i-1]
112 dG2 = ka2*g1[i-1]- kd2*g2[i-1]
113 else:
114 dG1 = 0
115 dG2 = 0
116 if (abs(g1[i]) > 999999999): dG1 = 0
117 if (abs(g2[i]) > 999999999): dG2 = 0
118
119 g1[i] = g1[i-1] + dG1 * (t[i] - t[i-1])
120 g2[i] = g2[i-1] + dG2 * (t[i] - t[i-1])
121 g[i] = g1[i] + g2[i]
122
123 return g
124
125 ##### End of twostate model class definition.