Skip to content

Commit

Permalink
Merge pull request #327 from SciML/fix_hill_via_symbolics
Browse files Browse the repository at this point in the history 
Implement system functions (Hill and michaelis-menten) via Symbolics.jl
  • Loading branch information
@TorkelE
TorkelE committed Apr 26, 2021
2 parents de16846 + 1a84043 commit 0c41fdd
Show file tree
Hide file tree
Showing 8 changed files with 264 additions and 64 deletions.
4 changes: 3 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,16 +11,18 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

[compat]
DataStructures = "0.18"
DocStringExtensions = "0.8"
Latexify = "0.14, 0.15"
MacroTools = "0.5.5"
ModelingToolkit = "5.14.0"
ModelingToolkit = "5.15.0"
Parameters = "0.12"
Reexport = "0.2, 1.0"
Requires = "1.0"
Symbolics = "0.1.22"
julia = "1.5"

[extras]
Expand Down
4 changes: 4 additions & 0 deletions src/Catalyst.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ const DEFAULT_IV = (@parameters t)[1]
@reexport using ModelingToolkit
import MacroTools
import Base: (==), merge!, merge
using Symbolics
using Latexify, Requires

# as used in Catlab
Expand All @@ -33,6 +34,9 @@ include("reaction_network.jl")
# reaction network macro
export @reaction_network, @add_reactions

# registers CRN specific functions using Symbolics.jl
include("registered_functions.jl")

# functions to query network properties
include("networkapi.jl")
export species, params, reactions, speciesmap, paramsmap, numspecies, numreactions, numparams
Expand Down
23 changes: 19 additions & 4 deletions src/latexify_recipes.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,17 @@
# Implements handling of registered functions.
mm_names = ([:mm, :Mm, :MM])
mmR_names = ([:mmR, :MmR, :MMR])
hill_names = ([:hill, :Hill])
hillR_names = ([:hillR, :HillR])
hillAR_names = ([:hillC, :hillAR, :HillC, :HillAR])

make_mm_exp(expr::Expr) = :($(expr.args[3])*$(expr.args[2])/($(expr.args[4])+$(expr.args[2])))
make_mmR_exp(expr::Expr) = :($(expr.args[3])*$(expr.args[4])/($(expr.args[4])+$(expr.args[2])))
make_hill_exp(expr::Expr) = :($(expr.args[3])*($(expr.args[2])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))
make_hillR_exp(expr::Expr) = :($(expr.args[3])*($(expr.args[4])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))
make_hillAR_exp(expr::Expr) = :($(expr.args[4])*($(expr.args[2])^$(expr.args[6]))/($(expr.args[5])^$(expr.args[6])+$(expr.args[2])^$(expr.args[6])+$(expr.args[3])^$(expr.args[6])))


#Recursively traverses an expression and removes things like X^1, 1*X. Will not actually have any effect on the expression when used as a function, but will make it much easier to look at it for debugging, as well as if it is transformed to LaTeX code.
function recursive_clean!(expr)
(expr isa Symbol) && (expr == :no___noise___scaling) && (return 1)
Expand Down Expand Up @@ -27,10 +41,11 @@ function recursive_clean!(expr)
end
if expr.head == :call
(expr.args[1] == :/) && (expr.args[3] == 1) && (return expr.args[2])
in(expr.args[1],hill_name) && return hill(expr)
in(expr.args[1],hillR_name) && return hillR(expr)
in(expr.args[1],mm_name) && return mm(expr)
in(expr.args[1],mmR_name) && return mmR(expr)
in(expr.args[1],mm_names) && return make_mm_exp(expr)
in(expr.args[1],mmR_names) && return make_mmR_exp(expr)
in(expr.args[1],hill_names) && return make_hill_exp(expr)
in(expr.args[1],hillR_names) && return make_hillR_exp(expr)
in(expr.args[1],hillAR_names) && return make_hillAR_exp(expr)
(expr.args[1] == :binomial) && (expr.args[3] == 1) && return expr.args[2]
#@isdefined($(expr.args[1])) || error("Function $(expr.args[1]) not defined.")
end
Expand Down
16 changes: 0 additions & 16 deletions src/reaction_network.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,25 +268,9 @@ function recursive_expand_functions!(expr::ExprValues)
(typeof(expr)!=Expr) && (return expr)
foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]), 1:length(expr.args))
if expr.head == :call
in(expr.args[1],hill_name) && return hill(expr)
in(expr.args[1],hillR_name) && return hillR(expr)
in(expr.args[1],mm_name) && return mm(expr)
in(expr.args[1],mmR_name) && return mmR(expr)

