Negative concentrations

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
Tuoma
Posts: 21
Joined: Mon Apr 23, 2018 6:59 am

Negative concentrations

Post by Tuoma »

Is there a way to force the concentrations to be non-negative? I'm using CVode and adaptive time step in a system of >100 reactions, and very often end up with negative concentrations for some of the species, although initial concentrations are all non-negative. Consequently (at least so it seems), the simulation often terminates with an error "CVode-- At t = XXX and h = 1e-10, the error test failed repeatedly or with |h| = hmin." I've tried different values of tolerance, hmin and hmax: sometimes choosing the right hmax and/or tolerance avoids the problem but in most cases doesn't.
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: Negative concentrations

Post by ramcdougal »

Not directly. If you're not using mass-action, it's possible to write reactions that mathematically demand concentrations go negative. Can you tell us more about your reaction specification?

Assuming some are not mass action, have you checked to make sure the reaction rates are not such that the concentration "should" go negative given the mathematics?

For a trivial example, suppose I have an rxd.Rate that clears a molecule at a constant non-zero rate. If that molecule is never produced, then mathematically, the concentration should go negative.

e.g. Don't do

Code: Select all

clear_ip3 = rxd.Rate(ip3, -1)
Instead make clearance depend on the concentration (this is biologically plausible: an enzyme can't break down a molecule if it can't find it):

Code: Select all

clear_ip3 = rxd.Rate(ip3, -ip3)
Obviously the rate could be governed by other kinetics (Hill-type, etc...)
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: Negative concentrations

Post by ramcdougal »

Also: do the concentrations go negative with fixed-step integration?
Tuoma
Posts: 21
Joined: Mon Apr 23, 2018 6:59 am

Re: Negative concentrations

Post by Tuoma »

Thanks for the reply. The system was a mass-action reaction system, and the rates of all injected species were positive. But I finally managed to solve the problem by making the tolerance way smaller (down to 3e-9) than it used to be. The downside is of course that now my simulations are much slower, but I guess there's no silver bullet for that problem. The fixed time step method was even slower and gave negative concentrations as well.
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Negative concentrations

Post by ted »

In a system with so many reactions it seems likely that at least a few will involve state variables that are much smaller than others. For such cases, it is often useful to use state-specific error tolerance scale multipliers. NEURON's GUI simplifies this task tremendously.

1. Make sure you have imported gui. After model setup is complete, in the NEURON Main Menu toolbar click on
Tools/VariableStepControl
2. In the VariableTimeStep panel, verify that the "Use variable dt" box shows a checkmark.
3. Click on the "Atol Scale Tool" button.
4. The bottom panel in the Absolute Tolerance Scale Factors panel shows the state names in the first column and the error tolerance scale factors in the second column (default value 1).
Click on the "Analysis Run" button and see what happens to the values in the second column.
After the end of the simulation, click on the "Rescale" button.
5. Now see if you can relax the value of the absolute error tolerance (see numeric field next to the Absolute Tolerance button in the VariableTimeStep panel), i.e. increase it from 1e-9 to something more reasonable, without getting negative concentrations.

If that helps, you can save the configured VariableTimeStep tool to its own .ses file for future re-use.
1. Click on
NEURON Main Menu/Window/Print & File Window Manager
then examine the "screen map" (red rectangle on the left side of the Print & File Window Manager) to identify which of the numbered blue rectangles corresponds to the VariableTimeStep panel on your desktop. To confirm that you have found the correct one, use your mouse to drag the VariableTimeStep panel to a new position, and release the mouse button. The blue outline that moved is the one you want.
2. Select that rectangle by clicking inside it, and you should see it and its number appear in the PFWM's "page layout" area (red rectangle on the right side of the PFWM).
3. Click on "Session/Save selected" in the PFWM and a new tool will pop up that asks you to tell it what name you want the .ses file to have. Something like varstep.ses would be rather mnemonic. Then click on the Save to file button.

After exiting NEURON you can insert the following statement into your own program, after model setup code is complete and after h.cvode_active(1) (the statement you should be using to activate CVODE), but before executing h.run() or whatever you're doing to launch simulations.
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: Negative concentrations

Post by ramcdougal »

I second Ted's point about the need for state-specific error tolerance scale multipliers.

Unfortunately, the rxd state variables are not integrated into the VariableStepControl GUI's analysis at this time.

You can, however, specify a species-specific atolscale. The rxd.Species constructor prototype is:

Code: Select all

 rxd.Species(regions, d=0, name=None, charge=0, initial=None, atolscale=1)
As an example of how to do this, consider calcium. A "typical" calcium concentration is around 100 nM or 1e-4 mM (recall: NEURON's concentrations are measured in mM). With the default absolute tolerance settings of 1e-3 and atolscale=1, you'd expect zero digits of precision in calcium concentration. That's not good and could easily lead to negative concentrations. If, on the other hand, you set calcium's atolscale=1e-4 (the "typical" value), you'd expect 3 digits of precision. If there are large swings across orders of magnitude, you have to think about the tradeoff between performance and precision, but this should give you some hints.
Post Reply