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

TorkelE · 2021-04-20T11:39:42Z

Previously, when the function recursive_expand_functions!() scanned the input expressions, it searched for a couple of keywords (like :hill), and if encountering an appropriate function call, replacing it with the correct expression (from a manually written list). This PR will:

Remove the old implementation.
Add a new implementation via Symbolics.
Catalyst.jl will now also have Symbolics.jl as a listed dependency.

There are now 5 functions specially implemented:

The Michaelis-Menten function, mm(X,v,K) = v*X/(X+K). Keywords: mm, Mm, and MM.
The repressive Michaelis-Menten function, mmR(X,v,K) = v*K/(X+K). Keywords: mmR, MmR, and MMR.
The Hill function, hill(X,v,K,n) = v*X^n/(X^n+K^n). Keywords: hill and Hill.
The repressive Hill function, hillR(X,v,K,n) = v*K^n/(X^n+K^n). Keywords: hillR and HillR.
The activating/repressing (or combined) Hill function, hillR(X,Y,v,K,n) = v*X^n/(X^n+Y^n+K^n). Keywords: hillAR, hillC, HillAR, and HillC.
The keywords are the recognised calls for the respective function.

The implementation looks like this (example for the Hill function):

# 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

(I used Wolfram Alpha to double check the derivatives)

It should note that this is currently a breaking change. Previously we allowed some more obscure function calls to be recognized (these were recognised for the repressing hill: ([:hill_repressor, :hillr, :hillR, :HillR, :hR, :hR, :Hr, :HR, :HILLR])). Initially I thought this change might implement these as normal function in people's environment, and figured it might be good to be conservative. I removed some options.

I don't think anyone uses these. But it would be possible to implement these using Symbolics, or as a compromise, implement them using the old methods. However, this might also be a breaking change we could accept to entirely scrap the old system.

TorkelE · 2021-04-20T12:30:19Z

The error seems to be due to two systems that should be "identical" no longer are so. This is basically the test we run:

using Catalyst
using ModelingToolkit

rn1 = @reaction_network begin
    mm(X2,v1,K1), ∅ → X1
    mm(X3,v2,K2), ∅ → X2
    (k1,k2), X1 ⟷ X3
    (k3,k4), X3 + X2 ⟷ X4 + X1
    d, (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 k1 k2 k3 k4 d

@parameters v1 K1 v2 K2 k1 k2 k3 k4 k5 p d t
@variables X1(t) X2(t) X3(t) X4(t) X5(t)

rxs_1 = [Reaction(v1*X2/(K1+X2), nothing, [X1], nothing, [1]),
       Reaction(v2*X3/(K2+X3), nothing, [X2], nothing, [1]),
       Reaction(k1, [X1], [X3], [1], [1]),
       Reaction(k2, [X3], [X1], [1], [1]),
       Reaction(k3, [X3,X2], [X4,X1], [1,1], [1,1]),
       Reaction(k4, [X4,X1], [X3,X2], [1,1], [1,1]),
       Reaction(d, [X1], nothing, [1], nothing),
       Reaction(d, [X2], nothing, [1], nothing),
       Reaction(d, [X3], nothing, [1], nothing),
       Reaction(d, [X4], nothing, [1], nothing)]
rn2 = ReactionSystem(rxs_1 , t, [X1,X2,X3,X4], [v1,K1,v2,K2,k1,k2,k3,k4,d])

rn1.eqs[1].rate    # = Catalyst.mm(X2(t), v1, K1)
rn2.eqs[1].rate    # = v1*X2(t)*((K1 + X"(t))^-1)

isequal(rn1.eqs[1].rate,rn2.eqs[1].rate)    # = false

I tried writing

rxs_1 = [Reaction(mm(X2,v1,K1), nothing, [X1], nothing, [1]),
       Reaction(mm(X3,v2,K2), nothing, [X2], nothing, [1]),
       Reaction(k1, [X1], [X3], [1], [1]),
       Reaction(k2, [X3], [X1], [1], [1]),
       Reaction(k3, [X3,X2], [X4,X1], [1,1], [1,1]),
       Reaction(k4, [X4,X1], [X3,X2], [1,1], [1,1]),
       Reaction(d, [X1], nothing, [1], nothing),
       Reaction(d, [X2], nothing, [1], nothing),
       Reaction(d, [X3], nothing, [1], nothing),
       Reaction(d, [X4], nothing, [1], nothing)]

but this yields a ERROR: UndefVarError: mm not defined

I'm trying to remove this test and rerun.

isaacsas · 2021-04-20T13:03:06Z

Looks good to me! We can always add back those other names if someone complains by defining them too.

isaacsas · 2021-04-20T13:03:20Z

~~We should figure out what is going on in that test.~~ I see, that test just no longer makes sense since we are now registering the new functions. Maybe just expand it to test if the value of the generated ODE LHS is the same (using isapprox)?

TorkelE · 2021-04-20T13:42:18Z

I'm thinking about removing that test (since as you say, it maybe just doesn't make sense anymore).

