Looking for advice: Copasi, steady state, eigenvalues, Jacobian

[Thomas Ligon]

I am looking for some advice in using Copasi for calculating steady state and examining the results, including eigenvalues and the Jacobian.

In Copasi, steady state succeeds if I set Resolution to 1e-5, but the results are not completely convincing, because some species show a (final) concentration different from what I see in the time course. I also have some eigenvalues of 0.

-> How can I find out which reactions are responsible for the zero eigenvalues?

As far as I can tell, the Jacobian looks OK (all diagonal values are negative).

-> Is there more information that I am overlooking?

When I run a time course, it looks like everything has settled into a steady state by a time of 2e+5. However, if I try to run the time course to 1e20, it fails at about 2e+19, even though the plot and the results look OK. In comparison, when I run this in AMICI, the time course fails at about 3e+199 and steady state fails completely.

Copasi error running time course to 1e20:

>EXCEPTION 2021-05-24T16:12:18<

  CTrajectoryMethod (6): Deterministic integration failed. LSODA reported:

  DLSODA-  At T (=R1) and step size H (=R2), the    

        corrector convergence failed repeatedly     

        or with ABS(H) = HMIN   

  In above message, R1 = '2.06262e+19', R2 = '7.92357e+08'

  Please see result for indications of numerical instability.

AMICI error tying to run steady state:

[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : At t = 3.18384e+199, the right-hand side routine failed in an unrecoverable manner. 

[Warning] AMICI:simulation: AMICI simulation failed:

Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.

Error occurred in:

stacktrace not available on windows platforms

[Frank Bergmann]

I think it would help to get some more information: 

* what is the acceptance criterium for the steady state (we added a simple 'rate' criterium in recent versions, that looks at the rate of changes for accepting steady states. This can be combined with the distance criterium and with both active, we hardly see false positives anymore). 

* it would be interesting to see the result in the protocol tab

* instead of looking at the time course concentration results and the particle number plots, i think it would be helpful to also look at the rates of changes plot / report values.  (in the output assistant they are created with ("Concentration Rates, Volume Rates, and Global Quantity Rates" / "Time, Concentration Rates, Volume Rates, and Global Quantity Rates")

best

Frank

[sven....@gmail.com]

Hello Thomas,

while Steady States are a simple concept, they are surprisingly difficult from a numerical point of view since there are several issues that can go wrong independently.

The most problematic point is often (as Frank mentioned) the criterium that COPASI uses to decide if a Steady state was found. On the other hand the error message wrt the integration is difficult to interprete without knowing the model.

When you say that you have Eigenvalues of 0, is that in the reduced or the full jacobian? Some 0 Eigenvalues are the result of the structure of the model (rather than the specific kinetics) and COPASI should be able to deal with them, in this case the 0 Eigenvalues should not be present in the "reduced" Jacobian.

In any case, if you would send the model (or a version of it that shows the problem) to me or Frank (we would obviously treat it confidentially) we could help much better.

best,

Sven

[Mendes,Pedro]

Thomas,

it is legitimate for dynamical systems to have zero eigenvalues; when this happens the system is described as "marginally stable". Perturbations to that steady state do not grow (ie it is not unstable) but also do not decay back to the original steady state, but rather end up in a nearby steady state. The most well known system with this property is the Lotka-Volterra model.

This is not to say that those zero eigenvalues are necessarily correct, as Sven noted, you get zero eigenvalues in the full Jacobian matrix if you have mass conservation relations. However in the reduced Jacobian (ie the one where mass conservation was removed) if there is a zero, it is either a very small number that appears to be zero or a true zero (there is not much difference between these, a tiny eigenvalue is a system that takes a very very long time to decay back to the original steady state, so effectively almost the same as marginal stability)

The one thing that I find interesting in your post is that you say that COPASI finds the steady state when you set the tolerance to 1e-5. Does it not find it if you have it set to 1e-6 ? In that case it may be not a true steady state, but just a very very slow transient state. But the criterion used (distance, rate, or both) is also important here.

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.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/7dc456e3-643c-490a-ab93-7db54588dde7n%40googlegroups.com.

-- 
Pedro Mendes, PhD
Professor and Director,
Richard D. Berlin Center for Cell Analysis and Modeling
University of Connecticut School of Medicine
group website: http://www.comp-sys-bio.org

[Upi Bhalla]

Dear Pedro, Frank et al,

I'm also running into difficulties with the steady-state calculation in this simulation. It fails and the result section is full of nans and infs. I just used the default settings. The same model runs fine for time-series, but also fails when I try to do model reduction by time scale separation. I'm not really sure where to start with debugging this. I'm running COPASI 4.30 build 240 on Ubuntu 18.04.

Thank you,

Upi

[Mendes,Pedro]

Dear Upi,

Steady state is surprisingly difficult to properly define numerically. Since version 4.28 (build 226) we have added a the option for the user to select among three different criteria for steady state convergence: "distance", "rate" and "distance and rate".  The default criterion is the most strict one "distance and rate" in order to avoid false positives. However this sometimes misses states that could be considered steady states. So I would advise you to try changing the criterion and try again. Of course the "Resolution" is also important and the smallest this value is the hardest it is to find a steady state.

COPASI's manual page includes a description of how it calculates the steady state: http://copasi.org/Support/User_Manual/Methods/Steady_State_Calculation/
However, unfortunately, we have not yet added the convergence criteria to COPASI's manual, but it should come shortly. These are:

Criterion "Rate": the largest of all the right-hand side of the ODEs needs to be smaller than the Resolution specified (in concentration units)

Criterion "Distance": the Newton method increment needs to be smaller than the Resolution specified
 
Criterion "Distance and rate": both of the above need to be satisfied

By the way, there is a new version available with a few bug fixes (4.31).

Pedro

To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/cdf8abdf-f552-4dc1-99b2-74e1c878a9a9n%40googlegroups.com.

[Thomas Ligon]

Dear Pedro,

thanks again for the help. Now I have a steady state, but there is still an integration error.

When I run Steady State and set Target Criterion to Rate, it runs to completion quickly and shows very plausible results. In addition, all eigenvalues are negative (no zeros). When I look at Result / Species, I can see that the Transition Time starts at about 6e-6, with many values around 1e+5 and a maximum of 6e+7 (except for two!). This seems plausible. However, two species have Transition Times of 6e+19 and 3.7e+20. Let's look at just this last one:

This species (PTEN) is a protein, used as a modifier in a separate reaction, which shouldn't affect this time. It is produced by a constant production, and degraded by a standard degradation reaction. If I set the initial concentration to zero, the time course looks exactly like what I expect, beginning with a linear growth, then leveling off. After about 2e5 seconds, it reaches almost a maximum of about 6.5e5 and continues to grow asymptotically.

When I run a Time Course to a time of 1e+20, it terminates with an integration error at 1e+18 or 1e+19:

>EXCEPTION 2021-06-02T09:07:59<

  CTrajectoryMethod (6): Deterministic integration failed. LSODA reported:

  DLSODA-  At T (=R1) and step size H (=R2), the    

        corrector convergence failed repeatedly     

        or with ABS(H) = HMIN   

  In above message, R1 = '7.89558e+18', R2 = '1.15365e+09'

  Please see result for indications of numerical instability.

When I run the time course in AMICI, it has an integration error at around 1e+199.

These integration errors are the problem that I am currently trying to solve, because they interfere with PETab/pyPESTO/AMICI.

In summary, the integration error is a problem, but the error doesn't tell me anything about which species or reaction in causing it. The very large transition time reported by Copasi Steady State looks like a very useful hint. Now, I don't understand why this species should have such a high transition time, and I don't see if or why this species is causing an integration error. Looking at the plots of all concentrations, all rates, or PTEN alone hasn't uncovered the problem yet.

[sven....@gmail.com]

Dear all,

having looked at the model with simulations of moderate dureation (e.g. 1e5), it seems the model cannot be in a steady state because degraded_protein is produced but never consumed. It seems to grow forever. The rate of production is indeed smalll (in the order of 1e-9) so that you can force COPASI to "find" a "steady state" by relaxing the resolution of the steady state criterium. Actually all the other species have a much smaller rate of change than "degraded_protein" (by many orders of magnitude; so I think the model is in a steady state except this one species. If I set the degraded_protein to "fixed" I can find a steady state with reasonable resolution.

As for the simulation, I cannot just now reproduce the problem since my computer is quite slow and it would need to run a long time to reach near 1e20s. However, I would not be surprised if the integration eventually fails if you simulate a model where most of the variables reach a steady state quickly but one variable keeps rising infinitely for an astronomical time scale. At some point the ratio between the growing variable and the rest would overwhelm the computers floating point accuracy (or confuse the tolerance adjustment of the integrator).

In any case, for "debugging" or analsing a situation like this I would usually simulate for a reasonable long (but not astronomical) time, and then look the the actual rates of the concentration changes, rather than the transition times. This allows you to find out which variables might be the problem for finding a steady state. If simulation for extremely long times is necessary to settle the model in steady state, it indicates that something with the time scales of the model is off.

Hope this helps,

Sven

[Thomas Ligon]

Dear Sven,

thanks! Your help is greatly appreciated. However, I can't find the "degraded_protein" you are referring to. I searched in Copasi and in SBML.

In fact, this is a problem that I had in the past, and it was normally very visible in the concentration plot, clearly showing a species that grows in an unbounded fashion to + or - infinity.

For this, there is a small problem that I had but fixed. For a constant production or a degradation, in Copasi, I can just leave the "empty set" species out, but in Cell Designer, I need to include it, so I have one such species for each compartment. Here is one example in Cell Designer:

DEGRADED s978 nil mc inside Amount 0.0 count true true 0 true

In SBML, this is:

<species metaid="s978" id="s978" name="nil" compartment="mc" initialAmount="0" substanceUnits="count" hasOnlySubstanceUnits="true" boundaryCondition="true" charge="0" constant="true">

Setting constant="true" should prevent the problem.

By the way, in order to run a time course to an astronomical time, I just set Intervals to 100, and then it runs fast.

[sven....@gmail.com]

Dear Thomas,

it seems I got the models mixed up. I was referring to the model Upi posted, where there is a species called "degraded_protein".

Sven

[Upi Bhalla]

Dear Sven,

Thanks for the diagnosis of the problem with steady-states in my model, indeed you are right. I was also able to get a steady-state solution once I fixed the degraded_protein term and put in a small relaxation of the criterion.

Best,

Upi

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/I1B5uKi2LiQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to copasi-user-fo...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/06715c85-a23f-47b3-8494-990ab43c73c9n%40googlegroups.com.