!isdefined(Catalyst,expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
end
return expr
end

#Hill function made avaiable (activation and repression).
hill_name = Set{Symbol}([:hill, :Hill, :h, :H, :HILL])
hill(expr::Expr) = :($(expr.args[3])*($(expr.args[2])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))
hillR_name = Set{Symbol}([:hill_repressor, :hillr, :hillR, :HillR, :hR, :hR, :Hr, :HR, :HILLR])
hillR(expr::Expr) = :($(expr.args[3])*($(expr.args[4])^$(expr.args[5]))/($(expr.args[4])^$(expr.args[5])+$(expr.args[2])^$(expr.args[5])))

#Michaelis-Menten function made available (activation and repression).
mm_name = Set{Symbol}([:MM, :mm, :Mm, :mM, :M, :m])
mm(expr::Expr) = :($(expr.args[3])*$(expr.args[2])/($(expr.args[4])+$(expr.args[2])))
mmR_name = Set{Symbol}([:mm_repressor, :MMR, :mmr, :mmR, :MmR, :mMr, :MR, :mr, :Mr, :mR])
mmR(expr::Expr) = :($(expr.args[3])*$(expr.args[4])/($(expr.args[4])+$(expr.args[2])))

106 changes: 106 additions & 0 deletions src/registered_functions.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
using Symbolics

# Registers the Michaelis-Menten function.
mm(X,v,K) = v*X / (X + K)
Mm(X,v,K) = v*X / (X + K)
MM(X,v,K) = v*X / (X + K)

@register mm(X,v,K); @register Mm(X,v,K); @register MM(X,v,K);

Symbolics.derivative(::typeof(mm), args::NTuple{3,Any}, ::Val{1}) = (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(mm), args::NTuple{3,Any}, ::Val{2}) = args[1]/(args[1]+args[3])
Symbolics.derivative(::typeof(mm), args::NTuple{3,Any}, ::Val{3}) = - args[2]*args[1]/(args[1]+args[3])^2

Symbolics.derivative(::typeof(Mm), args::NTuple{3,Any}, ::Val{1}) = (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(Mm), args::NTuple{3,Any}, ::Val{2}) = args[1]/(args[1]+args[3])
Symbolics.derivative(::typeof(Mm), args::NTuple{3,Any}, ::Val{3}) = - args[2]*args[1]/(args[1]+args[3])^2

Symbolics.derivative(::typeof(MM), args::NTuple{3,Any}, ::Val{1}) = (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(MM), args::NTuple{3,Any}, ::Val{2}) = args[1]/(args[1]+args[3])
Symbolics.derivative(::typeof(MM), args::NTuple{3,Any}, ::Val{3}) = - args[2]*args[1]/(args[1]+args[3])^2

# Registers the repressing Michaelis-Menten function.
mmR(X,v,K) = v*K / (X + K)
MmR(X,v,K) = v*K / (X + K)
MMR(X,v,K) = v*K / (X + K)

@register mmR(X,v,K); @register MmR(X,v,K); @register MMR(X,v,K);

Symbolics.derivative(::typeof(mmR), args::NTuple{3,Any}, ::Val{1}) = - (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(mmR), args::NTuple{3,Any}, ::Val{2}) = args[3]/(args[1]+args[3])
Symbolics.derivative(::typeof(mmR), args::NTuple{3,Any}, ::Val{3}) = args[2]*args[1]/(args[1]+args[3])^2

Symbolics.derivative(::typeof(MmR), args::NTuple{3,Any}, ::Val{1}) = - (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(MmR), args::NTuple{3,Any}, ::Val{2}) = args[3]/(args[1]+args[3])
Symbolics.derivative(::typeof(MmR), args::NTuple{3,Any}, ::Val{3}) = args[2]*args[1]/(args[1]+args[3])^2

Symbolics.derivative(::typeof(MMR), args::NTuple{3,Any}, ::Val{1}) = - (args[2]*args[3]) / (args[1]+args[3])^2
Symbolics.derivative(::typeof(MMR), args::NTuple{3,Any}, ::Val{2}) = args[3]/(args[1]+args[3])
Symbolics.derivative(::typeof(MMR), args::NTuple{3,Any}, ::Val{3}) = args[2]*args[1]/(args[1]+args[3])^2


# Registers the Hill function.
hill(X,v,K,n) = v*(X^n) / (X^n + K^n)
Hill(X,v,K,n) = v*(X^n) / (X^n + K^n)

@register hill(X,v,K,n); @register Hill(X,v,K,n);

Symbolics.derivative(::typeof(hill), args::NTuple{4,Any}, ::Val{1}) = args[2] * args[4] * (args[3]^args[4]) * (args[1]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(hill), args::NTuple{4,Any}, ::Val{2}) = (args[1]^args[4]) / (args[1]^args[4] + args[3]^args[4])
Symbolics.derivative(::typeof(hill), args::NTuple{4,Any}, ::Val{3}) = - args[2] * args[4] * (args[1]^args[4]) * (args[3]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(hill), args::NTuple{4,Any}, ::Val{4}) = args[2] * (args[1]^args[4]) * (args[3]^args[4]) * (log(args[1])-log(args[3])) / (args[1]^args[4] + args[3]^args[4])^2

Symbolics.derivative(::typeof(Hill), args::NTuple{4,Any}, ::Val{1}) = args[2] * args[4] * (args[3]^args[4]) * (args[1]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(Hill), args::NTuple{4,Any}, ::Val{2}) = (args[1]^args[4]) / (args[1]^args[4] + args[3]^args[4])
Symbolics.derivative(::typeof(Hill), args::NTuple{4,Any}, ::Val{3}) = - args[2] * args[4] * (args[1]^args[4]) * (args[3]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(Hill), args::NTuple{4,Any}, ::Val{4}) = args[2] * (args[1]^args[4]) * (args[3]^args[4]) * (log(args[1])-log(args[3])) / (args[1]^args[4] + args[3]^args[4])^2


# Registers the repressing hill function (alterantive to negative n).
hillR(X,v,K,n) = v*(K^n) / (X^n + K^n)
HillR(X,v,K,n) = v*(K^n) / (X^n + K^n)

@register hillR(X,v,K,n); @register HillR(X,v,K,n);

Symbolics.derivative(::typeof(hillR), args::NTuple{4,Any}, ::Val{1}) = - args[2] * args[4] * (args[3]^args[4]) * (args[1]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(hillR), args::NTuple{4,Any}, ::Val{2}) = (args[3]^args[4]) / (args[1]^args[4] + args[3]^args[4])
Symbolics.derivative(::typeof(hillR), args::NTuple{4,Any}, ::Val{3}) = args[2] * args[4] * (args[1]^args[4]) * (args[3]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(hillR), args::NTuple{4,Any}, ::Val{4}) = args[2] * (args[1]^args[4]) * (args[3]^args[4]) * (log(args[3])-log(args[1])) / (args[1]^args[4] + args[3]^args[4])^2

Symbolics.derivative(::typeof(HillR), args::NTuple{4,Any}, ::Val{1}) = - args[2] * args[4] * (args[3]^args[4]) * (args[1]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(HillR), args::NTuple{4,Any}, ::Val{2}) = (args[3]^args[4]) / (args[1]^args[4] + args[3]^args[4])
Symbolics.derivative(::typeof(HillR), args::NTuple{4,Any}, ::Val{3}) = args[2] * args[4] * (args[1]^args[4]) * (args[3]^(args[4]-1)) / (args[1]^args[4] + args[3]^args[4])^2
Symbolics.derivative(::typeof(HillR), args::NTuple{4,Any}, ::Val{4}) = args[2] * (args[1]^args[4]) * (args[3]^args[4]) * (log(args[3])-log(args[1])) / (args[1]^args[4] + args[3]^args[4])^2


# Registers the activation/repressing hill function.
hillC(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)
hillAR(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)
HillC(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)
HillAR(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)

@register hillC(X,Y,v,K,n); @register hillAR(X,Y,v,K,n); @register HillC(X,Y,v,K,n); @register HillAR(X,Y,v,K,n);

Symbolics.derivative(::typeof(hillC), args::NTuple{5,Any}, ::Val{1}) = args[3] * args[5] * (args[1]^(args[5]-1)) * (args[2]^args[5]+args[4]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillC), args::NTuple{5,Any}, ::Val{2}) = - args[3] * args[5] * (args[2]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillC), args::NTuple{5,Any}, ::Val{3}) = (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])
Symbolics.derivative(::typeof(hillC), args::NTuple{5,Any}, ::Val{4}) = - args[3] * args[5] * (args[3]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillC), args::NTuple{5,Any}, ::Val{5}) = args[3] * (args[1]^args[5]) * (log(args[1])*(args[2]^args[5] + args[4]^args[5]) - (args[2]^args[5])*log(args[2]) - (args[4]^args[5])*log(args[4])) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2