In a different file we do test whenever two such systems evaluate to similar things:

### Fetch required packages and reaction networks ###
using DiffEqBase, Catalyst, Random, Test
using ModelingToolkit: get_states, get_ps

using StableRNGs
rng = StableRNG(12345)

### Tests various cutom made functions ###
new_hill(x, v, k, n) = v*x^n/(k^n+x^n)
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

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

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]
    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())
end

This mostly worked, but there were occasional small errors, especially for larger input (I think especially having the exponent 2 in hill functions can amplify small numerical differences).

I modified this file by:

Increasing the error margin to 10e-6
removing the test for the largest factor (1e3)
and now ti seems to run fine.

TorkelE · 2021-04-20T13:57:59Z

Ok, I switched out the network which which didn't make sense to another (so we no longer test if the expression mm(...) is equal to v*X...).

isaacsas · 2021-04-20T14:02:53Z

test/custom_functions.jl

-    @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())
-end
+    @test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 10e-6)


Why such a loose tolerance now?

Some test fails when I do 1e-8. I'm quite sure that it's not the functions that are written wrongly (they pas lots of these tests), but that it is something which happens with the computer when it calculates large numbers. Don't actually know much about how this works, but looking into creating a minimal example.

Ok, turns out by making each reaction depend on their own X and Y value, things become a bit better (somehow the errors compounded in these variables across all reactions). Still can only push it to 1e-10 reliably, but it is still much better.

isaacsas · 2021-04-20T15:51:47Z

Probably should bump the ModelingToolkit version. The error is likely due to needing to regenerate the Latexify tests...

TorkelE · 2021-04-21T10:11:57Z

Turned out Latexify used some of the previous functions to detect and insert hill functions etc. That's fixed now, and things should (hopefully) run fine now.

There's a still a problem when using equation form:

using Catalyst
using Latexify
rn = @reaction_network begin
    mm(X,v,K), 0 --> X
end v K
latexify(convert(ODESystem, rn))

which doesn't unfold the mm(...) in the equation. However, the no-unfolding was already a case, and is probably due to that code not being in Catalyst. I'll register this as an issue in Latexify once this is done just to keep track on it.

TorkelE · 2021-04-21T11:27:17Z

Ok, just running latexify in a notebook, things look fine to me. However, there seems to be some difference in the actual codes used for comparison in the tests, e.g.

r = @reaction_network begin
    hillR(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
    hill(X5,v2,K2,n2), ∅ → X2
    hill(X3,v3,K3,n3), ∅ → X3
    hillR(X1,v4,K4,n4), ∅ → X4
    hill(X2,v5,K5,n5), ∅ → X5
    (k1,k2), X2 ⟷ X1 + 2X4
    (k3,k4), X4 ⟷ X3
    (k5,k6), 3X5 + X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5

# Latexify.@generate_test latexify(r)
@test latexify(r) == replace(
raw"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[K1^{n1} v1^{2} \left( \mathrm{X4}\left( t \right) \right)^{n1} \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1} \right) \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1} \right)] X1}\\
\ce{ \varnothing &->[v2 \left( \mathrm{X5}\left( t \right) \right)^{n2} \mathrm{inv}\left( K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2} \right)] X2}\\
\ce{ \varnothing &->[v3 \left( \mathrm{X3}\left( t \right) \right)^{n3} \mathrm{inv}\left( K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3} \right)] X3}\\
\ce{ \varnothing &->[v4 K4^{n4} \mathrm{inv}\left( K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4} \right)] X4}\\
\ce{ \varnothing &->[v5 \left( \mathrm{X2}\left( t \right) \right)^{n5} \mathrm{inv}\left( K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5} \right)] X5}\\
\ce{ X2 &<=>[k1][k2] X1 + 2 X4}\\
\ce{ X4 &<=>[k3][k4] X3}\\
\ce{ 3 X5 + X1 &<=>[k5][k6] X2}\\
\ce{ X1 &->[d1] \varnothing}\\
\ce{ X2 &->[d2] \varnothing}\\
\ce{ X3 &->[d3] \varnothing}\\
\ce{ X4 &->[d4] \varnothing}\\
\ce{ X5 &->[d5] \varnothing}
\end{align}
", "\r\n"=>"\n")

errors. The actual text now produced by latexify(r) is

