networks to add as correctness tests

[isaacsas]

This issue is to collate references that have good examples we can add as correctness tests:

Variable rate tests:

• https://pubs.aip.org/aip/jcp/article/143/5/054104/194056/Simulation-of-biochemical-reactions-with-time look good.

Hybrid model tests:

• jump-SDE (https://doi.org/10.1063/1.1835951)
• jump-SDE and jump-ODE (https://doi.org/10.1063/1.1505860)
• jump-ODE with non-negativity focus (https://doi.org/10.1186/s12859-019-2836-z)

[isaacsas]

For these examples we may want to setup a system where we simulate them in an alternative external tool we trust, and then the test compares against saved solutions from that tool. (Since in many cases they will not have analytic solution values we can compare against.)

[isaacsas]

This is the variable rate DNA example from the RSSA paper above:

    rn = @reaction_network gene_model begin
        @parameters begin
            r[1:10] = [0.043, 0.0007, 0.0715, 0.0039, 0.0199, 0.4791, 0.00019, 0.8765, 0.083, 0.5]
            k = -log(2) / 30
        end
        r[1], RNA --> RNA + M
        r[2], M --> ∅
        r[3], DNAD --> RNA + DNAD
        r[4], RNA --> 0
        r[5] * exp(k*t), DNA + D --> DNAD
        r[6], DNAD --> DNA + D
        r[7] * exp(k*t), DNAD + D --> DNA2D  
        r[8], DNA2D --> DNA + D
        r[9] * exp(k*t), 2M --> D
        r[10], D --> 2M
    end
    u0 = [:DNA => 10, :M => 10, :D => 30, :RNA => 0, :DNAD => 0, :DNA2D => 0]
    jinputs = JumpInputs(rn, u0, (0.0, 120.0); save_positions = (false, false))
    jprob = JumpProblem(jinputs; rng, save_positions = (false, false))
    Nsims = 1600
    tsave = range(0.0, 120.0, length = 121)
    Mmean = zeros(length(tsave))
    Dmean = zeros(length(tsave))
    for n in 1:Nsims
        sol = solve(jprob, Tsit5(); saveat = 1.0)
        Mmean .+= sol[:M]
        Dmean .+= sol[:D]
    end
    Mmean ./= Nsims
    Dmean ./= Nsims

It might make a nice benchmark example, and should be a good example for testing variable rate aggregator correctness (but we'd need to simulate it with another tool to then check consistency in a test as I mention above).

[ChrisRackauckas]

https://docs.sciml.ai/SciMLBenchmarksOutput/stable/HybridJumps/Synapse/ should get updated with the new variable rate jumps too

[isaacsas]

May need Romain's help for that given the PDMP errors. But we could expand it to test using the new aggregators too. I think it should actually be possible to just write the whole model in Catalyst now (aside from the bounds that are needed for CoEvolve).

[ChrisRackauckas]

I just bumped on dev and it looked fine. So it's now setup for someone to do the next part.

I think it should actually be possible to just write the whole model in Catalyst now (aside from the bounds that are needed for CoEvolve).

That would be good to do anyways.

[isaacsas]

However, I have a suspicion that getting the new method to outperform our old approach could be very dependent on both the ODE integrator being used and how the quadrature is handled within the callback.