Homework 5

Env setup

Julia-1.4.1-DE
Download as Docker image from:
Copy
This image was imported from: docker.nextjournal.com/environment@sha256:8a2137d6683e4b488fb15c3ad5be0cc0a2ed276f6d9d4c66fb50c9f6414c2635
"$VERSION"
0.0s
Julia
Julia-1.4.1-DE
"1.4.1"
using Distributions, StatsBase, Plots, DiffEqBiological, DifferentialEquations
using InteractiveUtils
145.0s
Julia
Julia-1.4.1-DE

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 comma

  • The 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 = @reaction_network 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());
40.7s
Julia
Julia-1.4.1-DE
plot(sol, lw=1, label=["A" "R" "C"])
28.8s
Julia
Julia-1.4.1-DE
# Remember to convert the data type of state variable to real number (Float64) for the ODE problem
u0f = convert.(Float64, u0)
1.4s
Julia
Julia-1.4.1-DE
3-element Array{Float64,1}: 0.0 10.0 35.0
odeprob = ODEProblem(AR_model, u0f, tspan, p);
odesol = solve(odeprob, progress=true);
43.4s
Julia
Julia-1.4.1-DE
plot(odesol, label=["A" "R" "C"])
3.0s
Julia
Julia-1.4.1-DE

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 if length(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)
1.0s
Julia
Julia-1.4.1-DE
plot(ts, us, label=["NA" "NR" "NC"])
Shift+Enter to run
Julia
Julia-1.4.1-DE
Runtimes (1)