L"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[\frac{v1 \left( \mathrm{X4}\left( t \right) \right)^{n1}}{K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1}} \frac{v1 K1^{n1}}{K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1}}] X1}\\
\ce{ \varnothing &->[\frac{v2 \left( \mathrm{X5}\left( t \right) \right)^{n2}}{K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2}}] X2}\\
\ce{ \varnothing &->[\frac{v3 \left( \mathrm{X3}\left( t \right) \right)^{n3}}{K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3}}] X3}\\
\ce{ \varnothing &->[\frac{v4 K4^{n4}}{K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4}}] X4}\\
\ce{ \varnothing &->[\frac{v5 \left( \mathrm{X2}\left( t \right) \right)^{n5}}{K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5}}] X5}\\
\ce{ X2 &<=>[k1][k2] X1 + 2 X4}\\
\ce{ X4 &<=>[k3][k4] X3}\\
\ce{ 3 X5 + X1 &<=>[k5][k6] X2}\\
\ce{ X1 &->[d1] \varnothing}\\
\ce{ X2 &->[d2] \varnothing}\\
\ce{ X3 &->[d3] \varnothing}\\
\ce{ X4 &->[d4] \varnothing}\\
\ce{ X5 &->[d5] \varnothing}
\end{align}
"

(haven't gotten Latexify.@generate_test latexify(r) to work)

@korsbo would you have time to have a look at this sometime in the near future? I've been unable to display the old text in Jupyetr/LaTeX without errors, so not sure exactly what is going on. We could just replace the old text with the new one (but it might be useful to have someone actually understand what is happening). Alternatively we can comment out the test, or put this PR on hold until you can check that we don't introduce something bad.

korsbo · 2021-04-21T12:04:24Z

@TorkelE I don't know if there are other problems but what you write in the last comment seems to be fine but that you just need to generate new tests.

For Latexify.@generate_test to work you need to have either xsel or xclip installed. IIRC then you're on ubuntu. apt-get install xsel.

Once you have that the Latexify.@generate_test call should work. Note though that is does not produce any visible output. That output it redirected to your clipboard, so run the command and then hit Ctrl + v in your editor to paste the test.

TorkelE · 2021-04-21T12:18:28Z

That's great!

I regenerated the tests and just pasted them into the doc, they now pass locally on my computer.

TorkelE · 2021-04-23T09:25:00Z

Added some additional tests, hoping to avoid the large drop in coverage. After this it should be good to go.

isaacsas · 2021-04-26T19:01:52Z

Looks great, thanks!

initial_upload

f33b862

ChrisRackauckas approved these changes Apr 20, 2021

View reviewed changes

TorkelE mentioned this pull request Apr 20, 2021

Add Double Hill Activator/Repressor function #315

Closed

skip_problem_test

56ada92

TorkelE added 2 commits April 20, 2021 15:36

modify_tests

deeaf7a

add_hillAR_test

46ea8da

swap out network which didn't make sense

4784520

isaacsas reviewed Apr 20, 2021

View reviewed changes

reduce_test_limit

419a6ca

TorkelE added 2 commits April 20, 2021 17:54

requires_mtk_5_15

39882c5

enable latexify to read functions

8fa0d1d

see_if_it_works_wothout_latexify_tests

b3aae31

re-generate latexify tests

b7b78b2

additional_tests

1a84043

TorkelE merged commit 0c41fdd into master Apr 26, 2021

isaacsas deleted the fix_hill_via_symbolics branch April 26, 2021 19:01

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

TorkelE commented Apr 20, 2021 •

edited

Loading

TorkelE commented Apr 20, 2021

isaacsas commented Apr 20, 2021

isaacsas commented Apr 20, 2021 •

edited

Loading

TorkelE commented Apr 20, 2021

TorkelE commented Apr 20, 2021

isaacsas Apr 20, 2021

TorkelE Apr 20, 2021

TorkelE Apr 20, 2021

isaacsas commented Apr 20, 2021

TorkelE commented Apr 21, 2021

TorkelE commented Apr 21, 2021

korsbo commented Apr 21, 2021

TorkelE commented Apr 21, 2021

TorkelE commented Apr 23, 2021

isaacsas commented Apr 26, 2021

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

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

Conversation

TorkelE commented Apr 20, 2021 • edited Loading

TorkelE commented Apr 20, 2021

isaacsas commented Apr 20, 2021

isaacsas commented Apr 20, 2021 • edited Loading

TorkelE commented Apr 20, 2021

TorkelE commented Apr 20, 2021

isaacsas Apr 20, 2021

Choose a reason for hiding this comment

TorkelE Apr 20, 2021

Choose a reason for hiding this comment

TorkelE Apr 20, 2021

Choose a reason for hiding this comment

isaacsas commented Apr 20, 2021

TorkelE commented Apr 21, 2021

TorkelE commented Apr 21, 2021

korsbo commented Apr 21, 2021

TorkelE commented Apr 21, 2021

TorkelE commented Apr 23, 2021

isaacsas commented Apr 26, 2021

TorkelE commented Apr 20, 2021 •

edited

Loading

isaacsas commented Apr 20, 2021 •

edited

Loading