-
Notifications
You must be signed in to change notification settings - Fork 3
/
runKai.hs
223 lines (180 loc) · 7.74 KB
/
runKai.hs
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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
import CpiLib
import CpiParser
import CpiSemantics
import CpiODE(timePoints,speciesIn,timeSeries)
import CpiLogic
import CpiPlot
import CpiTest(tEnv,tProc)
import CpiMatlab
import qualified TestKai as K
import Text.ParserCombinators.Parsec
import System.IO
import Data.List as L
import qualified Data.Map as Map
import qualified Numeric.GSL as GSL
import qualified Numeric.LinearAlgebra as LA
import qualified Graphics.Plot as Plot
import Control.Concurrent(forkIO)
import System.TimeIt
-- time points:
tps = (720,(0,72))
ts = timePoints (fst tps) (snd tps)
----------------
-- species:
----------------
findSpec str = maybe Nil id (lookupSpecName K.env str)
cc0 = findSpec "CC0"
cc1 = findSpec "CC1"
cc2 = findSpec "CC2"
cc3 = findSpec "CC3"
cc4 = findSpec "CC4"
cc5 = findSpec "CC5"
cc6 = findSpec "CC6"
c0 = findSpec "C0"
c1 = findSpec "C1"
c2 = findSpec "C2"
c3 = findSpec "C3"
c4 = findSpec "C4"
c5 = findSpec "C5"
c6 = findSpec "C6"
kaib = findSpec "B"
kaia = findSpec "A"
----------------------
-- formulae:
----------------------
-- F(([C6]+[CC6])>0)
-- Eventually we get some fully phosphorylated KaiC hexamers:
fullPhos = Pos (0,infty) (ValGT (Plus (Conc c6) (Conc cc6)) (R 0.0))
-- Phosphorylation level:
-- (([C0]+[CC0])*0)+...+(([C6]+[CC6])*6)
-- -------------------------------------
-- ([C0]+[CC0])+...+([C6]+[CC6])
phosL = Quot
(Plus (Times (Plus (Conc c0) (Conc cc0)) (R 0))
(Plus (Times (Plus (Conc c1) (Conc cc1)) (R 1))
(Plus (Times (Plus (Conc c2) (Conc cc2)) (R 2))
(Plus (Times (Plus (Conc c3) (Conc cc3)) (R 3))
(Plus (Times (Plus (Conc c4) (Conc cc4)) (R 4))
(Plus (Times (Plus (Conc c5) (Conc cc5)) (R 5))
(Times (Plus (Conc c6) (Conc cc6)) (R 6))))))))
(Plus (Plus (Conc c0) (Conc cc0))
(Plus (Plus (Conc c1) (Conc cc1))
(Plus (Plus (Conc c2) (Conc cc2))
(Plus (Plus (Conc c3) (Conc cc3))
(Plus (Plus (Conc c4) (Conc cc4))
(Plus (Plus (Conc c5) (Conc cc5))
(Plus (Conc c6) (Conc cc6))))))))
-- F(phosL >= x)
-- Phosphorylation level reaches x (in {0..6}).
reaches x = Pos (0,infty) (ValGE phosL (R x))
-- F(G(phosL >= x))
-- Phos. level reaches x and remains >=x indefinitely
remains x = Pos (0,infty) (Nec (0,infty) (ValGE phosL (R x)))
{-- Oscillation around phos level x?
-- G(((phosL==x)==>F(phosL/=x))AND((phosL/=x)==>F(phosL==k)))
osc x = Nec (Conj (Impl (ValEq phosL (R x)) (Pos (ValNEq phosL (R x))))
(Impl (ValNEq phosL (R x)) (Pos (ValEq phosL (R x)))))
-}
-- Oscillation around phos level x
-- F(¬(G(phosL < x)OR(G(phosL > y))))
oscB x = Pos (0,infty) (Neg (Disj (Nec (0,infty) (ValLT phosL (R x))) (Nec (0,infty) (ValGT phosL (R x)))))
-- Oscillation of a species concentration around conc x:
-- (formula as for oscB)
oscS spec x = Pos (0,infty) (Neg (Disj (Nec (0,infty) (ValLT (Conc spec) (R x))) (Nec (0,infty) (ValGT (Conc spec) (R x)))))
-- [0.56]A |> oscB(3)
-- Introducing some A gives oscillation?
gtee1 = Gtee pA (oscB 3)
where
pA = Process [(kaia,"0.56")] (AffNet [])
-- ¬([0.5]Inhib(in):{a-in@2e5} |> oscB(3))
-- Introducing Inhib kills osc.
-- Inib binds KaiB
ginhib = Neg (Gtee pIn (oscB 3))
pIn = Process [(Def "Inhib" ["in"],"0.5")] (AffNet [Aff (("in","a"),"2e5")])
-- Inhib(in) = {x-u@5.0} in<x>.u.Inhib(in)
inhib = SpeciesDef "Inhib" ["in"]
(New (AffNet [Aff (("x","u"),"5.0")])
(Sum [(Comm "in" ["x"] [],Sum [(Comm "u" [] [],Def "Inhib" ["in"])])]))
-- G(¬([0.5]Inhib(in):{a-in@2e5} |> oscB(3)))
-- nested gtee: Whenever we introduce Inhib it kills the osc.?
ginhibG = Nec (0,infty) (Neg (Gtee pIn (oscB 3)))
-- New time-bounded Osc formula: (Oscillates between x and y (amplitude >= y-x) for at least t with period <= p.
oscT x y t p = Nec (0,t) (Conj (Pos (0,(t+p)) (ValLT phosL (R x))) (Pos (0,(t+p)) (ValGT phosL (R y))))
-- Introducing x Inhib kills the oscillation (using time bounded formula)?
gInhibOscT x = Gtee (pIn' x) (oscT 2 4 48 24)
pIn' x = Process [(Def "Inhib" ["in"],(d2s x))] (AffNet [Aff (("in","a"),"2e5")])
--------------
-- computation
--------------
main = do
-- get time series:
putStrLn "Solving KaiABC model..."
let soln = solveODEoctave K.env K.kai K.dpdt tps
timeIt $ putStrLn ("...done: " ++ show (LA.cols soln) ++ " species.")
let ss = speciesIn K.env K.dpdt
ss' = speciesInProc K.kai
trace = timeSeries ts soln ss
{- -- plot time series
putStrLn "Plotting..."
plotTimeSeriesToFileFiltered ts soln ss ss' "KaiABC-plot.pdf"
-}
-- We oscillate to begin with?
putStrLn ("Checking: OscT(2,4,48,24): ")
let oscBasic = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (oscT 2 4 48 24)
timeIt $ putStrLn ("...done: " ++ show oscBasic)
-- Introducing varying concs of Inhib kills the osc?
putStrLn ("Checking: [0.2]Inhib |> OscT(2,4,48,24): ")
let oscInhib2 = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (gInhibOscT 0.2)
timeIt $ putStrLn ("...done: " ++ show oscInhib2)
putStrLn ("Checking: [0.3]Inhib |> OscT(2,4,48,24): ")
let oscInhib2 = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (gInhibOscT 0.3)
timeIt $ putStrLn ("...done: " ++ show oscInhib2)
putStrLn ("Checking: [0.4]Inhib |> OscT(2,4,48,24): ")
let oscInhib2 = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (gInhibOscT 0.4)
timeIt $ putStrLn ("...done: " ++ show oscInhib2)
{-
-- We get some fully phos KaiC?
putStrLn ("Checking :" ++ pretty fullPhos)
let fp = modelCheck K.env solveODEoctave (Just trace) K.kai tps (fullPhos)
timeIt $ putStrLn ("...done: " ++ show fp)
-- Phosphorylation level:
putStrLn ("Checking: F(phosL>=6)")
let r6 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (reaches 6)
timeIt $ putStrLn ("...done: " ++ show r6)
putStrLn ("Checking: F(phosL>=5)")
let r5 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (reaches 5)
timeIt $ putStrLn ("...done: " ++ show r5)
-- Phos level reaches and remains above x
putStrLn ("Checking: F(G(PhosL>=5))")
let rem5 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (remains 5)
timeIt $ putStrLn ("...done: " ++ show rem5)
{-- Phos level oscillates around 3:
putStrLn ("Checking: Osc(3): ") -- ++ pretty (osc 3)) {- not pretty! -}
let osc3 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (osc 3)
timeIt $ putStrLn ("...done: " ++ show osc3)-}
-- Phos level oscillates between around 3:
putStrLn ("Checking: OscB(3): ")
let oscB3 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (oscB 3)
timeIt $ putStrLn ("...done: " ++ show oscB3)
-- Conc of KaiB oscillates around 1M
putStrLn ("Checking: OscS([B],1): ")
let oscSB = modelCheck K.env solveODEoctave (Just trace) K.kai tps (oscS kaib 1)
timeIt $ putStrLn ("...done: " ++ show oscSB)
-- Conc of KaiB oscillates around 1M
putStrLn ("Checking: OscS([A],0.1): ")
let oscSA = modelCheck K.env solveODEoctave (Just trace) K.kai tps (oscS kaia 0.1)
timeIt $ putStrLn ("...done: " ++ show oscSA)
{-- Introducing some KaiA causes oscillation?
putStrLn ("Checking: [0.56]A |> OscB(3): ")
let gt1 = modelCheck K.env solveODEoctave (Just trace) K.kai tps (gtee1)
timeIt $ putStrLn ("...done: " ++ show gt1)-}
-- Introducing some Inhib kills the osc.?
putStrLn ("Checking: ¬([0.5]Inhib |> OscB(3)): ")
let gin = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (ginhib)
timeIt $ putStrLn ("...done: " ++ show gin)
{-- Introducing some Inhib at any time kills the osc.?
-- Takes way too long to check right now....
putStrLn ("Checking: G(¬([0.5]Inhib |> OscB(3))): ")
let ginG = modelCheck (inhib : K.env) solveODEoctave (Just trace) K.kai tps (ginhibG)
timeIt $ putStrLn ("...done: " ++ show ginG)-}
-}