-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathModel_PC12_48.py
174 lines (141 loc) · 4.88 KB
/
Model_PC12_48.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 19 19:28:35 2014
@author: Andrew
"""
import numpy as np
from modeling import BaseModel, DAECalc, Composition, ParallelModel
calc = DAECalc.DAECalc("__PC12_48__","__dPC12_48__","__d2PC12_48__")
class PC12_48(BaseModel):
def __init__(self, ts, data = 0, weights = 1.0):
self.ts = ts
self.weights = weights
self.data = data
self.calc = calc #DAECalc.DAECalc("__PC12_48__","__dPC12_48__","__d2PC12_48__")
BaseModel.__init__(self,len(ts) * 15, 69, "PC12_48") #15 dVars, 69 Parameters
self.calc.kwargs['max_steps'] = 25000
def r(self, x): #residual function
return (self.calc.evaluate(x, self.ts).flatten() - self.data)*self.weights #A matrix flattened into a 300 value vector
def v(self, x, v):
return self.calc.evaluate_derivative(x, v, self.ts).flatten()*self.weights
class Expt(BaseModel):
def __init__(self, x0):
self.x0 = x0 # x0 contains the default values for the experiment
BaseModel.__init__(self, 69, 48, "Expt")
def r(self, x): #residuals
return np.append(self.x0, x)
expts = np.loadtxt("Expts.txt")
class Prior(BaseModel):
def __init__(self, x0 = 1, weights = 25.0):
self.x0 = x0
self.weights = weights
BaseModel.__init__(self,48,48, "LinearPrior")
def r(self, x):
return np.log(x/(self.x0))*self.weights #(x - self.x0) * self.weights
def v(self, x, v):
return (v/x) * self.weights
# # expt = expts[0]
# expt = np.ones(21)
# expt[5] = 0.0 #125.025 # EGF
# expt[6] = 45.6
# expt[-1] = 0 # Zero
#########################
### Setup base model
##########################
ts = np.linspace(0, 120, 100)
basemodel = PC12_48(ts)
#########################
### First Apgar Expt
#########################
expt1 = np.ones(21)
expt1[-1] = 0 # Zero
expt1[5] = 125.025*1e-2 # EGF
expt1[6] = 45.6*1e2 # NGF
expt1[9] *= 100 # Sos
expt1[11] *= 100 # Ras
expt1[18] *= 100 # C3G
Expt1Model = Expt(expt1)
basemodel1 = PC12_48(ts)
model1 = Composition(basemodel1, Expt1Model)
#########################
### Second Apgar Expt
#########################
expt2 = np.ones(21)
expt2[-1] = 0 # Zero
expt2[5] = 125.025*1e-6 # EGF
expt2[6] = 45.6*1e-4 # NGF
expt2[14] *= 100 # Mek
expt2[15] *= 100 # Erk
expt2[4] *= 0 # Raf1PPtase
Expt2Model = Expt(expt2)
basemodel2 = PC12_48(ts)
model2 = Composition(basemodel2, Expt2Model)
#########################
### Third Apgar Expt
#########################
expt3 = np.ones(21)
expt3[-1] = 0 # Zero
expt3[5] = 0.0 # EGF
expt3[6] = 45.6*1e0 # NGF
expt3[13] *= 100 # BRaf
expt3[19] *= 100 # Rap1
expt3[2] *= 0 # RapGap
Expt3Model = Expt(expt3)
basemodel3 = PC12_48(ts)
model3 = Composition(basemodel3, Expt3Model)
#########################
### Fourth Apgar Expt
#########################
expt4 = np.ones(21)
expt4[-1] = 0 # Zero
expt4[5] = 125.025*1e-6 # EGF
expt4[6] = 45.6*1e2 # NGF
expt4[10] *= 100 # P90Rsk
expt4[16] *= 100 # PI3K
expt4[17] *= 100 # Akt
Expt4Model = Expt(expt4)
basemodel4 = PC12_48(ts)
model4 = Composition(basemodel4, Expt4Model)
#########################
### Fifth Apgar Expt
#########################
expt5 = np.ones(21)
expt5[-1] = 0 # Zero
expt5[5] = 125.025*1e-4 # EGF
expt5[6] = 45.6*1e-2 # NGF
expt5[12] *= 100 # Raf1
expt5[1] *= 0 # RasGap
Expt5Model = Expt(expt5)
basemodel5 = PC12_48(ts)
model5 = Composition(basemodel5, Expt5Model)
x = np.ones(48) #np.log(np.loadtxt("TrueParameters.txt"))
#expx = np.exp(x)
r = model1.r(x) #expx)
y1 = r.reshape(len(ts), 15)
sigma1 = np.maximum(0.1*y1.flatten(), 1e-2) ## 10 percent error bar
basemodel1.weights = 1.0/np.loadtxt("data/sigma1t.txt").flatten()
basemodel1.data = np.loadtxt("data/y1t.txt").flatten()
r = model2.r(x) #expx)
y2 = r.reshape(len(ts), 15)
sigma2 = np.maximum(0.1*y2.flatten(), 1e-2)
basemodel2.weights = 1.0/np.loadtxt("data/sigma2t.txt").flatten()
basemodel2.data = np.loadtxt("data/y2t.txt").flatten()
r = model3.r(x) #expx)
y3 = r.reshape(len(ts), 15)
sigma3 = np.maximum(0.1*y3.flatten(), 1e-2)
basemodel3.weights = 1.0/np.loadtxt("data/sigma3t.txt").flatten()
basemodel3.data = np.loadtxt("data/y3t.txt").flatten()
r = model4.r(x) #expx)
y4 = r.reshape(len(ts), 15)
sigma4 = np.maximum(0.1*y4.flatten(), 1e-2)
basemodel4.weights = 1.0/np.loadtxt("data/sigma4t.txt").flatten()
basemodel4.data = np.loadtxt("data/y4t.txt").flatten()
r = model5.r(x) #expx)
y5 = r.reshape(len(ts), 15)
sigma5 = np.maximum(0.1*y5.flatten(), 1e-2)
basemodel5.weights = 1.0/np.loadtxt("data/sigma5t.txt").flatten()
basemodel5.data = np.loadtxt("data/y5t.txt").flatten()
PriorModel = Prior(weights = 1.0)
from Exponential import Exponential
model = Composition(ParallelModel([model1, model2, model3, model4, model5]), Exponential(48))
#x = np.log(np.loadtxt("TrueParameters.txt"))