Reaction Event Fluxes

[artur....@posteo.de]

I am currently having a look into reaction flux statistics.  I want to use Copasi to do the stochastic simulations for me.  I expected the “Direct Method” for the stochastic time course to simulate every single reaction event individually.

However, I do not understand, what exactly copasi exports while writing the report in the form of “Time and Reaction Event Fluxes”.  Chosing one second as a time interval I get non-integer results written into the report.  Does this really mean, there was a non-integer number of reaction events during that time interval of one second, while keeping integer values for the particle numbers?  Does Copasi do internal averaging or estimation of the stochastic fluxes, before writing the report?  

Is there any more detailed information/documentation than the user manual?  The latter does not really speak about reaction event fluxes, anyhow.

[juergen]

Hi Artur,

Copasi's stochastic simulation methods (Direct Method, Gibson + Bruck) internally simulate individual reaction events on discrete particle numbers as you would expect.

Currently, there is no way to output reaction propensities (e.g. from a stochastic perspective on the model) or the number of stochastic firings of each reaction from the graphical user interface directly (e.g. without programming using the C++ interface or other language bindings). What is meant by reaction event flux in the GUI is the rate of reaction events per time unit (e.g. from a deterministic/ODE perspective on the model). That is why you see fractional numbers in your report.

If you just want to follow the number of firings of each reaction, one easy hack would be to add one pseudo-product to each of the reactions. It would start from 0 and would then be incremented by 1 at each firing. These pseudo-species, in other words, would represent the accumulated number of firings of each reaction.

Anyway, we are currently working on a different internal data model that will most probably also allow direct access to the propensities and possibly the number of stochastic firings in the future.

Hope that helps!
Juergen

[artur....@posteo.de]

Single reaction events would be great!  Is there up to date documentation of the python bindings?  I wanted to integrate copasi to my automated command line workflow, anyhow.  However, I am a total beginner with using anything like this :-(

Somehow I am not sure I fully understand what you mean by

rate of reaction events per time unit (e.g. from a deterministic/ODE perspective on the model)

It might be that I am interested exactly in the differences of the fluxes as seen from stochastic or from determintistic models.  So I would prefer not to use these numbers – as I do not understand what is being done internally. 

[Stefan Hoops]

Hello Artur,

The fluxes are reaction events per time independent from the
integration method you choose. Thus non integer values are common even
for stochastic simulations if your kinetic parameters are non integers.
The difference between deterministic and stochastic integration is that
the changes in the fluxes are continuous vs step functions.

COPASI does not keep count of the number of actually realized reaction
events. However you can use Juergens trick to keep count of them by
adding a counting species as product. Comparing deterministic and
stochastic values will give you the difference between accumulated
fluxes.

Thanks,
Stefan

On Tue, 14 Apr 2015 09:24:24 -0700 (PDT)
artur....@posteo.de wrote:

> It might be that I am interested exactly in the differences of the
> fluxes as seen from stochastic or from determintistic models.  So I
> would prefer not to use these numbers – as I do not understand what
> is being done internally. 

--
Stefan Hoops, Ph.D.
Senior Project Associate
Virginia Bioinformatics Institute
Virginia Tech
1015 Life Science Circle (0477)
Blacksburg, Va 24061, USA

Phone: (540) 231-1799
Fax: (540) 231-2606
Email: sho...@vbi.vt.edu

[juergen]

Hi Artur,

just to clarify this a bit more: The flux of reaction events per time in a deterministic interpretation of the model is of course closely related to but not always the same as the average flux of reaction events per time in a stochastic interpretation of the model.

One example is the reaction 2 * A -> B with concentration [A] = 0.34 nmol in a volume of 5e-18 l, where the average stochastic flux is exactly zero and the deterministic flux is positive. This distinction can be very important for extinction events. Internally, Copasi handles this correctly during simulations.

However, as Stefan has already said, what Copasi WRITES OUT as "Reaction Event Fluxes" is always the reaction events per time in a deterministic interpretation of the model, no matter what simulation method the user has selected (e.g. no matter what interpretation of the model, deterministic or stochastic, Copasi internally uses for a simulation at that time).

Best wishes,
Juergen

[artur....@posteo.de]

All in all this sounds like I should not even try to look at histograms of the reaction event fluxes – or at least not interpret them as being the "true" distributions of the stochastic fluxes.

[Stefan Hoops]

Hello Artur,

I have not seen you model and therefore do not know what reaction
kinetics you use. Especially I do not know whether the
"Apply Stochastic Correction to Rate Laws"
flag is set on the model. If it is set the stochastic and deterministic
flux differ for small particle numbers in higher order kinetics as in
Juergens example. If it not set they are the same in all cases.

Thanks,
Stefan

[juergen]

Hi Artur,

well, I guess it depends on what you try to do. It might be the case that reaction event fluxes are still useful, or there might be other, less obvious ways of achieving what you need with Copasi.

Anyway, if this becomes a bit too technical or if you don't want to discuss your current research in public we can also take this offline. If you want, just contact me directly, and I would be happy to talk with you about your case a bit more in detail.

Best,
Juergen

[Pedro Mendes]

I am finding this discussion interesting and important too. So if you
want to continue in public, it won't be lost (plus will can be found
through search engines)