Symbolics.derivative(::typeof(hillAR), args::NTuple{5,Any}, ::Val{1}) = args[3] * args[5] * (args[1]^(args[5]-1)) * (args[2]^args[5]+args[4]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillAR), args::NTuple{5,Any}, ::Val{2}) = - args[3] * args[5] * (args[2]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillAR), args::NTuple{5,Any}, ::Val{3}) = (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])
Symbolics.derivative(::typeof(hillAR), args::NTuple{5,Any}, ::Val{4}) = - args[3] * args[5] * (args[3]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(hillAR), args::NTuple{5,Any}, ::Val{5}) = args[3] * (args[1]^args[5]) * (log(args[1])*(args[2]^args[5] + args[4]^args[5]) - (args[2]^args[5])*log(args[2]) - (args[4]^args[5])*log(args[4])) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2

Symbolics.derivative(::typeof(HillC), args::NTuple{5,Any}, ::Val{1}) = args[3] * args[5] * (args[1]^(args[5]-1)) * (args[2]^args[5]+args[4]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillC), args::NTuple{5,Any}, ::Val{2}) = - args[3] * args[5] * (args[2]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillC), args::NTuple{5,Any}, ::Val{3}) = (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])
Symbolics.derivative(::typeof(HillC), args::NTuple{5,Any}, ::Val{4}) = - args[3] * args[5] * (args[3]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillC), args::NTuple{5,Any}, ::Val{5}) = args[3] * (args[1]^args[5]) * (log(args[1])*(args[2]^args[5] + args[4]^args[5]) - (args[2]^args[5])*log(args[2]) - (args[4]^args[5])*log(args[4])) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2

