Learn practical skills, build real-world projects, and advance your career
#Load all the packages.
using DifferentialEquations,Plots,DiffEqBiological
using Latexify
#constants
global growthrate = 1. #growth rate
global N_0 = 10^20 #population cap


@reaction_func growth(x) = growthrate*x*(1 - (A+B+B1+C)/N_0)


model = @reaction_network begin
    

    
    #model equations
    growth(A), 0 → A
    growth(B), 0 → B
    growth(B1), 0 → B1
    growth(C), 0 → C
    
    0.01, A → B
    0.02, B → A
    
     0.01, A → B1
    0.02, B1 → A
    
    0.001,A → 0
    0.001,B → G
    0.001,C → 0
    0.001,B1 → G
    0.002,G → 0
    
    0.05, C + G → B1
    
    end

┌ Warning: The RegularJump interface has changed to be matrix-free. See the documentation for more details. └ @ DiffEqJump /home/ved/.julia/packages/DiffEqJump/sEB5p/src/jumps.jl:43
(::reaction_network) (generic function with 2 methods)
latexify(model,env=:chemical)
\begin{align} \require{mhchem} \ce{ \varnothing &->[growthrate \cdot A \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right)] A}\\ \ce{ \varnothing &->[growthrate \cdot B \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right)] B}\\ \ce{ \varnothing &->[growthrate \cdot B1 \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right)] B1}\\ \ce{ \varnothing &->[growthrate \cdot C \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right)] C}\\ \ce{ A &<=>[0.01][0.02] B}\\ \ce{ A &<=>[0.01][0.02] B1}\\ \ce{ A &->[0.001] \varnothing}\\ \ce{ B &->[0.001] G}\\ \ce{ C &->[0.001] \varnothing}\\ \ce{ B1 &->[0.001] G}\\ \ce{ G &->[0.002] \varnothing}\\ \ce{ C + G &->[0.05] B1} \end{align}
latexify(model)
\begin{align} \frac{dA(t)}{dt} =& growthrate \cdot A \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right) - 0.01 \cdot A + 0.02 \cdot B - 0.01 \cdot A + 0.02 \cdot B1 - 0.001 \cdot A \\ \frac{dB(t)}{dt} =& growthrate \cdot B \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right) + 0.01 \cdot A - 0.02 \cdot B - 0.001 \cdot B \\ \frac{dB1(t)}{dt} =& growthrate \cdot B1 \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right) + 0.01 \cdot A - 0.02 \cdot B1 - 0.001 \cdot B1 + 0.05 \cdot C \cdot G \\ \frac{dC(t)}{dt} =& growthrate \cdot C \cdot \left( 1 - \frac{A + B + B1 + C}{N_{0}} \right) - 0.001 \cdot C - 0.05 \cdot C \cdot G \\ \frac{dG(t)}{dt} =& 0.001 \cdot B + 0.001 \cdot B1 - 0.002 \cdot G - 0.05 \cdot C \cdot G \end{align}
tspan = (0.,500.)
initial = [0.,1.0 * 10^6,0.,1.0*10^4,0]
prob = ODEProblem(model,initial,tspan)
sol = solve(prob,maxiters = 1e10,saveat = 0.1);

timepoints = [0.0:0.1:500.0]


a = [sol[i][1] for i = 1:length(sol)];
b = [sol[i][2] for i = 1:length(sol)];
b1 = [sol[i][3] for i = 1:length(sol)];
c = [sol[i][4] for i = 1:length(sol)];
g = [sol[i][5] for i = 1:length(sol)];