Skip to content

Commit

Permalink
simplify ISO sys more
Browse files Browse the repository at this point in the history
  • Loading branch information
@sosiristseng
sosiristseng committed Dec 13, 2024
1 parent 7373d91 commit 5a33e3f
Show file tree
Hide file tree
Showing 2 changed files with 208 additions and 108 deletions.
166 changes: 163 additions & 3 deletions docs/iso-sens.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ plot(iso, extract(sim, sys.PKACI / sys.RItot); lab="PKACI", ylabel="Activation f
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
# ## Fitting active PKACI
@. model(x, p) = p[1] * x / (x + p[2]) + p[3]
xdata = iso
ydata = extract(sim, sys.PKACI / sys.RItot)
Expand All @@ -72,7 +72,7 @@ p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot],
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="PKACI error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Least-square fitting of PKACII activity
# ## Fitting active PKACII
xdata = iso
ydata = extract(sim, sys.PKACII / sys.RIItot)
p0 = [0.4, 0.01μM, 0.2]
Expand Down Expand Up @@ -115,7 +115,7 @@ p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot],
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="PP1 error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Least-square fitting of PLBp
# ## Fitting PLBp
xdata = iso
ydata = extract(sim, sys.PLBp / sys.PLBtotBA)
plot(xdata, ydata, title="PLBp fraction", lab=false; xopts...)
Expand All @@ -128,6 +128,11 @@ fit = curve_fit(model_plb, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("PLBp")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
Expand Down Expand Up @@ -176,3 +181,158 @@ plot(p1, p2, layout=(2, 1), size=(600, 600))

#---
println("RMSE: ", rmse(fit))

# ## Fitting PLMp
xdata = iso
ydata = extract(sim, sys.PLMp / sys.PLMtotBA)
plot(xdata, ydata, title="PLMp fraction", lab=false; xopts...)

#---
@. model_plm(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.5, 1e-9μM, 1.0, 0.0]
fit = curve_fit(model_plm, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("PLMp")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
ypred = model_plm.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="PLMp", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="PLMp error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Fitting TnIp
xdata = iso
ydata = extract(sim, sys.TnIp / sys.TnItotBA)
plot(xdata, ydata, title="TnIp fraction", lab=false; xopts...)

#---
@. model_tni(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.1, 1e-9μM, 1.0, 0.0]
fit = curve_fit(model_tni, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("TnIp")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
ypred = model_tni.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="TnIp", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="TnIp error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Fitting LCCap
xdata = iso
ydata = extract(sim, sys.LCCap / sys.LCCtotBA)
plot(xdata, ydata, title="LCCap fraction", lab=false; xopts...)

#---
@. model_lcc(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.1, 1e-9μM, 0.5, 0.0]
fit = curve_fit(model_lcc, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("LCCap")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
ypred = model_lcc.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="LCCap", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="LCCap error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Fitting LCCbp
xdata = iso
ydata = extract(sim, sys.LCCbp / sys.LCCtotBA)
plot(xdata, ydata, title="LCCbp fraction", lab=false; xopts...)

#---
@. model_lcc(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.1, 1e-9μM, 0.5, 0.0]
fit = curve_fit(model_lcc, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("LCCbp")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
ypred = model_lcc.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="LCCbp", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="LCCbp error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Fitting KURp
xdata = iso
ydata = extract(sim, sys.KURp / sys.IKurtotBA)
plot(xdata, ydata, title="LCCbp fraction", lab=false; xopts...)

#---
@. model_kur(x, p) = p[1] * hil(x, p[2], p[3]) + p[4]
p0 = [0.8, 1e-2μM, 1.0, 0.1]
lb = [0.1, 1e-9μM, 0.5, 0.0]
fit = curve_fit(model_kur, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("KURp")
println("Basal activity: ", pestim[4])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("Hill coefficient: ", pestim[3])
println("RMSE: ", rmse(fit))

#---
ypred = model_lcc.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="KURp", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="KURp error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))

# ## Fitting RyRp
xdata = iso
ydata = extract(sim, sys.RyR_PKAp)
plot(xdata, ydata, title="RyRp fraction", lab=false; xopts...)

#---
@. model(x, p) = p[1] * x / (x + p[2]) + p[3]
p0 = [0.3, 1e-2μM, 0.1]
lb = [0.0, 1e-9μM, 0.0]
fit = curve_fit(model, xdata, ydata, p0; lower=lb, autodiff=:forwarddiff)
pestim = coef(fit)

#---
println("RyRp")
println("Basal activity: ", pestim[3])
println("Activated activity: ", pestim[1])
println("Michaelis constant: ", pestim[2], " μM")
println("RMSE: ", rmse(fit))

#---
ypred = model.(xdata, Ref(pestim))
p1 = plot(xdata, [ydata ypred], lab=["Full model" "Fitted"], line=[:dash :dot], title="RyRp", legend=:topleft; xopts...)
p2 = plot(xdata, (ypred .- ydata) ./ ydata .* 100; title="RyRp error (%)", lab=false, xopts...)
plot(p1, p2, layout=(2, 1), size=(600, 600))
150 changes: 45 additions & 105 deletions src/isoproterenol.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,8 +211,9 @@ function get_bar_sys(ATP=5000μM, ISO=0μM; name=:bar_sys, simplify=false)
add_rate!(rates, kf_AC_Gsa, [AC, GsaGTP], kr_AC_Gsa, [AC_GsaGTP]) # AC + GsaGTP <--> AC_GsaGTP
add_rate!(rates, k_G_hyd, [AC_GsaGTP], 0, [AC, GsaGDP]) # AC_GsaGTP --> AC + GsaGDP
add_rate!(rates, k_PKA_PDE * PKACII, [PDE], k_PP_PDE, [PDEp]) # PDE <--> PDEp
v = k_AC_Gsa * AC_GsaGTP * hil(ATP, Km_AC_Gsa) + k_AC_basal * AC * hil(ATP, Km_AC_basal) - (k_cAMP_PDE * PDE + k_cAMP_PDEp * PDEp) * hil(cAMP, Km_PDE_cAMP)
add_raw_rate!(rates, v, [], [cAMP]) # 0 <--> cAMP
vf = k_AC_Gsa * AC_GsaGTP * hil(ATP, Km_AC_Gsa) + k_AC_basal * AC * hil(ATP, Km_AC_basal)
vr = (k_cAMP_PDE * PDE + k_cAMP_PDEp * PDEp) * hil(cAMP, Km_PDE_cAMP)
add_raw_rate!(rates, vf - vr, [], [cAMP]) # 0 <--> cAMP
# cAMP disinhibition of PKA
add_rate!(rates, kf_RC_cAMP, [RC_I, cAMP], kr_RC_cAMP, [RCcAMP_I]) # RC_I + cAMP <--> RCcAMP_I
add_rate!(rates, kf_RC_cAMP, [RC_II, cAMP], kr_RC_cAMP, [RCcAMP_II]) # RC_II + cAMP <--> RCcAMP_II
Expand Down Expand Up @@ -258,58 +259,9 @@ function get_bar_sys(ATP=5000μM, ISO=0μM; name=:bar_sys, simplify=false)
return sys
end

"Algebraic beta-adrenergic system."
"Algebraic fitted beta-adrenergic system"
function get_bar_sys_reduced(ISO=0μM; name=:bar_sys)
@parameters begin
PP1totBA = 890nM
PKACItot = 1.18μM
PKACIItot = 118nM
PKACII_LCCtotBA = 25nM
PKAIItot = 59nM
PKAII_KURtot = 25nM
PP1_KURtot = 25nM
LCCtotBA = 25nM
PLBtotBA = 106μM
PLMtotBA = 48μM
TnItotBA = 70μM
IKurtotBA = 25nM
RyRtotBA = 135nM # Total RyR content for the BAR module
PKAII_RyRtot = 34nM
epsilon = 10
k_PKA_PLB = 54Hz / μM
Km_PKA_PLB = 21μM
k_PP1_PLB = 8.5Hz / μM
Km_PP1_PLB = 7.0μM
k_PKA_PLM = 54Hz / μM
Km_PKA_PLM = 21μM
k_PP1_PLM = 8.5Hz / μM
Km_PP1_PLM = 7.0μM
PP1_LCC = 0.025μM
PP2A_LCC = 0.025μM
k_PKA_LCC = 54Hz / μM
Km_PKA_LCC = 21μM
k_PP1_LCC = 8.52Hz / μM
Km_PP1_LCC = 3μM
k_PP2A_LCC = 10.1Hz / μM
Km_PP2A_LCC = 3μM
PP2A_TnI = 0.67μM
k_PKA_TnI = 54Hz / μM
Km_PKA_TnI = 21μM
k_PP2A_TnI = 10.1Hz / μM
Km_PP2A_TnI = 4.1μM
k_pka_KUR = 54Hz / μM
Km_pka_KUR = 21μM
k_pp1_KUR = 8.52Hz / μM
Km_pp1_KUR = 7μM
PP1_RyR = 34nM
PP2A_RyR = 34nM
kcat_pka_RyR = 54Hz / μM
Km_pka_RyR = 21μM
kcat_pp1_RyR = 8.52Hz / μM
Km_pp1_RyR = 7μM
kcat_pp2a_RyR = 10.1Hz / μM
Km_pp2a_RyR = 4.1μM
## Fitted PKACI, PKACII, and PP1 parameters
PKACI_basal = 0.0734 ## basal activity
PKACI_activated = 0.1995
PKACI_KM = 0.0139μM
Expand All @@ -319,72 +271,60 @@ function get_bar_sys_reduced(ISO=0μM; name=:bar_sys)
PP1_basal = 0.8927
PP1_activated = 0.0492
PP1_KI = 0.00637μM
PLBp_basal = 0.0824
PLBp_activated = 0.7961
PLBp_KM = 0.00597μM
PLBp_nHill = 1.8167
PLMp_basal = 0.1172
PLMp_activated = 0.6645
PLMp_KM = 0.00823μM
PLMp_nHill = 1.35784
TnIp_basal = 0.0669
TnIp_activated = 0.7524
TnIp_KM = 0.007913μM
TnIp_nHill = 1.6736
LCCap_basal = 0.2205
LCCap_activated = 0.2339
LCCap_KM = 0.00726μM
LCCap_nHill = 0.9932
LCCbp_basal = 0.2517
LCCbp_activated = 0.2461
LCCbp_KM = 0.00695μM
LCCbp_nHill = 0.9934
KURp_basal = 0.4390
KURp_activated = 0.2563
KURp_KM = 0.00557μM
KURp_nHill = 0.9931
RyRp_basal = 0.2054
RyRp_activated = 0.2399
RyRp_KM = 0.0075135μM
end

vs = @variables begin
PLBp(t)
PLMp(t)
TnIp(t)
LCCap(t)
LCCbp(t)
KURp(t)
RyRp(t) = 0μM
PLB(t)
PLM(t)
TnI(t)
LCCa(t)
LCCb(t)
KURn(t)
RyRn(t)
LCCa_PKAp(t)
LCCb_PKAp(t)
fracPKACI(t)
fracPKACII(t)
fracPP1(t)
fracPLBp(t)
fracPLMp(t)
TnI_PKAp(t)
IKUR_PKAp(t)
RyR_PKAp(t)
PKACI(t)
PKACII(t)
PP1(t)
end

## Solve x for Vf * x / (x + k1) = Vr * (1 - x) / (1 - x + k2)
## x is the steady-state phosphorylation fraction
function _phos_fraction(Vf, Vr, k1, k2)
A = 1 - (Vf / Vr)
B = (Vf / Vr) * (1 + k2) + (k1 - 1)
C = -k1
n = (-B + sqrt(B^2 - 4 * A * C)) / 2A
return 1 - n
end

## Fitted activities
eqs = [
PLBtotBA ~ PLB + PLBp,
PLMtotBA ~ PLM + PLMp,
TnItotBA ~ TnI + TnIp,
LCCtotBA ~ LCCa + LCCap,
LCCtotBA ~ LCCb + LCCbp,
IKurtotBA ~ KURn + KURp,
RyRtotBA ~ RyRn + RyRp,
PLBp ~ PLBtotBA * fracPLBp,
PLMp ~ PLMtotBA * fracPLMp,
TnIp ~ TnItotBA * TnI_PKAp,
LCCap ~ LCCtotBA * LCCa_PKAp,
LCCbp ~ LCCtotBA * LCCb_PKAp,
KURp ~ IKurtotBA * IKUR_PKAp,
RyR_PKAp ~ RyRp / RyRtotBA,
## Fitted activities
PKACI ~ PKACItot * (PKACI_basal + PKACI_activated * hil(ISO, PKACI_KM)),
PKACII ~ PKACIItot * (PKACII_basal + PKACII_activated * hil(ISO, PKACII_KM)),
PP1 ~ PP1totBA * (PP1_basal + PP1_activated * hilr(ISO, PP1_KI)),
fracPLBp ~ _phos_fraction(k_PKA_PLB * PKACI, k_PP1_PLB * PP1, Km_PKA_PLB / PLBtotBA, Km_PP1_PLB / PLBtotBA),
fracPLMp ~ _phos_fraction(k_PKA_PLM * PKACI, k_PP1_PLM * PP1, Km_PKA_PLM / PLMtotBA, Km_PP1_PLM / PLMtotBA),
TnI_PKAp ~ _phos_fraction(k_PKA_TnI * PKACI, k_PP2A_TnI * PP2A_TnI, Km_PKA_TnI / TnItotBA, Km_PP2A_TnI / TnItotBA),
LCCa_PKAp ~ _phos_fraction(k_PKA_LCC * (PKACII_LCCtotBA / PKAIItot) * PKACII, k_PP2A_LCC * PP2A_LCC, Km_PKA_LCC / LCCtotBA / epsilon, Km_PP2A_LCC / LCCtotBA / epsilon),
LCCb_PKAp ~ _phos_fraction(k_PKA_LCC * (PKACII_LCCtotBA / PKAIItot) * PKACII, k_PP1_LCC * PP1_LCC, Km_PKA_LCC / LCCtotBA / epsilon, Km_PP1_LCC / LCCtotBA / epsilon),
IKUR_PKAp ~ _phos_fraction(k_pka_KUR * (PKAII_KURtot / PKAIItot) * PKACII, PP1_KURtot * k_pp1_KUR, Km_pka_KUR / IKurtotBA / epsilon, Km_pp1_KUR / IKurtotBA / epsilon),
# _phos_fraction is not applicable to RyR
D(RyRp) ~ kcat_pka_RyR * (PKAII_RyRtot / PKAIItot) * PKACII * hil(RyRn * epsilon, Km_pka_RyR) - kcat_pp1_RyR * PP1_RyR * hil(RyRp * epsilon, Km_pp1_RyR) - kcat_pp2a_RyR * PP2A_RyR * hil(RyRp * epsilon, Km_pp2a_RyR)
fracPKACI ~ PKACI_basal + PKACI_activated * hil(ISO, PKACI_KM),
fracPKACII ~ PPKACII_basal + PKACII_activated * hil(ISO, PKACII_KM),
fracPP1 ~ PP1_basal + PP1_activated * hilr(ISO, PP1_KI),
fracPLBp ~ PLBp_basal + PLBp_activated * hil(ISO, PLBp_KM, PLBp_nHill),
fracPLMp ~ PLMp_basal + PLMp_activated * hil(ISO, PLMp_KM, PLMp_nHill),
TnI_PKAp ~ TnIp_basal + TnIp_activated * hil(ISO, TnIp_KM, TnIp_nHill),
LCCa_PKAp ~ LCCap_basal + LCCap_activated * hil(ISO, LCCap_KM, LCCap_nHill),
LCCb_PKAp ~ LCCbp_basal + LCCbp_activated * hil(ISO, LCCbp_KM, LCCbp_nHill),
IKUR_PKAp ~ KURp_basal + KURp_activated * hil(ISO, KURp_KM, KURp_nHill),
RyR_PKAp ~ RyRp_basal + RyRp_activated * hil(ISO, RyRp_KM),
]

return ODESystem(eqs, t; name)
Expand Down

0 comments on commit 5a33e3f

Please sign in to comment.