Symbolics.derivative(::typeof(HillAR), args::NTuple{5,Any}, ::Val{1}) = args[3] * args[5] * (args[1]^(args[5]-1)) * (args[2]^args[5]+args[4]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillAR), args::NTuple{5,Any}, ::Val{2}) = - args[3] * args[5] * (args[2]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillAR), args::NTuple{5,Any}, ::Val{3}) = (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])
Symbolics.derivative(::typeof(HillAR), args::NTuple{5,Any}, ::Val{4}) = - args[3] * args[5] * (args[3]^(args[5]-1)) * (args[1]^args[5]) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
Symbolics.derivative(::typeof(HillAR), args::NTuple{5,Any}, ::Val{5}) = args[3] * (args[1]^args[5]) * (log(args[1])*(args[2]^args[5] + args[4]^args[5]) - (args[2]^args[5])*log(args[2]) - (args[4]^args[5])*log(args[4])) / (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
130 changes: 111 additions & 19 deletions test/custom_functions.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
### Fetch required packages and reaction networks ###
using DiffEqBase, Catalyst, Random, Test
using ModelingToolkit: get_states, get_ps

Expand All @@ -11,32 +10,125 @@ new_poly(x,p1,p2) = .5*p1*x^2
new_exp(x,p) = exp(-p*x)

custom_function_network_1 = @reaction_network begin
hill(X,v1,K1,2), X + Y --> Z1
mm(X,v2,K2), X + Y --> Z2
p1*X^2+p2, X + Y --> Z3
exp(-p3*Y), X + Y --> Z4
hillR(X,v3,K3,2), X + Y --> Z5
mmR(X,v4,K4), X + Y --> Z6
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4
hill(X1,v1,K1,2), X1 + Y1--> Z1
mm(X2,v2,K2), X2 + Y2 --> Z2
p1*X3^2+p2, X3 + Y3 --> Z3
exp(-p3*Y4), X4 + Y4 --> Z4
hillR(X5,v3,K3,2), X5 + Y5 --> Z5
mmR(X6,v4,K4), X6 + Y6 --> Z6
hillAR(X7,Y7,v5,K5,2), X7 + Y7 --> Z7
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4 v5 K5

custom_function_network_2 = @reaction_network begin
new_hill(X,v1,K1,2), X + Y --> Z1
v2*X/(X+K2), X + Y --> Z2
2*new_poly(X,p1,p2)+p2, X + Y --> Z3
new_exp(Y,p3), X + Y --> Z4
v3*(K3^2)/(K3^2+X^2), X + Y --> Z5
v4*K4/(X+K4), X + Y --> Z6
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4
new_hill(X1,v1,K1,2), X1 + Y1 --> Z1
v2*X2/(X2+K2), X2 + Y2 --> Z2
2*new_poly(X3,p1,p2)+p2, X3 + Y3 --> Z3
new_exp(Y4,p3), X4 + Y4--> Z4
v3*(K3^2)/(K3^2+X5^2), X5 + Y5 --> Z5
v4*K4/(X6+K4), X6 + Y6 --> Z6
v5*(X7^2)/(K5^2+X7^2+Y7^2), X7 + Y7 --> Z7
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4 v5 K5

f1 = ODEFunction(convert(ODESystem,custom_function_network_1),jac=true)
f2 = ODEFunction(convert(ODESystem,custom_function_network_2),jac=true)
g1 = SDEFunction(convert(SDESystem,custom_function_network_1))
g2 = SDEFunction(convert(SDESystem,custom_function_network_2))
for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2]
u0 = factor*rand(rng,length(get_states(custom_function_network_1)))
p = factor*rand(rng,length(get_ps(custom_function_network_2)))
t = rand(rng)
@test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 100*eps())
@test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 100*eps())
@test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 100*eps())
@test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 10e-10)
@test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 10e-10)
@test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 10e-10)
end

