feasibility questions

[Jan Meyer]

Dear developers,

I am currently thinking about how to use kappa/Kasim for modelling certain polymerizations that I am interested in. For this purpose a few questions concerning feasibility of my ideas arised. I would be very thankful if you could comment on them:

Exothermic reactions:

In the reactions that I investigate the produced heat of reaction can be quite significant and therefore lead to increase in temperature. This of course affects the reaction rates again (assuming f.e. an Arrhenius like behavior). Is there  any option to include the heat build up leading to an increase in temperature. Am I for instance able to track which reactions have been performed and update the kinetic constants over the course of the reaction?

 

Chain length dependent reactions:

For polymerization processes it is common to scale the reaction rates according to the length (or mol weight) of the chains reacting. F.e let’s assume we have only one species of Monomers “x-A-x” having two equal reactive sites”x”. So A can react with itself to form chains. I would like to take into account that

x-A-x + x-A-x -> x-A-y-A-x might have a different reaction constant than the corresponding reaction of two dimers: x-A-y-A-x + x-A-y-A-x -> x-A-y-A -y-A -y-A -x 

Apart from comments to my questions any advice from you is very welcome!

Many thanks in advance.

Best regards,

Jan Meyer

[Walter Fontana]

You could use reaction-specific tokens  (Manual section 2.4.5). 

