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

[Rainer M. <ra...@tb...>]

Andrew

>
> I think its too early to do this.
> We need to have a rock solid understanding of the math first.

i agree. this can all get quite tricky.

>
> S1 in C  /* declared with conctration units */
> S1->junk : k[S1]
>       /* this rate law is a function of species concentrations
> returning in a rate of change of substance i.e. k has units litre per
> sec */

then k is not the original constant for a unimolecular reaction 
[ 1 / sec ] but has already been corrected for some specific volume!

i.e. it is a volume-dependent parameter!!!

i think you won't find such type of rate constants outside of the SBML 
specific kinetic law constructs!

and it actually shouldn't be done this way but the kinetic law should 
explicitly contain the conversion, i.e.:

k*V*S1

the reason is not only that a) that k in litre/sec is quite an ugly 
construct but also that b) this model would simply become wrong if someone 
changes the size of the compartment, i.e. uses a volume different from the 
one once used for correction of the forward rate k!!

this is also touching a very old discussion on rate `constants' in 
general. mathematicians and physicists always rolled their eyes when 
chemists called them `constants' because they are actually functions of 
temperature, pH  etc.

using volume-dependent parameters would take this imprecision to a really 
ugly level!

i think this is quite dangerous and there should actually be big 
warnings in the kinetic law and/or compartment sections of the SBML specs, 
that people should not convert rate constants but explicitly include the 
volume in the kinetic law.



i will answer on the rest next week :)

Rainer



On Fri, 16 Jun 2006, Andrew Finney wrote:



> Rainer
>
> you wrote:
>
>> 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
>>
>
> I think its too early to do this.
> We need to have a rock solid understanding of the math first.
>
> For example:
>
> given the following SBML system:
>
> S1 in C  /* declared with conctration units */
> S1->junk : k[S1]
> 	/* this rate law is a function of species concentrations
> returning in a rate of change of substance i.e. k has units litre per
> sec */
>
> dC/dt = 1
>
> The current niave forumlation is:
>
> [S1] = [S1(0)]
> d[S1]/dt = k[S1]/C
> dC/dt = 1
>
> Your proposed solution I think (which should reserve the units of k) is:
>
> S1 = [S1(0)] * C
> dS1/dt = k.(S1/C)
> 	/* species amounts have to be converted back to concetrations
> for the constant*/
> dC/dt = 1
>
> I think someone who knows needs to verify that this is correct.
> I don't believe it is but I'm not qualified to comment.
>
>>
>> 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?
>
> event assignments are going to be fine in the existing scheme
> but the rest are I agree problematic.
>
>>
>>
>> *) 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?
>>
>
> I don't think there is a need to worry:
>
> The code for constructing the minimal sets of assignment rules
> is not aware of the type of variables (compartment, species or
> parameter)
> it is dealing with.  This
> code is called after the initial set of odes and a 'maximal' set of
> assignment rules have been constructed.
>
> If you get the ODEs and 'maximal' assignment rule set right then the
> code should
> do the right thing.  If there is a bug in the minimization code
> resist the urge to hack it so that it treats symbols differently
> depending on their type.
>
>> 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?
>>
>
> hmm, don't understand what you mean
>
>>
>> *) 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 ?
>>
>
> I think its best if species are output in the units in which they
> were defined on the sbml species element i.e. by default concentration
> units
>
>>
>>
>> *) 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.
>>
>
> Yuck
>
>>
>> Any suggestions, opinions or comments ???
>>
>
> I feel out of my depth!
>
>>
>> Please add other problems here:
>>
>>
>
> Yes, we need to check that the SBML semantic test suite has got the
> right results
> I can't remember if variable compartments are in the current version
>
> yours Andrew
>
>> -----Original Message-----
>> From: sbm...@li...
>> [mailto:sbm...@li...] On
>> Behalf Of Rainer Machne
>> Sent: 16 June 2006 16:56
>> To: SOSlib Development
>> Subject: Re: [SOSlib-devel] variable compartment bug :
>> correction[Scanned]
>>
>> 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
>>>
>>
>>
>> _______________________________________________
>> 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
>