Pedro

> --
> You received this message because you are subscribed to the Google
> Groups "COPASI User Forum" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to copasi-user-fo...@googlegroups.com
> <mailto:copasi-user-fo...@googlegroups.com>.
> To post to this group, send email to copasi-u...@googlegroups.com
> <mailto:copasi-u...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/copasi-user-forum/bf021bd2-71d3-444e-88b3-52734a529651%40googlegroups.com
> <https://groups.google.com/d/msgid/copasi-user-forum/bf021bd2-71d3-444e-88b3-52734a529651%40googlegroups.com?utm_medium=email&utm_source=footer>.
> For more options, visit https://groups.google.com/d/optout.

--
Pedro Mendes
Professor of Computational Systems Biology
School of Computer Science
Manchester Centre for Integrative Systems Biology
University of Manchester

Manchester Institute of Biotechnology
131 Princess Street
Manchester, M1 7DN, U.K.

[Tolomè Pettinati]

Hi, I'm working with Artur. What we are interested in are the net number a reactions occurred along a randomly generated trajectory, and I agree that the easy hack that was proposed above should make this number accessible without affecting the reaction kinetics.

That's what we refer to as the reaction flux, and its definition does not invoke any underlying deterministic model. So, I made the following hypothesis. Consider reaction X + Y -> Z with mass-action law. What COPASI might be calculating is the product of the number of particles X*Y (which are random variables) multiplied by the rate constant of the reaction, which is a noninteger number. This is very different from actually counting the number of times reaction X + Y -> Z has happened.

Is that it?

[Juergen Pahle]

Hi Tolomè,

> Hi, I'm working with Artur. What we are interested in are the net number
> a reactions occurred along a randomly generated trajectory, and I agree
> that the easy hack that was proposed above should make this number
> accessible without affecting the reaction kinetics.

OK, good to hear that this hack is actually exactly what you need.

> That's what we refer to as the reaction flux, and its definition does
> not invoke any underlying deterministic model. So, I made the following
> hypothesis. Consider reaction X + Y -> Z with mass-action law. What
> COPASI might be calculating is the product of the number of particles
> X*Y (which are random variables) multiplied by the rate constant of the
> reaction, which is a noninteger number. This is very different from
> actually counting the number of times reaction X + Y -> Z has happened.
>
> Is that it?

If you do a stochastic simulation in Copasi it uses a Gillespie-type
algorithm internally that simulates each and every single, discrete
reactive event (with exponentially distributed time steps and reaction
probabilities proportional to their propensities etc.)

However, if you request the output of a "reaction event flux", e.g. if
you define a corresponding plot or have it in one of your reports, what
Copasi will calculate is exactly what you described. It is the -possibly
fractional- number of reaction events per time unit in a deterministic
interpretation of the model, which for the type of reaction you used in
your example is exactly the same as the average number of reaction
events per time unit in a stochastic interpretation.

Please keep in mind though:
- Copasi calculates these values only at the points in time where an
output is needed.
- These calculations don't really affect the stochastic simulation, e.g.
they are only for outputting purposes.
- Copasi actually only uses discrete particle numbers during stochastic
simulations. If you let Copasi output the reaction event fluxes during a
deterministic simulation then it will internally use fractional particle
numbers/concentrations for the calculation of reaction event fluxes.

Hope that helps!
Juergen

[matteoeo]

I think this definitively solves our problem. We are going to make a few tests on some simple models to check all of the above facts.

As a side comment, it would be extremely interesting to implement in COPASI the possibility of counting the reaction events as a further possible output, not reconstructing fluxes from the deterministic model on top of the stochastic populations of species (this should also not affect the way the stochastic simulation is conducted).

The theoretical reason behind this is that fluxes and populations obey different statistics and have different "large deviations" (statistical behavior in the long-time limit). This is extremely important for the community working on stochastic thermodynamics (rather than stochastic kinetimatics), since thermodynamic free-energy increases are defined along reaction events and cannot be reconstructed from knowledge of the populations.

