Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement system functions (Hill and michaelis-menten) via Symbolics.jl #327

Merged
merged 11 commits into from
Apr 26, 2021
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ 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"
Expand All @@ -21,6 +22,7 @@ ModelingToolkit = "5.14.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
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
16 changes: 0 additions & 16 deletions src/reaction_network.jl
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
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