'r1'  A(x[.]), A([y[.]) ->  A(x[1]), A([y[1])  |  + r1_token @ 'rate_constant'

You would define 'rate_constant' as a function of the heat generated per reaction event and of r1_token, which gives you the number of times the reaction has fired.  

Check the manual to understand how tokens work, and how to declare and initialize them.

Careful about rules with molecular ambiguity (section 3.7); e.g. 'r1' above will also form rings, not just linear polymers. Depending on your intentions, you may need to associate two rate constants with ‘r1’ (and set the unimolecular rate constant to zero to prevent ring formation); see the Manual.

wf

--
You received this message because you are subscribed to the Google Groups "kappa-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to kappa-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/2114958f-986b-4338-872b-6f57e5614f0b%40googlegroups.com.

[Jan Meyer]

Thank you very much for your response. I will definetly have a look at the manual.

Do you think that the token functionality could also be a solution to the second point I raised? My first approach towards this problem would be the usage of counters. Giving every monomer a counter that is raised according to the size of the oligomer that the monomer might be part of. Than I would use this counter value to rescale the kinetics. I am also a bit worried about the overhead this might produce. Do you maybe have any experience with that?

Thanks in advance.

Have a nice weekend!

Best regards,

Jan

To unsubscribe from this group and stop receiving emails from it, send an email to kappa...@googlegroups.com.

[Walter Fontana]

Your second point (length-dependent rate constants): In principle, counters are an option, but, as you say, you would incur substantial overhead. I’m not sure, however, what the latest implementation of counters looks like. Pierre?

This said, if you need such a fine-grained length-dependence that the fusion between any two lengths has a specific rate constant, you seem no longer to deal with a world characterized by rules but with a world described by reactions. In that case you might be better off writing a dedicated stochastic process model yourself or using some available Gillespie simulator. If you just want a quick exploration, you could “abuse” rule-based approaches by programmatically generating reaction rules (i.e. fully refined rules that completely specify molecular species) for each length combination. Even if you have 10,000 rules (100*100), you should be fine in principle, but I don’t know off the top of my head whether Pierre Boutillier and Jean Krivine’s elegant “domain structure” built upon initialization (in KaSim) is able to handle this in reasonable time (RAM issues) or whether it can be switched off. Pierre? (I know Pierre is busy right now, but maybe he can get back to the group in a week.)

I’m not up to date on BNGL; maybe BNGL has a solution that is efficient. Anyone here?

On the theory side: If you are aiming for an equilibrium distribution, the abundance of each polymer length will only depend on its energy content, not on how it has been assembled. In your case of homopolymers, each bond would have an affinity (an exponential of free energy) that depended on how it was formed. I have difficulties understanding how that is possible at the level of the agent abstraction used in rule-based models. There would have to be contrived configurational aspects to retain memory of how each bond was formed, in which case you are well outside the abstraction level of geometry-less Kappa. An exception is the formation of rings: the ring closure reaction will depend on the length of the ring, meaning that one bond in the ring contributes a different free energy — but it doesn't matter which one it is because of “rotational” symmetry. So the only “memory” retained in the free energy comes from the simple fact that the polymer is a ring (of a given size), which is perfectly within the agent abstraction of Kappa (or BNGL). Happy to take this offline, if needed.

Best,

w

To unsubscribe from this group and stop receiving emails from it, send an email to kappa-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/8dd3325c-32c4-448a-8669-3a20d0c26ce0%40googlegroups.com.

[John Bachman]

Walter alluded to BNGL--with respect to changing the rate as a function of polymer length, you may be able to use local functions in BNG/NFSim. I have not used them myself, but the NFSim paper describes an example of using a local function that is evaluated over a complex to determine the rate of a reaction rule:

From Sneddon et al. 2011 (https://www.ncbi.nlm.nih.gov/pubmed/21186362):

Experimental and theoretical studies have shown that the activ- ity of each signaling team can be represented in a coarse-grained manner using a Monod-Wyman-Changeux model12,33. In NFsim this is achieved by defining a local function, ‘pOn’, that depends on the methylation and ligand binding state of each receptor in a given signaling team (Fig. 5b). The argument ‘x’ refers to the particular complex over which pOn must be evaluated. The pOn function can then be used to represent how the rate of methyla- tion of one receptor depends on the activity of its signaling team. This is illustrated in the block of code in Figure 5b, where the rate of the methylation reaction is multiplied by pOn(x). The ‘%x’ prefixed to ‘Rec(m~?)’ tells NFsim to evaluate the function pOn(x) on the entire complex (signaling team) that is connected to the indicated receptor. The local-function syntax also allows the arrangement of receptors in a signaling team to be modified with- out changing rule definitions.

and a more recent review of BNGL (https://academic.oup.com/bioinformatics/article/32/21/3366/2415298)

Local functions enable a single rule to specify many reactions with rates that depend on the specific properties of the reacting species

Hope that is helpful!

John

 

To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/BD533FAA-7279-4CD1-AFEB-129C19EC41E7%40gmail.com.

[Bill Hlavacek]

Hi Jan,

As Walter mentioned...

There's a chance - possibly small - that BioNetGen/NFsim (https://github.com/RuleWorld/bionetgen) and BioNetGen language (BNGL) could be helpful.

You could use a "global function" to define a rate law incorporating temperature dependence. But you'd need some sort of model for how temperature increases as your exothermic reactions occur. I'm not aware of an example of this type of thing in BNGL-formatted models I've seen.

You could use a "local function" to query properties of reactants matching a pattern in a rule to adjust the rate of reaction based on reactant properties. I'm reasonably confident that you could get the desired length-dependent rate law. However, functions, IMHO, are a very poorly documented feature of BNGL. Others might disagree.

A quick overview is presented here:

https://drive.google.com/file/d/0B2lPm2_GUE01X3ZaamZxUl80NTA/edit

A better resource might be Justin Hogg's PhD thesis from Univeristy of Pittsburgh. But I'm not sure about that.

If you want to give this a try, it could be helpful to consult someone familiar with BNGL

Another suggestion... you might want to look at Kappa/KaSim- and BioNetGen-like software tools developed specifically for non-biological chemical systems, like NetGen, RMG (https://rmg.mit.edu), etc. There's list of tools given in this recent review paper from Linda Broadbelt:

https://aiche.onlinelibrary.wiley.com/doi/10.1002/aic.16663 

But... I suspect that these tools can't handle polymerization reactions, except by truncation of network generation.

Good luck!

--Bill

--

You received this message because you are subscribed to the Google Groups "kappa-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to kappa-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/2114958f-986b-4338-872b-6f57e5614f0b%40googlegroups.com.

[Bill Hlavacek]

I'm not sure I would recommend looking at the NFsim paper for guidance on functions. This paper *did* introduce functions, but the way functions work has changed since then.

To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/CAAs77FDZoJUO5mMg_%3DP0GfCAEKHrOyZKSMWM8LO%3Dko7ZLJWcMw%40mail.gmail.com.

[Jan Meyer]

Dear all,
Thank you so much for all your responses and the expertise you shared with me. I will have a deeper look at the suggested papers in the upcoming week.

Maybe a few words to my situation to clarify certain things. I work in the R&D Department of a chemical company. I deal with molecular modeling of polymeric systems on different length scales. However, one of my responsibilities is an in-house kmc code. The Problem with this code is, that it has been developed decades ago and it is not easy to access for modifications. Therefore I was just curious if there might be more state of the art Software packages that I could use. But of course I wanted to make sure that I wont get bottlenecked right from the start.

While I did not really find to many solutions in the chemical community the Bio-chemical community did offer much more :).

For the Chain length dependance we usually asume something like 1/n to 1/n^2 as a scaling factor where n is the number of monomers.

For the ring building I wanted to use md-simulations to predict for which length of a chain it is reasonable to form a ring.

Again, I thank you all for your Input. Please feel free to give me any other advice!

Have a nice Weekend.
Best regards from Germany,
Jan

[Pierre Boutillier]

Hi all,

Extremely late reply indeed, Sorry.

- Dealing with heat is indeed doable with tokens. Despite fighting with the syntax and having to be careful with counter intuitive behaviors of tokens (For example, nothing prevent them to become negative and a negative activity for a rule is a weird thing :-)), there should be no difficulty. Complains as soon as you don't find/understand something.

- Enumerating complexes up to length 100 (even 1000) is piece of cake for the simulator but I'm not sure it is of any help if the question is really to deal with arbitrary complexes.

Technical details for the insiders and the interested: Problems with the sharing structure arise when you have 1 agent type with 17 phosphorylation sites and you enumerate all the possibilities because by default we are enumerating the 17! possible paths to discover all these possibilities (because sites can be visited in any order...) To avoid that, you can indeed launch the simulator with option `-sharing None` or wait that I finally take the time to debug the real solution to the problem which is stuck at 85% completion in the github branches `tons_of_rules` for 2 months... The case of discovering a linear polymer is not problematic because at each step of your exploration, you only have 2 possibilities : discovering a new agent because the only site you haven't inspected yet is bound to someone or completing the discovery of a complex because this site is free! So you end up with only a few hundreds (thousands) elements.

- There is no way for a Kappa agent to introspect the size of the complex it is part of. I don't see any way to encode this knowledge in Kappa.

Technical details for the insiders and the interested: Yes, at the price of changing any bits of the simulator, we could add this knowledge with a computational complexity as bad as it is already when you deal with rule involving "molecular ambiguity". The real price to pay is in the increase in complexity of the engine. First, although we discuss that already for other use, there is no syntax to give a name to an agent in a rule yet. That means there are some work to do in order to refer in the rule rate to "(the size of the complex in which appears) the second agent of type EGF mentioned in the left hand side". Moreover, all instance of a rule has the same activity for now so the data-structure used to store and pick at random a instance is uniform. If now the rate depends on the instance, we must introduce the notion of weight in between occurrences. Again nothing is impossible here, but it reprensent a substantial surgery :-)

Best,

Pierre B.

To view this discussion on the web visit https://groups.google.com/d/msgid/kappa-users/BD533FAA-7279-4CD1-AFEB-129C19EC41E7%40gmail.com.