Thanks,

tolomé

--
You received this message because you are subscribed to a topic in the Google Groups "COPASI User Forum" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/copasi-user-forum/fI2Fqn30Kn8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to copasi-user-forum+unsubscribe@googlegroups.com.
To post to this group, send email to copasi-user-forum@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/553012E4.8000308%40pahle.de.

[Pedro Mendes]

Dear Tolomè,

thanks very much for your suggestion, we will indeed investigate how we
can implement this. I agree with your view that this would be important.

Hopefully the suggested "hack" of adding a dummy product to each
reaction can work for now.

For others following this thread, in order to count the number of actual
reaction events that happened for some reaction R1, such as
A + B -> C
is to rewrite the reaction with an additional product, (eg "r1count"):
A + B -> C + r1count
then make sure this species "r1count" has initial particle number of 0,
and make sure it is not used in any other reaction. That way, r1count
will hold the number of reaction events that have happened since time 0
until the present (it any time point). This is not very elegant, but
should work fine.

Pedro

> https://groups.google.com/d/__topic/copasi-user-forum/__fI2Fqn30Kn8/unsubscribe
> <https://groups.google.com/d/topic/copasi-user-forum/fI2Fqn30Kn8/unsubscribe>.

> To unsubscribe from this group and all its topics, send an email to

> copasi-user-forum+unsubscribe@__googlegroups.com
> <mailto:copasi-user-forum%2Bunsu...@googlegroups.com>.

> To post to this group, send email to

> copasi-user-forum@__googlegroups.com
> <mailto:copasi-u...@googlegroups.com>.

> To view this discussion on the web visit

> https://groups.google.com/d/__msgid/copasi-user-forum/__553012E4.8000308%40pahle.de
> <https://groups.google.com/d/msgid/copasi-user-forum/553012E4.8000308%40pahle.de>.
>
> For more options, visit https://groups.google.com/d/__optout
> <https://groups.google.com/d/optout>.
>
>
> --
> You received this message because you are subscribed to the Google

> Groups "COPASI User Forum" group.

> To unsubscribe from this group and stop receiving emails from it, send
> an email to copasi-user-fo...@googlegroups.com
> <mailto:copasi-user-fo...@googlegroups.com>.

> To post to this group, send email to copasi-u...@googlegroups.com
> <mailto:copasi-u...@googlegroups.com>.

> To view this discussion on the web visit

> https://groups.google.com/d/msgid/copasi-user-forum/CANirP0K-0pO2K0WXSC9ju9uG%3DkNis0kMvDjMVhczk_PHpMDiTQ%40mail.gmail.com
> <https://groups.google.com/d/msgid/copasi-user-forum/CANirP0K-0pO2K0WXSC9ju9uG%3DkNis0kMvDjMVhczk_PHpMDiTQ%40mail.gmail.com?utm_medium=email&utm_source=footer>.

> For more options, visit https://groups.google.com/d/optout.

[juergen]

Hi,

as Pedro already said, we are always grateful for suggestions and constructive input. Many thanks for that!

It is likely that there are a number of interesting features at the interface between kinetic models and thermodynamics that could be implemented in Copasi, as this seems to be a very active area of research currently.

Let's keep in touch,
Juergen

[Dante Kienigiel]

Hey Juergen!

I'm wondering if something like this was ever implemented to avoid the "hack"? I have the "more general" issue of wanting to know which is the event that happened at a given time (from which it would then be easy to count). The alternative to what is proposed with the dummy species here is to do a .diff() over the species at each timestep, see which changed and infer which was the reaction which made the changes. This shifts compute from the simulation to post-processing, but is also very inefficient.

Best,

Dante

[juergen]

Hi Dante,
Unfortunately, we haven't implemented it yet. I will bring this up in the next COPASI developer meeting, to discuss how something like this could be implemented and whether we have the capacity atm actually to do it. Thank you for the reminder though!
Best wishes,
Juergen

[matteoeo]

Hi, I' m not sure I belong to this conversation any more, given the great time lapse. Anyway I am still interested in the issue.

Best

Matteo Polettini 

--
You received this message because you are subscribed to a topic in the Google Groups "COPASI User Forum" group.

To unsubscribe from this topic, visit https://groups.google.com/d/topic/copasi-user-forum/fI2Fqn30Kn8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to copasi-user-fo...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/copasi-user-forum/961bfa0c-77a5-4660-a718-863a9a78c94bn%40googlegroups.com.