### Tests that the various notations gives identical results ###

# Michaelis-Menten function.
mm_network = @reaction_network begin
(1.,1.), 0 ↔ X
mm(X,v,K), 0 --> X1
Mm(X,v,K), 0 --> X2
MM(X,v,K), 0 --> X3
end v K
f_mm = ODEFunction(convert(ODESystem,mm_network),jac=true)

u0 = 10*rand(rng,length(get_states(mm_network)))
p = 10*rand(rng,length(get_ps(mm_network)))
t = 10*rand(rng)

f_mm_output = f_mm(u0,p,t)[2:end]
f_mm_jac_output = f_mm.jac(u0,p,t)[2:end,1]
@test (maximum(f_mm_output) - minimum(f_mm_output)) .< 100*eps()
@test (maximum(f_mm_jac_output) - minimum(f_mm_jac_output)) .< 100*eps()

# Repressing Michaelis-Menten function.
mmR_network = @reaction_network begin
(1.,1.), 0 ↔ X
mmR(X,v,K), 0 --> X1
MmR(X,v,K), 0 --> X2
MMR(X,v,K), 0 --> X3
end v K
f_mmR = ODEFunction(convert(ODESystem,mmR_network),jac=true)

u0 = 10*rand(rng,length(get_states(mmR_network)))
p = 10*rand(rng,length(get_ps(mmR_network)))
t = 10*rand(rng)

