Generating reactions programmatically using Reaction Systems

[yewalenikhil65]

Tutorial for Implementation of (Smoluchowski coagulation equation) population balance equation- by programmatically generating the reactions using ReactionSystem
Based on discussion at #292

[isaacsas]

@yewalenikhil65 I'll go through tomorrow morning and give you some comments. Got busy with meetings and class stuff today.

[yewalenikhil65]

walenikhil65 I'll go through tomorrow morning and give you some comments. Got busy with meetings and class stuff today.

No problem ,Sir

[isaacsas]

@yewalenikhil65 I started making some edits, but then ran into a couple questions. See below.

[isaacsas]

You should be able to git pull to sync my changes back to your local git repository.

[yewalenikhil65]

Hey @isaacsas updated according to your suggestions..

[isaacsas]

@yewalenikhil65 made some more updates. Slowly getting through it. (This week is pretty busy.)

I do have one question for you, see below.

[isaacsas]

Actually never mind; I forgot that I had figured out what I was confused about!

[yewalenikhil65]

Corrected a little mistake I spotted in writing kv for constant kernel case.. I like your changes in wording too

[isaacsas]

OK, I updated it one more time. I'd be happy to merge it now. Some notes:

• You weren't actually using ModelingToolkit / Catalyst in your previous code. When you created the MassActionJump and passed it to JumpProblem you were actually just using DiffEqJump directly.
• Along those lines, the ReactionSystem model used 2*k[n] for reactions with involving reactants that had the same number of monomers in them, however, what you were plotting and matching to the theory was just using k[n] for such reactions! I switched the ReactionSystem to use k[n] too since that is what made it match the theoretical formulas you used. Please double check that this is all correct though! You can see the new figure I made in the assets folder.

Let me know after you've taken a look, and if you approve I'll merge.

[yewalenikhil65]

OK, I updated it one more time. I'd be happy to merge it now. Some notes:

* You weren't actually using ModelingToolkit / Catalyst in your previous code. When you created the `MassActionJump` and passed it to `JumpProblem` you were actually just using `DiffEqJump` directly.

Yes, you are right. But i guess pushing reactions into for loop uses ModelingToolkit .

* Along those lines, the `ReactionSystem` model used `2*k[n]` for reactions with involving reactants that had the same number of monomers in them, however, what you were plotting and matching to the theory was just using `k[n]` for such reactions! I switched the `ReactionSystem` to use `k[n]` too since that is what made it match the theoretical formulas you used. Please double check that this is all correct though! You can see the new figure I made in the assets folder.

Again you are right.I plotted for just k[n] when I wrote the tutorial and uploaded that figure...But i have double checked using 2*k[n] and k[n]also. In both cases, multiplicative kernel doesn't match with analytical ..I do not know why!

Let me know after you've taken a look, and if you approve I'll merge.

We can merge it now. Only bit I am unsure is how to reason no-match of multiplicative kernel case, so i skipped commenting on it.

[isaacsas]

Hmm, perhaps we should just use the additive kernel then and remove the other two?

[yewalenikhil65]

Hmm, perhaps we should just use the additive kernel then and remove the other two?

Additive and constant kernels..Both are matching with analytical profiles..
I am not comfortable with removing multiplicative kernel even though its not matching. What do you think about calculation of the bulk volume of the system with which we are using in calculation of rate-constants ? Did we miss anything there ?

[yewalenikhil65]

Hey @isaacsas ..just realised something
This non-matching behaviour for 2*k[n] that you got is because of this 58a703e
I believe if you do add stoich construction as I did previously ,it matches in additive and constant kernel even if you use 2*k[n] as i mentioned earlier. Following part i mean

rx = [];              # empty-reaction vector
reactant_stoich = Array{Array{Pair{Int64,Int64},1},1}(undef, nr);
net_stoich = Array{Array{Pair{Int64,Int64},1},1}(undef, nr);
@time for n = 1:nr
    if (vᵢ[n] == vⱼ[n])    # checking the reactants
        push!(rx, Reaction(2*k[n], [ X[vᵢ[n]] ] ,[ X[sum_vᵢvⱼ[n]] ] ,[2],[1]));
        reactant_stoich[n] = [vᵢ[n] => 2];                                    # this
        net_stoich[n] = [vᵢ[n] => -2, sum_vᵢvⱼ[n] => 1];                 # this
    else
        push!(rx, Reaction(k[n], [ X[vᵢ[n]] , X[vⱼ[n]] ] ,[ X[sum_vᵢvⱼ[n]] ],
                                [1, 1],[1]));
        reactant_stoich[n] = [vᵢ[n] => 1 , vⱼ[n] => 1];                          # this
        net_stoich[n] = [vᵢ[n] => -1 , vⱼ[n] => -1 , sum_vᵢvⱼ[n] => 1];   # this
    end
end
rs = ReactionSystem(rx, t, X, k);
jumpsys = convert(JumpSystem, rs; combinatoric_ratelaws = true);
dprob = DiscreteProblem(jumpsys, u₀map, tspan, pₘₐₚ; parallel = true);
alg = Direct();                     stepper = SSAStepper();
mass_act_jump = MassActionJump(kₛ ,reactant_stoich, net_stoich);
jprob = @btime JumpProblem(dprob, alg ,mass_act_jump ,save_positions=(false, false));

I am confused why it is such though. If you remove this stoich construction, you don't get the match for 2*k[n] as you pointed

[isaacsas]

The MassActionJump you are building is using k[n] though not 2*k[n] for those reactions, right? When you use the mass_act_jump none of the ReactionSystem code is even getting used so it is irrelevant what you set as the rate in the Reactions.

[isaacsas]

Unfortunately, I haven't read this exact solution paper / method, so can't comment on why the multiplicative case does not match it.

[yewalenikhil65]

The MassActionJump you are building is using k[n] though not 2*k[n] for those reactions, right? When you use the mass_act_jump none of the ReactionSystem code is even getting used so it is irrelevant what you set as the rate in the Reactions.

Ahh.. understood now.. In that case, 58a703e seems correct and usage of k[n] seems more correct. I just latexified reactions into ODEs and checked if it matches with equation in wiki page.Usage of k[n] is appropriate one as you have done.
I think we can leave multiplicative case for now then. and merge only additive and constant kernels case ?

[isaacsas]

OK, can you remove the multiplicative case from the tutorial text and code? I'll then take a look and merge.

[yewalenikhil65]

OK, can you remove the multiplicative case from the tutorial text and code? I'll then take a look and merge.

@isaacsas ..Done!

[isaacsas]

All merged! Thanks!