Homework 5
Env setup
"$VERSION"
using Distributions, StatsBase, Plots, DiffEqBiological, DifferentialEquations
using InteractiveUtils
About @reaction_network
By default, the propensity is multiplied by the Number of reactants except for synthesis reactions (e.g. 0 --> A)
Parameters are listed here (without comma) at the end of
@reaction_network
, without commaThe initial conditions correspond to A, R, C, so it should be an array with a length of 3
I made the time span longer so you can see the (irregular) oscillations
# A + R -> C model
AR_model = begin
(γA/bA)*((α0+(A/KA))/(1+(A/KA))), 0 --> A
(γR/bR)*(((A/KR))/(1+(A/KR))), 0 --> R
δA, A --> 0
δR, R --> 0
kC, A + R --> C
δA, C --> R
end γA bA α0 KA γR bR KR δA δR kC
# Parameter list
p = (γA = 250, bA = 5, α0 = 0.1, KA = 0.5, γR = 50, bR = 10, KR = 1, δA = 1, δR = 0.1, kC = 100)
u0 = [0, 10, 35]
tspan = (0.0, 200.0)
# Not the problem in stochastic methods are divided in two steps
prob = DiscreteProblem(u0, tspan, p)
jumpProb = JumpProblem(prob, Direct(), AR_model)
sol = solve(jumpProb, SSAStepper());
plot(sol, lw=1, label=["A" "R" "C"])
# Remember to convert the data type of state variable to real number (Float64) for the ODE problem
u0f = convert.(Float64, u0)
odeprob = ODEProblem(AR_model, u0f, tspan, p);
odesol = solve(odeprob, progress=true);
plot(odesol, label=["A" "R" "C"])
As for DIY stochastic solver,
NA
means the number of A, not the rate (propensity) of synthesis of A,aA
. The program will tell you it cannot find the definition of NA since NA is dependent on itself in this line:NA = (γA/bA)*((α0+(NA/KA))/(1+(NA/KA)))
And NA is part of the state variable,
u
. Thus,NA = u[1], NR = u[2]
, and so on.The state variable,
u
, represent three components of A, R, and C. Check iflength(u) == 3
The reaction stoichiometry should match the reactions in this problem
ss = ([-1 1], [1 -1])
is for example of A<-->B only.
#=
Stochastic chemical reaction: Relaxation oscillator by José Vilar
MMSB 7.8.27 page 270.
=#
# Moved u0, tspan, Reaction stoichiometry on top of the list
function ssa(u0, tspan, ss; γA = 250, bA = 5, α0 = 0.1, KA = 0.5, γR = 50, bR = 10, KR = 1, δA = 1, δR = 0.1, kC = 100) # parameters u= current state
t = tspan[1] # time counter
ts = [t] # time recorder
u = copy(u0) # state variable counter
us = copy(u0) # state variable recorder
while t < tspan[2]
#
NA, NR, NC = u
aA = (γA/bA)*((α0+(NA/KA))/(1+(NA/KA))) # propensity of reaction: -> bA^A (Activator Synthetis)
aR = (γR/bR)*(((NA/KR))/(1+(NA/KR))) # propensity of reaction: -> bR^R (Repressor Synthetis)
Ad = δA*NA # propensity of reaction: A -> (Activator Decay)
Rd = δR*NR # propensity of reaction: R -> (Repressor Decay)
aC = kC*NA*NR #propensity of reaction: A + R -> C (Association)
Cd = δA*NC # propensity of reaction: C -> R (Dissociation and Decay)
aTot = aA + aR + Ad + Rd + aC + Cd # Total propensity
# Exponentially distributed waiting time
dt = rand(Exponential(1 / aTot))
# Draw which reaction to occur
draw = sample(Weights([aA, aR, Ad, Rd, aC, Cd]))
# Update state variables (u) and time (t)
u .+= ss[draw]
t += dt
# Record u and t
us = [us u]
ts = [ts; t]
end
# transpose the us to match the convention of one column = one series
return (ts, us')
end
u0 = [0, 10, 35]
tspan = (0.0, 100.0)
ts, us = ssa(u0, tspan, ss)
plot(ts, us, label=["NA" "NR" "NC"])