Model descriptions#
The Mitochondrial Retrograde (RTG) Signalling Model is a system of 11 coupled ODEs.
using RetroSignalModel
using ModelingToolkit
@named sys = RtgMTK(; simplify=false)
\[\begin{split} \begin{align}
\frac{\mathrm{d} \mathrm{Rtg2A}_{c}\left( t \right)}{\mathrm{d}t} =& v(t)_1 - v(t)_2 \\
\frac{\mathrm{d} \mathrm{Rtg2Mks}_{c}\left( t \right)}{\mathrm{d}t} =& v(t)_2 \\
\frac{\mathrm{d} \mathrm{BmhMks}\left( t \right)}{\mathrm{d}t} =& v(t)_3 \\
\frac{\mathrm{d} \mathrm{Rtg13A}_{c}\left( t \right)}{\mathrm{d}t} =& - v(t)_4 + v(t)_7 \\
\frac{\mathrm{d} \mathrm{Rtg13I}_{c}\left( t \right)}{\mathrm{d}t} =& v(t)_4 + v(t)_8 \\
\frac{\mathrm{d} \mathrm{Rtg3A}_{c}\left( t \right)}{\mathrm{d}t} =& - v(t)_{1 2} + v(t)_5 - v(t)_7 \\
\frac{\mathrm{d} \mathrm{Rtg3A}_{n}\left( t \right)}{\mathrm{d}t} =& v(t)_{1 2} - v(t)_6 - v(t)_9 \\
\frac{\mathrm{d} \mathrm{Rtg3I}_{n}\left( t \right)}{\mathrm{d}t} =& - v(t)_{1 0} + v(t)_{1 3} + v(t)_6 \\
\frac{\mathrm{d} \mathrm{Rtg1}_{n}\left( t \right)}{\mathrm{d}t} =& - v(t)_{1 0} + v(t)_{1 1} - v(t)_9 \\
\frac{\mathrm{d} \mathrm{Rtg13A}_{n}\left( t \right)}{\mathrm{d}t} =& v(t)_9 \\
\frac{\mathrm{d} \mathrm{Rtg13I}_{n}\left( t \right)}{\mathrm{d}t} =& v(t)_{1 0} \\
s\left( t \right) =& 0 \\
v(t)_1 =& \frac{\left( mul_{S} s\left( t \right) \right)^{n_{S}} ksV \mathrm{Rtg2I}_{c}\left( t \right)}{\left( mul_{S} s\left( t \right) \right)^{n_{S}} + ksD^{n_{S}}} - k2I \mathrm{Rtg2A}_{c}\left( t \right) \\
v(t)_2 =& - kn2M \mathrm{Rtg2Mks}_{c}\left( t \right) + k2M \mathrm{Mks}\left( t \right) \mathrm{Rtg2A}_{c}\left( t \right) \\
v(t)_3 =& - knBM \mathrm{BmhMks}\left( t \right) + kBM \mathrm{Bmh}\left( t \right) \mathrm{Mks}\left( t \right) \\
v(t)_4 =& \left( k13I + \frac{k13IV \mathrm{BmhMks}\left( t \right)}{k13ID + \mathrm{BmhMks}\left( t \right)} \right) \mathrm{Rtg13A}_{c}\left( t \right) \\
v(t)_5 =& k3A_{c} \mathrm{Rtg3I}_{c}\left( t \right) - k3I_{c} \mathrm{Rtg3A}_{c}\left( t \right) \\
v(t)_6 =& k3I_{n} \mathrm{Rtg3A}_{n}\left( t \right) \\
v(t)_7 =& - kn13_{c} \mathrm{Rtg13A}_{c}\left( t \right) + k13_{c} \mathrm{Rtg1}_{c}\left( t \right) \mathrm{Rtg3A}_{c}\left( t \right) \\
v(t)_8 =& - kn13_{c} \mathrm{Rtg13I}_{c}\left( t \right) + k13_{c} \mathrm{Rtg3I}_{c}\left( t \right) \mathrm{Rtg1}_{c}\left( t \right) \\
v(t)_9 =& - kn13_{n} \mathrm{Rtg13A}_{n}\left( t \right) + k13_{n} \mathrm{Rtg3A}_{n}\left( t \right) \mathrm{Rtg1}_{n}\left( t \right) \\
v(t)_{1 0} =& - kn13_{n} \mathrm{Rtg13I}_{n}\left( t \right) + k13_{n} \mathrm{Rtg3I}_{n}\left( t \right) \mathrm{Rtg1}_{n}\left( t \right) \\
v(t)_{1 1} =& k1in \mathrm{Rtg1}_{c}\left( t \right) - k1out \mathrm{Rtg1}_{n}\left( t \right) \\
v(t)_{1 2} =& k3inA \mathrm{Rtg3A}_{c}\left( t \right) - k3outA \mathrm{Rtg3A}_{n}\left( t \right) \\
v(t)_{1 3} =& k3inI \mathrm{Rtg3I}_{c}\left( t \right) - k3outI \mathrm{Rtg3I}_{n}\left( t \right) \\
\mathrm{Rtg13}_{n}\left( t \right) =& \mathrm{Rtg13A}_{n}\left( t \right) + \mathrm{Rtg13I}_{n}\left( t \right) \\
\mathrm{Rtg13}_{c}\left( t \right) =& \mathrm{Rtg13I}_{c}\left( t \right) + \mathrm{Rtg13A}_{c}\left( t \right) \\
\mathrm{Rtg3}_{n}\left( t \right) =& \mathrm{Rtg13}_{n}\left( t \right) + \mathrm{Rtg3I}_{n}\left( t \right) + \mathrm{Rtg3A}_{n}\left( t \right) \\
\mathrm{Rtg3}_{c}\left( t \right) =& \mathrm{Rtg3I}_{c}\left( t \right) + \mathrm{Rtg3A}_{c}\left( t \right) + \mathrm{Rtg13}_{c}\left( t \right) \\
{\Sigma}Rtg1_{n}\left( t \right) =& \mathrm{Rtg13}_{n}\left( t \right) + \mathrm{Rtg1}_{n}\left( t \right) \\
{\Sigma}Rtg1_{c}\left( t \right) =& \mathrm{Rtg1}_{c}\left( t \right) + \mathrm{Rtg13}_{c}\left( t \right) \\
{\Sigma}Rtg1 =& \mathrm{Rtg13I}_{c}\left( t \right) + \mathrm{Rtg1}_{c}\left( t \right) + \mathrm{Rtg13A}_{c}\left( t \right) + \mathrm{Rtg13A}_{n}\left( t \right) + \mathrm{Rtg1}_{n}\left( t \right) + \mathrm{Rtg13I}_{n}\left( t \right) \\
{\Sigma}Rtg2 =& \mathrm{Rtg2Mks}_{c}\left( t \right) + \mathrm{Rtg2I}_{c}\left( t \right) + \mathrm{Rtg2A}_{c}\left( t \right) \\
{\Sigma}Rtg3 =& \mathrm{Rtg13I}_{c}\left( t \right) + \mathrm{Rtg3I}_{c}\left( t \right) + \mathrm{Rtg3I}_{n}\left( t \right) + \mathrm{Rtg13A}_{c}\left( t \right) + \mathrm{Rtg3A}_{n}\left( t \right) + \mathrm{Rtg13A}_{n}\left( t \right) + \mathrm{Rtg3A}_{c}\left( t \right) + \mathrm{Rtg13I}_{n}\left( t \right) \\
{\Sigma}Mks =& \mathrm{BmhMks}\left( t \right) + \mathrm{Rtg2Mks}_{c}\left( t \right) + \mathrm{Mks}\left( t \right) \\
{\Sigma}Bmh =& \mathrm{BmhMks}\left( t \right) + \mathrm{Bmh}\left( t \right)
\end{align}
\end{split}\]
Parameter estimation#
Simulated annealing (SAMIN) with bounds was used to find a set of parameters that fit exprerimental conditions. The following code is just a demonstration. In peactice we ran the code for much longer.
using RetroSignalModel
using Optim
optim_params(targetratio=3, optimoptions=Optim.Options(iterations=100, show_trace=true, show_every=10))
Iter Function value Gradient norm
0 4.771213e-01 NaN
* time: 14.417295932769775
10 4.771212e-01 NaN
* time: 19.610579013824463
20 7.503527e-01 NaN
* time: 23.518472909927368
30 9.499060e-01 NaN
* time: 27.427124977111816
40 4.885749e-01 NaN
* time: 31.3331880569458
50 6.239914e-01 NaN
* time: 35.16519498825073
60 5.767704e-01 NaN
* time: 39.068016052246094
70 4.962482e-01 NaN
* time: 43.00966000556946
80 4.199331e-01 NaN
* time: 46.942384004592896
90 4.283133e-01 NaN
* time: 50.758256912231445
================================================================================
SAMIN results
NO CONVERGENCE: MAXEVALS exceeded
Obj. value: 0.38137
parameter search width
95.94731 999.99900
4.80265 5.00000
5.37021 999.99900
429.77213 999.99900
6.01627 999.99900
867.13359 999.99900
37.79521 999.99900
53.86821 999.99900
51.31396 999.99900
514.77603 999.99900
195.73608 999.99900
255.56565 999.99900
337.39836 999.99900
5.76125 999.99900
232.27003 999.99900
722.48373 999.99900
85.34188 999.99900
368.15902 999.99900
784.47771 999.99900
332.50347 999.99900
757.54210 999.99900
251.63085 999.99900
612.34513 999.99900
187.25967 999.99900
================================================================================
(res = * Status: failure (reached maximum number of iterations)
* Candidate solution
Final objective value: 3.813664e-01
* Found with
Algorithm: SAMIN
* Convergence measures
|x - x'| = Inf ≰ 1.0e-06
|x - x'|/|x'| = NaN ≰ 0.0e+00
|f(x) - f(x')| = Inf ≰ 1.0e-12
|f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
|g(x)| = NaN ≰ 0.0e+00
* Work counters
Seconds run: 54 (vs limit Inf)
Iterations: 100
f(x) calls: 100
∇f(x) calls: 0
, parammap = Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}(k2I => 95.94731033304312, kn2M => 6.016274705090665, k2M => 867.1335921544978, k13ID => 514.7760348265062, kn13_c => 232.270025781651, k1out => 784.4777143398859, k3outI => 187.25966765242913, knBM => 37.79521150738747, kn13_n => 85.34188348844316, k3inI => 612.34513249146…))
Valid parameter sets#
using RetroSignalModel
using CSV
using DataFrames
using Plots
filename = joinpath(dirname(pathof(RetroSignalModel)), "data", "solution_rtgMTK_optim.csv")
dfoptim = let
df = CSV.read(filename, DataFrame)
df[!, Not(:n_S)] .= log10.(df[!, Not(:n_S)])
end
res = map(sort(names(dfoptim))) do lab
histogram(dfoptim[!, lab], label=lab)
end
plot(res..., layout=(6, 4), size=(1024, 1024))
# savefig("optimparams.png")
filename = joinpath(dirname(pathof(RetroSignalModel)), "data", "solution_rtgM4.csv")
dforig = let
df = CSV.read(filename, DataFrame)
df[!, Not(:n_S)] .= log10.(df[!, Not(:n_S)])
end
res = map(sort(names(dforig))) do lab
histogram(dforig[!, lab], label=lab)
end
plot(res..., layout=(6, 4), size=(1024, 1024))
# savefig("randomparams.png")