Re: [SOSlib-devel] variable compartment bug : correction

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

hi All

the easiest way to solve the variable compartment problem without 
having a second data->values array might probably be to:

a) leave away the division by volume in Species_odeFromReactions
    to get ODEs in substance/time instead of concentration/time and convert
    all initial values to substance instead of concentration values

b) use data->values for substance instead of concentration values

c) calculate concentrations directly in evaluateAST:
    the easiest way would be to further extend our AST indexing with an
    boolean field isConcentration,
    and a second index for the data->values field which contains the
    the value of compartment by which this substance value must be divided

However, there might be several problems:

*) as already pointed out by Andrew: we will have to 
"...  block the combination of rate rules for species in concentration 
units with rate rules for compartment volumes"

i guess, the same is true for variable compartments in general, i.e. 
compartments changed by assignment, rate rules and by event assignments

right?

*) One problem that I can think of now is the current 
`observables' implementation, i.e. that assignment rules are only 
calculated if the value is required for other formulae.

The volume might however not be explicitly present in formulae and thus 
the `observables' code might not detect that a potential assignment rule 
for the volume must always be calculated. Andrew, what do you think?

on the other hand, one could argue that the reaction system is formulated 
wrong, because a variable compartment MUST explicitly appear in kinetic 
laws. It is wrong to include the compartment correction into a parameter, 
if the compartment is variable, right? So in a correct 
variable-compartment model an assignment rule for a variable compartment 
will always be in the list of assignment rules that have to be evaluated.

Or are there cases where this might not be true?

*) The next problem is the output of results, e.g. should cvodeResults and 
sbmlResults contain the substance or the concentration values?
Should we have a function like 
IntegratorInstance_getVariableValueAsConcentration ?

*) We could make it even more complicated by allowing two ODE versions,
i.e. calculate the ODEs in substance/time only if variable compartments 
appear in the model.

Any suggestions, opinions or comments ???

Please add other problems here:

Rainer

On Wed, 14 Jun 2006, Rainer Machne wrote:

> Stefan
>
>
>> dX/dt = fX() / V0 * V
>> that's all!
>> of course, via fX the concs enter the odes for substances.
>
>
> That means we would divide the ODE by the initial volume and multiply
> again by the current volume?
>
> Why that? Why can't be just leave the volumes away?
>
> How to account for different volumes in one ODE, which happens for
> multiple compartment and species involved in transport reactions?
>
> e.g. if
> compartment V1 contains A; in substance units, with conc. A/V1 = [A] = a
> compartment V2 contains B and C
>
> A->B, would be the transport reaction with a kinetic law:  kt*a*V1
>
> B->C, is an irreversible conversion with kinetic law: kf*b*V2
>
> then our concentration ODE would be
>
> db/dt = (kt*a*V1 - kf*b*V2) / V2
>
> right?
>
> Can't we just leave away the last / V2 to have the correct substance ODE?
>
> dB/dt = kt*a*V1 - kf*b*V2
>
> where V1 and V2 can be variable?
>
> Rainer
>
>
> On Wed, 14 Jun 2006, Stefan Mueller wrote:
>
>> andrew, lukas, rainer!
>>
>> sorry guys, the problem of compartments is so difficult :)
>> that i made an error, too.
>>
>> given the law for concentrations dx/dt = k*y*z for a fixed volume V,
>> the "law" for substances X=x*V is correct even for a time-dependent volume:
>> dX/dt = (k*V)*y*z = (k/V)*Y*Z
>> whereas the law for concentrations has to be corrected:
>> dx/dt = k*y*z - x (dV/dt) / V
>>
>> if we write down odes for substances instead of concs, we can calculate V
>> either by an ode or by an assignment rule.
>> if we used odes for concs this would be more complicated.
>>
>> so with the sbml definition of a kinetic law for a fixed volume V0
>> dX/dt_V0 = V0*dx/dt = V0*k*x*y = fX(concs)
>> the odes for substances can be written as:
>> dX/dt = fX() / V0 * V
>> that's all!
>> of course, via fX the concs enter the odes for substances.
>>
>> so. i hope that was the last error.
>> anyway, let's discuss it on friday!
>>
>> cheers,
>> stefan.
>>
>>
>> Am Dienstag 13 Juni 2006 19:01 schrieb Stefan Mueller:
>>> andrew, lukas, rainer!
>>>
>>> the problem of compartments (multiple or even variable) is not that easy.
>>>
>>> the definition of kinetic laws in sbml is conceptionally bad.
>>> ok. the sbml people wanted to ease the treatment of multiple compartments.
>>> and it's working...
>>> but among other things, it complicates the treatment of variable
>>> compartments. anyway, we can "work around" the current definition. (see
>>> below.)
>>>
>>> first, how should it be:
>>> in chemistry, rate laws are defined for concs,
>>> e.g. dx/dt=k*x*y.
>>> and there is a reason for this: the law is independent of the volume V!
>>> if you formulate "laws" for substances X=x*V,
>>> you get dX/dt=(k*V)*x*y=(k/V)*X*Y.
>>> this is not a law of chemistry, but depends also on the volume V.
>>> moreover, for time-dependent V the above "law" is wrong!
>>>
>>> having this in mind (?), sbml people specified dX/dt for a given volume V0,
>>> let's call it dX/dt_V0=(k*V0)*x*y.
>>> by simple calculus, we can derive dX/dt from X=x*V:
>>> ----------------------------------------------------
>>> dX/dt = dX/dt_V0 / V0 * V + X * dV/dt / V
>>> ----------------------------------------------------
>>> dV/dt can be given explicitly or derived (!) from an assignment rule.
>>> as rainer already suspected, we have to use odes for substances to allow
>>> for variable compartments. (although the kinetic laws still depend on
>>> concs).
>>>
>>> i hope i clarified the situation a little bit.
>>>
>>> cheers,
>>> stefan.
>>>
>>> Am Dienstag 13 Juni 2006 15:25 schrieb Andrew Finney:
>>>> Rainer
>>>>
>>>> Slow down!!
>>>>
>>>> I don't actually understand the argument:
>>>> why can't you have the ODEs as ODEs of concentrations:
>>>>
>>>> if S1 is in a compartment with volume C then
>>>>
>>>> d[S1]/dt = (f1() + f2() + ... fn()) / C
>>>>
>>>> f1..fn are the reaction rates which are
>>>> expressions that have units substance/time
>>>>
>>>> f1...fn are functions of the concentrations of the
>>>> species [S1]..[Sn]
>>>>
>>>> I thought that's what we'd implemented.
>>>>
>>>> Why do you need two sets of values?
>>>>
>>>> I would strongly urge we go for the approach is having
>>>> one set of values and carefully converting the units as is appropriate.
>>>>
>>>> fundamentally what's the problem?
>>>>
>>>> In the SBML you sent the only fly in the ointment
>>>> is that the initial value for the species is
>>>> expressed in substance units via the initialAmount attribute rather than
>>>> the initialConcentration
>>>> attribute which is in concentration units.
>>>>
>>>> All you need to do is ensure that the initialAmount
>>>> values are used to calculate the initialConcentration values again by
>>>> dividing by the volume.  Have
>>>> we implemented that?
>>>>
>>>> Note the initial volume of the compartment is not specified in the xml.
>>>>
>>>> yours Andrew
>>>>
>>>>> -----Original Message-----
>>>>> From: sbm...@li...
>>>>> [mailto:sbm...@li...] On
>>>>> Behalf Of Rainer Machne
>>>>> Sent: 13 June 2006 11:42
>>>>> To: SOSlib Development
>>>>> Subject: [SOSlib-devel] variable compartment bug[Scanned]
>>>>>
>>>>> Folks
>>>>>
>>>>> Lukas has identified another serious bug!
>>>>>
>>>>> SOSlib doesn't support variable compartments.
>>>>>
>>>>> We initialize sundials with ODEs for concentrations d[A]/dt,
>>>>> and sundials calculates [A] + (d[A]/dt)*h, where h is a small
>>>>> time step and [A] is the initial condition or the current value.
>>>>>
>>>>> If we have a variable, e.g time-dependent compartment, then
>>>>> d[A]/dt will change with changing compartment because the
>>>>> compartment size is part of the right hand side.
>>>>>
>>>>> However, not only d[A]/dt is dependent on the volume, but
>>>>> also [A], i.e.
>>>>> the current value in data->value!
>>>>>
>>>>> The attached example file exemplifies this.
>>>>>
>>>>> So we actually should use an ODE system in amount/time
>>>>> instead of concentration/time, and have a second data->value
>>>>> array for concentrations to be used in formula evaluation.
>>>>>
>>>>> Didn't we discuss this once? I remember that I had started to
>>>>> do this (ODEs in amount/time), but then it got more
>>>>> complicated then I thought and something else came up.
>>>>>
>>>>> Can you think of any problems we will run into? Andrew, how
>>>>> about the compiling version?
>>>>>
>>>>> Rainer
>>>>
>>>> _______________________________________________
>>>> sbmlsolver-devel mailing list
>>>> sbm...@li...
>>>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>>
>>
>> _______________________________________________
>> sbmlsolver-devel mailing list
>> sbm...@li...
>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>>
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>