f_mmR_output = f_mmR(u0,p,t)[2:end]
f_mmR_jac_output = f_mmR.jac(u0,p,t)[2:end,1]
@test (maximum(f_mmR_output) - minimum(f_mmR_output)) .< 100*eps()
@test (maximum(f_mmR_jac_output) - minimum(f_mmR_jac_output)) .< 100*eps()

# Hill function.
hill_network = @reaction_network begin
(1.,1.), 0 ↔ X
hill(X,v,K,2), 0 --> X1
Hill(X,v,K,2), 0 --> X2
end v K
f_hill = ODEFunction(convert(ODESystem,hill_network),jac=true)

u0 = 10*rand(rng,length(get_states(hill_network)))
p = 10*rand(rng,length(get_ps(hill_network)))
t = 10*rand(rng)

f_hill_output = f_hill(u0,p,t)[2:end]
f_hill_jac_output = f_hill.jac(u0,p,t)[2:end,1]
@test (maximum(f_hill_output) - minimum(f_hill_output)) .< 100*eps()
@test (maximum(f_hill_jac_output) - minimum(f_hill_jac_output)) .< 100*eps()

# Repressing Hill function.
hillR_network = @reaction_network begin
(1.,1.), 0 ↔ X
hillR(X,v,K,2), 0 --> X1
HillR(X,v,K,2), 0 --> X2
end v K
f_hillR = ODEFunction(convert(ODESystem,hillR_network),jac=true)

u0 = 10*rand(rng,length(get_states(hillR_network)))
p = 10*rand(rng,length(get_ps(hillR_network)))
t = 10*rand(rng)

f_hillR_output = f_hillR(u0,p,t)[2:end]
f_hillR_jac_output = f_hillR.jac(u0,p,t)[2:end,1]
@test (maximum(f_hillR_output) - minimum(f_hillR_output)) .< 100*eps()
@test (maximum(f_hillR_jac_output) - minimum(f_hillR_jac_output)) .< 100*eps()

# Activation/repressing Hill function.
hillAR_network = @reaction_network begin
(1.,1.), 0 ↔ (X,Y)
hillC(X,Y,v,K,2), 0 --> X1
HillC(X,Y,v,K,2), 0 --> X2
hillAR(X,Y,v,K,2), 0 --> X3
HillAR(X,Y,v,K,2), 0 --> X4
end v K
f_hillAR = ODEFunction(convert(ODESystem,hillAR_network),jac=true)

u0 = 10*rand(rng,length(get_states(hillAR_network)))
p = 10*rand(rng,length(get_ps(hillAR_network)))
t = 10*rand(rng)

f_hillAR_output = f_hillAR(u0,p,t)[3:end]
f_hillAR_jac_output = f_hillAR.jac(u0,p,t)[3:end,1]
@test (maximum(f_hillAR_output) - minimum(f_hillAR_output)) .< 100*eps()
@test (maximum(f_hillAR_jac_output) - minimum(f_hillAR_jac_output)) .< 100*eps()
Loading

0 comments on commit 0c41fdd

Please sign in to comment.