fix activated activity (#44)

* wip

* wip

* pkg

* pkg

* wip

* pkg

* fix PP1 errors

* fix activated activity

Manifest.toml

@@ -854,9 +854,9 @@ version = "1.0.2"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
-git-tree-sha1 = "6c22309e9a356ac1ebc5c8a217045f9bae6f8d9a"
+git-tree-sha1 = "627fcacdb7cb51dc67f557e1598cdffe4dda386d"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "1.10.13"
+version = "1.10.14"
 
 [[deps.HarfBuzz_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"]
@@ -937,9 +937,9 @@ weakdeps = ["Dates", "Test"]
     InverseFunctionsTestExt = "Test"
 
 [[deps.InvertedIndices]]
-git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
+git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
@@ -1210,9 +1210,9 @@ version = "2.38.0"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea"
+git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.28"
+version = "0.3.29"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -1316,9 +1316,9 @@ version = "1.11.0"
 
 [[deps.ModelingToolkit]]
 deps = ["AbstractTrees", "ArrayInterface", "BlockArrays", "Combinatorics", "CommonSolve", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "DiffRules", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "DynamicQuantities", "EnumX", "ExprTools", "Expronicon", "FindFirstFunctions", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "Graphs", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "NaNMath", "NonlinearSolve", "OffsetArrays", "OrderedCollections", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SCCNonlinearSolve", "SciMLBase", "SciMLStructures", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"]
-git-tree-sha1 = "4194c3e88a4dd91526150a8b0fb8456eb4056279"
+git-tree-sha1 = "d3fd8c745a20d1548652c3cdb1e2f0f03541cde0"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "9.56.0"
+version = "9.58.0"
 
     [deps.ModelingToolkit.extensions]
     MTKBifurcationKitExt = "BifurcationKit"
@@ -2094,9 +2094,9 @@ version = "1.2.1"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "d0553ce4031a081cc42387a9b9c8441b7d99f32d"
+git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.7"
+version = "1.4.8"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -2206,9 +2206,9 @@ version = "0.3.9"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
+git-tree-sha1 = "64cca0c26b4f31ba18f13f6c12af7c85f478cfde"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.4.0"
+version = "2.5.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]

docs/iso-sens.jl

@@ -1,7 +1,7 @@
 # # Sensitivity to ISO
 using ModelingToolkit
 using OrdinaryDiffEq
-using SteadyStateDiffEq
+using DiffEqCallbacks
 using Plots
 using LsqFit
 using CaMKIIModel
@@ -13,33 +13,42 @@ Plots.default(lw=1.5)
 sys = get_bar_sys(ATP, ISO; simplify=true)
 
 #---
-prob = SteadyStateProblem(sys, [])
-alg = DynamicSS(Rodas5P())
+tend = 5000second
+prob = ODEProblem(sys, [], tend)
+alg = Rodas5P()
 
 # Log scale for ISO concentration
-iso = logrange(1e-4μM, 1μM, length=101)
+iso = logrange(1e-4μM, 3μM, length=1001)
 
 #---
 prob_func = (prob, i, repeat) -> begin
     remake(prob, p=[ISO => iso[i]])
 end
 trajectories = length(iso)
-sol = solve(prob, alg, abstol=1e-10, reltol=1e-10) ## warmup
-sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories, abstol=1e-10, reltol=1e-10)
+callback = TerminateSteadyState(1e-10, 1e-10)
+sol = solve(prob, alg; callback) ## warmup
+sim = solve(EnsembleProblem(prob; prob_func, safetycopy=false), alg; trajectories, callback)
 
 #---
 """Extract values from ensemble simulations by a symbol"""
-extract(sim, k) = map(s -> s[k], sim)
+extract(sim, k) = map(s -> s[k][end], sim)
 """Calculate Root Mean Square Error (RMSE)"""
 rmse(fit) = sqrt(sum(abs2, fit.resid) / length(fit.resid))
 
 #---
 xopts = (xlims=(iso[begin], iso[end]), minorgrid=true, xscale=:log10, xlabel="ISO (μM)",)
 plot(iso, extract(sim, sys.cAMP); lab="cAMP", ylabel="Conc. (μM)", legend=:topleft, xopts...)
 
+#---
+plot(iso, [extract(sim, sys.PKACI) extract(sim, sys.PKACII)], legend=:topleft, lab=["PKACI" "PKACII"]; xopts...)
+
+#---
+plot(iso, [extract(sim, sys.I1) extract(sim, sys.I1p) extract(sim, sys.I1p_PP1)], legend=:topleft, lab=["I1" "I1p" "I1p_PP1"]; xopts...)
+
 #---
 plot(iso, extract(sim, sys.PKACI / sys.RItot); lab="PKACI", ylabel="Activation fraction")
-plot!(iso, extract(sim, sys.PKACII / sys.RIItot), lab="PKACII", legend=:topleft; xopts...)
+plot!(iso, extract(sim, sys.PKACII / sys.RIItot), lab="PKACII")
+plot!(iso, extract(sim, sys.PP1 / sys.PP1totBA), lab="PP1", legend=:topleft; xopts...)
 
 # ## Least-square fitting of PKACI activity
 @. model(x, p) = p[1] * x / (x + p[2]) + p[3]
@@ -53,7 +62,7 @@ pkac1_coef = coef(pkac1_fit)
 #---
 println("PKACI")
 println("Basal activity: ", pkac1_coef[3])
-println("Maximal activity: ", pkac1_coef[1] + pkac1_coef[3])
+println("Activated activity: ", pkac1_coef[1])
 println("Michaelis constant: ", pkac1_coef[2], " μM")
 println("RMSE: ", rmse(pkac1_fit))
 
@@ -74,7 +83,7 @@ pkac2_coef = coef(pkac2_fit)
 #---
 println("PKACII")
 println("Basal activity: ", pkac2_coef[3])
-println("Maximal activity: ", pkac2_coef[1] + pkac2_coef[3])
+println("Activated activity: ", pkac2_coef[1])
 println("Michaelis constant: ", pkac2_coef[2], " μM")
 println("RMSE: ", rmse(pkac2_fit))
 
@@ -95,7 +104,7 @@ pp1_coef = coef(pp1_fit)
 
 #---
 println("PP1")
-println("Basal activity: ", pp1_coef[1] + pp1_coef[3])
+println("Repressible activity: ", pp1_coef[1])
 println("Minimal activity: ", pp1_coef[3])
 println("Repressive Michaelis constant: ", pp1_coef[2], " μM")
 println("RMSE: ", rmse(pp1_fit))
@@ -109,7 +118,7 @@ plot(p1, p2, layout=(2, 1), size=(600, 600))
 # ## Least-square fitting of PLBp
 xdata = iso
 ydata = extract(sim, sys.PLBp / sys.PLBtotBA)
-plot(xdata, ydata, title="PLBp activity", lab=false; xopts...)
+plot(xdata, ydata, title="PLBp fraction", lab=false; xopts...)
 
 # First try: Hill function
 @. model_plb(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
@@ -132,16 +141,16 @@ function plbp_analytic(iso)
     PLBtot = 106μM
     PKACItot = 1.18μM
     PP1tot = 0.89μM
-    k_PKA_PLB = 54Hz
+    k_PKA_PLB = 54Hz / μM
     Km_PKA_PLB = 21μM
-    k_PP1_PLB = 8.5Hz
+    k_PP1_PLB = 8.5Hz / μM
     Km_PP1_PLB = 7.0μM
-    PKACI_basal = 0.0831  ## basal activity
-    PKACI_activated = 0.25603
-    PKACI_KM = 0.0144μM
-    PP1_basal = 0.82365
-    PP1_activated = 0.1025
-    PP1_KI = 0.008465μM
+    PKACI_basal = 0.0734  ## basal activity
+    PKACI_activated = 0.1995
+    PKACI_KM = 0.0139μM
+    PP1_basal = 0.8927
+    PP1_activated = 0.0492
+    PP1_KI = 0.00637μM
 
     ## Solve for Vf * x / (x + k1) = Vr * (1 - x) / (1 - x + k2)
     PKACI = PKACItot * (PKACI_basal + PKACI_activated * hil(iso, PKACI_KM))