sbmlsolver-devel — developer's email list for the SBML ODE Solver Library

Re: [SOSlib-devel] implementing variable compartments[Scanned]

[Andrew F. <af...@ph...>]

Folks

So what your suggesting is pretty simple:
for any SBML system always transform the system to
create ODEs of rates of substance change
and assignment rules that derive concentrations from substance.

That's pretty straightforward I think.  I'm happy.

Note:

I don't think you need to actually substitute the x.amount/volume into
the
ODEs if you have the assignment rules.  In SOSlib

x2 =3D x2.amount/volume
dx1.amount/dt =3D k * x2
dx2.amount/dt =3D ???

should be equivalent to:

dx1.amount/dt =3D k * (x2.amount/volume)
dx2.amount/dt =3D ???

This works because the relevant assignment rules are evaluated in the
CVODE rhs function
before evaluating the ODEs themselves.

Or am I missing something? An issue with the jacobian???

If the above isn't true then the way assignment rules work
in the general case is wrong too.

yours Andrew

> -----Original Message-----
> From: sbm...@li...=20
> [mailto:sbm...@li...] On=20
> Behalf Of Stefan Mueller
> Sent: 21 June 2006 16:51
> To: sbm...@li...
> Subject: Re: [SOSlib-devel] implementing variable=20
> compartments[Scanned]
>=20
> rainer
>=20
> > Stefan, thanks for the clarification in your last email.
> >
> > i think i got it now :)
> >
> > you wrote:
> > > the ode for conc can be derived by differentiating x=3DX/V:
> > > dx/dt =3D (dX/dt)/V - x*(dV/dt)/V
> > > where dX/dt is given by the kinetic law above.
> >
> > dX/dt is just the sum of the +/- (stoichiometry*kineticLaw)=20
> > expressions for all reactions where X is a=20
> product/reactant, i.e. it=20
> > is a function of concentrations x, y, z ...
> >
> > Right? (If not forget, the rest ...)
>=20
> right.
>=20
> > So actually, the first part of this equation (dX/dt)/V is=20
> the current=20
> > form of ODEs in SOSlib!
> >
> > Could the implementation of variable compartments in SOSlib=20
> would be=20
> > as easy as:
> >
> > 1) add the expression " - x*(dV/dt)/V "
> > to the concentration ODEs for species in a variable compartment?
>=20
> yes.
>=20
> > 2a) if the compartment is defined by a rate rule, this rate=20
> rule could=20
> > be appended to all ODEs directly, i.e.
> >
> > newODE =3D currentODE - (x/V) * rateRule
>=20
> yes.
>=20
> > 2b) if the compartment is defined by an assignment rule, we=20
> would have=20
> > differentiate it towards time? For this we need to convert the=20
> > assignment rule to a form, where differentiateAST recognizes time=20
> > dependent variables that are in the assignment rule, e.g.=20
> x(t), which=20
> > appears only as ASTNode "x". Is this possible/easy?
>=20
> possible.
> easy? differentiateAST has to extended.
> there must be differentiation with respect to time, but only=20
> variables for which there is an ode are differentiated.
> (and the rhs of the ode is the result.)
> the chain rule does the rest of the work...
>=20
> > 2c) if the compartment is changed by events, we don't need=20
> to append=20
> > anything but just recalculate concentrations and re-initialize the=20
> > solver, just as we do now for events.
> >
> >
> > If above is right - but i guess it might just not be that easy - we=20
> > would have extracted two alternatives:
> >
> > A) use ODEs in substance/time, write assignment rules for=20
> x.amount =3D x=20
> > / volume and replace all x by x.amount in ODEs, as Andrew suggested
>=20
> you mean: x =3D x.amount / volume !?
>=20
> > or
> >
> > B) use above approach.
> >
> > Is that it?
>=20
> yes.
> i vote for alternative A :)
>=20
> reason:
> (i) in alternative B you have to do point 1 and 2a, the=20
> volume correction, for each ode.
> (ii) in alternative B you have to do point 2b, the additional=20
> differentiation.
> (iii) the use of assignment rules x =3D x.amount / volume are=20
> not a criterion.
> they can be used in both alternatives.
>=20
> cheers, s
>=20
>=20
> > Rainer
> > _______________________________________________
> > sbmlsolver-devel mailing list
> > sbm...@li...
> > https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>=20
>=20
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>=20

Re: [SOSlib-devel] mail for SBML list: variable Compartments in SBML[Scanned]

[Andrew F. <af...@ph...>]

That's better=20

> -----Original Message-----
> From: sbm...@li...=20
> [mailto:sbm...@li...] On=20
> Behalf Of Stefan Mueller
> Sent: 21 June 2006 14:36
> To: sbm...@li...
> Subject: Re: [SOSlib-devel] mail for SBML list: variable=20
> Compartments in SBML[Scanned]
>=20
> rainer,
> there are more mistakes!
> (please see below)
>=20
> > hi,
> >
> > find below a draft for bringing this issue to the SBML list, as=20
> > suggested by Andrew.
> >
> > could you review? is the problem stated clearly and correctly?
> >
> >
> > Dear SBML Community,
> >
> > As we recently wanted to model variable compartments we=20
> realized that=20
> > our tool - SOSlib - actually does not solve models with variable=20
> > compartments correctly.
> >
> > We wonder how other tools approach this topic.
> >
> > In summary:
> >
> > SOSlib constructs ODEs as recommended by the SBML specs, e.g.:
> >
> > species X1 in compartment V1, with concentration x1 =3D X1/V1=20
> species X2=20
> > and X3 in V2, with conc. as above
> >
> > Reactions:
> > X1->X2; k1*X1*V1 	// a translocation, k1 [1/sec]
> > X2->X3; k2*X2*V2    // a unimolecular reaction, k2 [1/sec]
> >
> > Derived ODEs:
> > dx1/dt =3D - ( k1 * X1* V1 ) / V1
> > dx2/dt =3D  ( k1 * X1* V1 - k2 * X2 * V2 ) / V2
> >
> > the kinetic laws are the SBML specific conversion of `classic' rate=20
> > law to substance/time units, instead of concentrations/time=20
> to allow=20
> > for easy software handling of multiple compartment models.
> >
> > If however, the volumes in the model are variable with time, i.e.=20
> > either an assignment rule or a rate rule, these equations=20
> would not be correct.
>=20
> the following section is not very clear and/or wrong:
> (please don't send this to the sbml list!)
>=20
> > A numerical integration routine would start with the initial=20
> > concentration
> > x(0) and add the changes during time t:
> >
> > x(t) =3D x(0) +  dx/dt * t
> >
> > i.e.
> >
> > x(t) =3D X(0)/V(0) + dx/dt * t
> >
> > dx/dt is calculated correctly from above ODE because it=20
> also accounts=20
> > for changing compartments by explicit inclusion of the=20
> volume in the=20
> > kinetic law.
> >
> > However, the change of the initial concentration x(0) is=20
> not accounted=20
> > for anywhere. The solver should rather return:
> >
> > x(t) =3D X(0)/V(t) + dx/dt * t
> >
> > If x is a function of time and of volume V, x(t,V(t)), then the=20
> > correct ODE would actually be
> >
> > dx/dt =3D dx/dt + dx/dV * dV/dt
> >
> > Right?
>=20
> no!
> please use the following math:
> (i already wrote it down in one of my previous mails :) )
>=20
> definitions:
> x(t) ... conc
> X(t) ... substance
> V(t) ... compartment size
>=20
> given a kinetic law for a reaction,
> e.g. X->Y: r*x*V or X+Y->Z: r*x*y*V,
> which explicitly includes the compartment size (!), the=20
> corresponding ode for substance X amounts to:
> dX/dt =3D - r*x*V or - r*x*y*V,
> or replacing conc by substances:
> dX/dt =3D - r*X or - r*X*Y/V
>=20
> the ode for conc can be derived by differentiating x=3DX/V:
> dx/dt =3D (dX/dt)/V - x*(dV/dt)/V
> where dX/dt is given by the kinetic law above.
> i.e. we would get a correction term for variable compartment size.
> this is the most important reason for using odes for substances!
>=20
> > We are currently discussing different solutions and their=20
> dangers for=20
> > SOSlib, i.e. how to correctly convert the ODE system to an=20
> ODE system=20
> > in amount/time instead of concentration/time or whether the above=20
> > partial differential equation could or should be=20
> implemented instead.
> >
> > It would be helpful for us to hear how other tools handle variable=20
> > compartments.
> >
> >
> > And BTW, it might be good to indicate it as good practice NOT to=20
> > convert parameters e.g. from [1/sec] to [litre/sec] to=20
> obtain correct=20
> > units for SBML kinetic laws, but to explicitly keep the=20
> volume in the=20
> > formula! This would allow changing the compartment size or=20
> even making=20
> > the variable, at least for tools which do it correctly.
> > E.g. the model in the specs, section 5.6 will not be valid=20
> anymore if=20
> > someone wanted to change the volume sizes!
> >
> > Such rather unconvential volume-dependent rate `constants' can be=20
> > quite dangerous!
> >
> >
> > Rainer
> >
> >
> >
> > Rainer Machne
> > Institute for Theoretical Chemistry
> > University of Vienna
> > Austria
> >
> >
> > _______________________________________________
> > sbmlsolver-devel mailing list
> > sbm...@li...
> > https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>=20
>=20
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>=20

Re: [SOSlib-devel] Compilation problems of libsbml 2.3.4

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

hi Eryk

i haven't seen this problem yet and never had problems with libsbml 2.3.4 
so far.

if you can't solve the problem and no one else here has a suggestion, you 
could also ask at the libsbml discussion list, see 
http://www.sbml.org/forums/

Rainer

On Wed, 21 Jun 2006, W Eryk Wolski wrote:

> Hi,
>
>
> Compilation of libsbml 2.3.4 fails with
>
> g++ -I.. -DPACKAGE_VERSION=\"2.3.4\"   -I. -DLINUX  -g -O2  -fPIC -MT
> SBMLHandler.o -MD -MP -MF ".deps/SBMLHandler.Po" -c -o SBMLHandler.o
> SBMLHandler.cpp
> SBMLSchemaInputSource.h:79: error: extra qualification
> 'SBMLSchemaInputSource::' on member 'makeStream'
> make[3]: *** [SBMLHandler.o] Error 1
> make[3]: Leaving directory
> `/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src/sbml'
> make[2]: *** [all-this-dir] Error 2
> make[2]: Leaving directory
> `/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src/sbml'
> make[1]: *** [sbml-recurse] Error 2
> make[1]: Leaving directory `/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src'
> make: *** [src-recurse] Error 2
>
>
> gcc --version
> gcc (GCC) 4.1.0 (SUSE Linux)
> Copyright (C) 2006 Free Software Foundation, Inc.
> This is free software; see the source for copying conditions.  There is NO
> warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
>
>
> Suggestions?
>
> Eryk
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>

[SOSlib-devel] Compilation problems of libsbml 2.3.4

[W E. W. <wew...@gm...>]

Hi,


Compilation of libsbml 2.3.4 fails with

g++ -I.. -DPACKAGE_VERSION=\"2.3.4\"   -I. -DLINUX  -g -O2  -fPIC -MT
SBMLHandler.o -MD -MP -MF ".deps/SBMLHandler.Po" -c -o SBMLHandler.o
SBMLHandler.cpp
SBMLSchemaInputSource.h:79: error: extra qualification
'SBMLSchemaInputSource::' on member 'makeStream'
make[3]: *** [SBMLHandler.o] Error 1
make[3]: Leaving directory
`/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src/sbml'
make[2]: *** [all-this-dir] Error 2
make[2]: Leaving directory
`/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src/sbml'
make[1]: *** [sbml-recurse] Error 2
make[1]: Leaving directory `/home/witek/devel/sbmlodesolve/libsbml-2.3.4/src'
make: *** [src-recurse] Error 2


gcc --version
gcc (GCC) 4.1.0 (SUSE Linux)
Copyright (C) 2006 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


Suggestions?

Eryk

Re: [SOSlib-devel] implementing variable compartments

[Stefan M. <ste...@oe...>]

rainer

> Stefan, thanks for the clarification in your last email.
>
> i think i got it now :)
>
> you wrote:
> > the ode for conc can be derived by differentiating x=X/V:
> > dx/dt = (dX/dt)/V - x*(dV/dt)/V
> > where dX/dt is given by the kinetic law above.
>
> dX/dt is just the sum of the +/- (stoichiometry*kineticLaw) expressions
> for all reactions where X is a product/reactant, i.e. it is a function of
> concentrations x, y, z ...
>
> Right? (If not forget, the rest ...)

right.

> So actually, the first part of this equation (dX/dt)/V is the current form
> of ODEs in SOSlib!
>
> Could the implementation of variable compartments in SOSlib would be as
> easy as:
>
> 1) add the expression " - x*(dV/dt)/V "
> to the concentration ODEs for species in a variable compartment?

yes.

> 2a) if the compartment is defined by a rate rule, this rate rule could be
> appended to all ODEs directly, i.e.
>
> newODE = currentODE - (x/V) * rateRule

yes.

> 2b) if the compartment is defined by an assignment rule, we would have
> differentiate it towards time? For this we need to convert the assignment
> rule to a form, where differentiateAST recognizes time dependent variables
> that are in the assignment rule, e.g. x(t), which appears only as ASTNode
> "x". Is this possible/easy?

possible.
easy? differentiateAST has to extended.
there must be differentiation with respect to time,
but only variables for which there is an ode are differentiated.
(and the rhs of the ode is the result.)
the chain rule does the rest of the work...

> 2c) if the compartment is changed by events, we don't need to append
> anything but just recalculate concentrations and re-initialize the solver,
> just as we do now for events.
>
>
> If above is right - but i guess it might just not be that easy -
> we would have extracted two alternatives:
>
> A) use ODEs in substance/time, write assignment rules for
> x.amount = x / volume and replace all x by x.amount in ODEs, as Andrew
> suggested

you mean: x = x.amount / volume !?

> or
>
> B) use above approach.
>
> Is that it?

yes.
i vote for alternative A :)

reason:
(i) in alternative B you have to do point 1 and 2a, the volume correction,
for each ode.
(ii) in alternative B you have to do point 2b, the additional differentiation.
(iii) the use of assignment rules x = x.amount / volume are not a criterion.
they can be used in both alternatives.

cheers, s


> Rainer
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel

[SOSlib-devel] implementing variable compartments

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

Stefan, thanks for the clarification in your last email.

i think i got it now :)

you wrote:

> the ode for conc can be derived by differentiating x=X/V:
> dx/dt = (dX/dt)/V - x*(dV/dt)/V
> where dX/dt is given by the kinetic law above.


dX/dt is just the sum of the +/- (stoichiometry*kineticLaw) expressions 
for all reactions where X is a product/reactant, i.e. it is a function of 
concentrations x, y, z ...

Right? (If not forget, the rest ...)


So actually, the first part of this equation (dX/dt)/V is the current form 
of ODEs in SOSlib!


Could the implementation of variable compartments in SOSlib would be as 
easy as:

1) add the expression " - x*(dV/dt)/V "
to the concentration ODEs for species in a variable compartment?


2a) if the compartment is defined by a rate rule, this rate rule could be 
appended to all ODEs directly, i.e.

newODE = currentODE - (x/V) * rateRule


2b) if the compartment is defined by an assignment rule, we would have
differentiate it towards time? For this we need to convert the assignment 
rule to a form, where differentiateAST recognizes time dependent variables 
that are in the assignment rule, e.g. x(t), which appears only as ASTNode 
"x". Is this possible/easy?

2c) if the compartment is changed by events, we don't need to append 
anything but just recalculate concentrations and re-initialize the solver, 
just as we do now for events.


If above is right - but i guess it might just not be that easy - 
we would have extracted two alternatives:

A) use ODEs in substance/time, write assignment rules for
x.amount = x / volume and replace all x by x.amount in ODEs, as Andrew 
suggested

or

B) use above approach.

Is that it?

Rainer

Re: [SOSlib-devel] mail for SBML list: variable Compartments in SBML

[Stefan M. <ste...@oe...>]

rainer,
there are more mistakes!
(please see below)

> hi,
>
> find below a draft for bringing this issue to the SBML list, as suggested
> by Andrew.
>
> could you review? is the problem stated clearly and correctly?
>
>
> Dear SBML Community,
>
> As we recently wanted to model variable compartments we realized that our
> tool - SOSlib - actually does not solve models with variable compartments
> correctly.
>
> We wonder how other tools approach this topic.
>
> In summary:
>
> SOSlib constructs ODEs as recommended by the SBML specs, e.g.:
>
> species X1 in compartment V1, with concentration x1 = X1/V1
> species X2 and X3 in V2, with conc. as above
>
> Reactions:
> X1->X2; k1*X1*V1 	// a translocation, k1 [1/sec]
> X2->X3; k2*X2*V2    // a unimolecular reaction, k2 [1/sec]
>
> Derived ODEs:
> dx1/dt = - ( k1 * X1* V1 ) / V1
> dx2/dt =  ( k1 * X1* V1 - k2 * X2 * V2 ) / V2
>
> the kinetic laws are the SBML specific conversion of `classic' rate law
> to substance/time units, instead of concentrations/time to allow for
> easy software handling of multiple compartment models.
>
> If however, the volumes in the model are variable with time, i.e. either
> an assignment rule or a rate rule, these equations would not be correct.

the following section is not very clear and/or wrong:
(please don't send this to the sbml list!)

> A numerical integration routine would start with the initial concentration
> x(0) and add the changes during time t:
>
> x(t) = x(0) +  dx/dt * t
>
> i.e.
>
> x(t) = X(0)/V(0) + dx/dt * t
>
> dx/dt is calculated correctly from above ODE because it also accounts for
> changing compartments by explicit inclusion of the volume in the kinetic
> law.
>
> However, the change of the initial concentration x(0) is not accounted for
> anywhere. The solver should rather return:
>
> x(t) = X(0)/V(t) + dx/dt * t
>
> If x is a function of time and of volume V, x(t,V(t)),
> then the correct ODE would actually be
>
> dx/dt = dx/dt + dx/dV * dV/dt
>
> Right?

no!
please use the following math:
(i already wrote it down in one of my previous mails :) )

definitions:
x(t) ... conc
X(t) ... substance
V(t) ... compartment size

given a kinetic law for a reaction,
e.g. X->Y: r*x*V or X+Y->Z: r*x*y*V,
which explicitly includes the compartment size (!),
the corresponding ode for substance X amounts to:
dX/dt = - r*x*V or - r*x*y*V,
or replacing conc by substances:
dX/dt = - r*X or - r*X*Y/V

the ode for conc can be derived by differentiating x=X/V:
dx/dt = (dX/dt)/V - x*(dV/dt)/V
where dX/dt is given by the kinetic law above.
i.e. we would get a correction term for variable compartment size.
this is the most important reason for using odes for substances!

> We are currently discussing different solutions and their dangers
> for SOSlib, i.e. how to correctly convert the ODE system to an ODE system
> in amount/time instead of concentration/time or whether the above partial
> differential equation could or should be implemented instead.
>
> It would be helpful for us to hear how other tools handle variable
> compartments.
>
>
> And BTW, it might be good to indicate it as good practice NOT to convert
> parameters e.g. from [1/sec] to [litre/sec] to obtain correct units for
> SBML kinetic laws, but to explicitly keep the volume in the formula! This
> would allow changing the compartment size or even making the variable, at
> least for tools which do it correctly.
> E.g. the model in the specs, section 5.6 will not be valid anymore if
> someone wanted to change the volume sizes!
>
> Such rather unconvential volume-dependent rate `constants' can be
> quite dangerous!
>
>
> Rainer
>
>
>
> Rainer Machne
> Institute for Theoretical Chemistry
> University of Vienna
> Austria
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel

Re: [SOSlib-devel] mail for SBML list: variable Compartments in SBML

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

I wrote:

> Reactions:
> X1->X2; k1*X1*V1      // a translocation, k1 [1/sec]
> X2->X3; k2*X2*V2    // a unimolecular reaction, k2 [1/sec]
>
> Derived ODEs:
> dx1/dt = - ( k1 * X1* V1 ) / V1
> dx2/dt =  ( k1 * X1* V1 - k2 * X2 * V2 ) / V2

should of course be (conc. x1 instead of subst. X1 in math):

> Reactions:
> X1->X2; k1*x1*V1      // a translocation, k1 [1/sec]
> X2->X3; k2*x2*V2    // a unimolecular reaction, k2 [1/sec]
>
> Derived ODEs:
> dx1/dt = - ( k1 * x1* V1 ) / V1
> dx2/dt =  ( k1 * x1* V1 - k2 * x2 * V2 ) / V2


the example could of course also be replaced by a simpler model with only 
one reaction in one variable compartment to make the point.

r


> hi,
>
> find below a draft for bringing this issue to the SBML list, as suggested
> by Andrew.
>
> could you review? is the problem stated clearly and correctly?
>
>
> Dear SBML Community,
>
> As we recently wanted to model variable compartments we realized that our
> tool - SOSlib - actually does not solve models with variable compartments
> correctly.
>
> We wonder how other tools approach this topic.
>
> In summary:
>
> SOSlib constructs ODEs as recommended by the SBML specs, e.g.:
>
> species X1 in compartment V1, with concentration x1 = X1/V1
> species X2 and X3 in V2, with conc. as above
>
> Reactions:
> X1->X2; k1*X1*V1 	// a translocation, k1 [1/sec]
> X2->X3; k2*X2*V2    // a unimolecular reaction, k2 [1/sec]
>
> Derived ODEs:
> dx1/dt = - ( k1 * X1* V1 ) / V1
> dx2/dt =  ( k1 * X1* V1 - k2 * X2 * V2 ) / V2
>
> the kinetic laws are the SBML specific conversion of `classic' rate law
> to substance/time units, instead of concentrations/time to allow for
> easy software handling of multiple compartment models.
>
>
> If however, the volumes in the model are variable with time, i.e. either
> an assignment rule or a rate rule, these equations would not be correct.
>
> A numerical integration routine would start with the initial concentration
> x(0) and add the changes during time t:
>
> x(t) = x(0) +  dx/dt * t
>
> i.e.

or

>
> x(t) = X(0)/V(0) + dx/dt * t
>
> dx/dt is calculated correctly from above ODE because it also accounts for
> changing compartments by explicit inclusion of the volume in the kinetic
> law.
>
> However, the change of the initial concentration x(0) is not accounted for
> anywhere. The solver should rather return:
>
> x(t) = X(0)/V(t) + dx/dt * t
>
>
> If x is a function of time and of volume V, x(t,V(t)),
> then the correct ODE would actually be
>
> dx/dt = dx/dt + dx/dV * dV/dt
>
> Right?

which would do the necessary adaption of the intial and current 
concentrations to a changing volume.

>
>
> We are currently discussing different solutions and their dangers
> for SOSlib, i.e. how to correctly convert the ODE system to an ODE system
> in amount/time instead of concentration/time or whether the above partial
> differential equation could or should be implemented instead.
>
> It would be helpful for us to hear how other tools handle variable
> compartments.
>
>
> And BTW, it might be good to indicate it as good practice NOT to convert
> parameters e.g. from [1/sec] to [litre/sec] to obtain correct units for
> SBML kinetic laws, but to explicitly keep the volume in the formula! This
> would allow changing the compartment size or even making the variable, at
> least for tools which do it correctly.
> E.g. the model in the specs, section 5.6 will not be valid anymore if
> someone wanted to change the volume sizes!
>
> Such rather unconvential volume-dependent rate `constants' can be
> quite dangerous!
>
>
> Rainer
>
>
>
> Rainer Machne
> Institute for Theoretical Chemistry
> University of Vienna
> Austria
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>

[SOSlib-devel] mail for SBML list: variable Compartments in SBML

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

hi,

find below a draft for bringing this issue to the SBML list, as suggested 
by Andrew.

could you review? is the problem stated clearly and correctly?


Dear SBML Community,

As we recently wanted to model variable compartments we realized that our 
tool - SOSlib - actually does not solve models with variable compartments 
correctly.

We wonder how other tools approach this topic.

In summary:

SOSlib constructs ODEs as recommended by the SBML specs, e.g.:

species X1 in compartment V1, with concentration x1 = X1/V1
species X2 and X3 in V2, with conc. as above

Reactions:
X1->X2; k1*X1*V1 	// a translocation, k1 [1/sec]
X2->X3; k2*X2*V2    // a unimolecular reaction, k2 [1/sec]

Derived ODEs:
dx1/dt = - ( k1 * X1* V1 ) / V1
dx2/dt =  ( k1 * X1* V1 - k2 * X2 * V2 ) / V2

the kinetic laws are the SBML specific conversion of `classic' rate law 
to substance/time units, instead of concentrations/time to allow for 
easy software handling of multiple compartment models.


If however, the volumes in the model are variable with time, i.e. either 
an assignment rule or a rate rule, these equations would not be correct.

A numerical integration routine would start with the initial concentration 
x(0) and add the changes during time t:

x(t) = x(0) +  dx/dt * t

i.e.

x(t) = X(0)/V(0) + dx/dt * t

dx/dt is calculated correctly from above ODE because it also accounts for 
changing compartments by explicit inclusion of the volume in the kinetic 
law.

However, the change of the initial concentration x(0) is not accounted for 
anywhere. The solver should rather return:

x(t) = X(0)/V(t) + dx/dt * t


If x is a function of time and of volume V, x(t,V(t)),
then the correct ODE would actually be

dx/dt = dx/dt + dx/dV * dV/dt

Right?


We are currently discussing different solutions and their dangers 
for SOSlib, i.e. how to correctly convert the ODE system to an ODE system 
in amount/time instead of concentration/time or whether the above partial 
differential equation could or should be implemented instead.

It would be helpful for us to hear how other tools handle variable 
compartments.


And BTW, it might be good to indicate it as good practice NOT to convert 
parameters e.g. from [1/sec] to [litre/sec] to obtain correct units for 
SBML kinetic laws, but to explicitly keep the volume in the formula! This 
would allow changing the compartment size or even making the variable, at 
least for tools which do it correctly.
E.g. the model in the specs, section 5.6 will not be valid anymore if 
someone wanted to change the volume sizes!

Such rather unconvential volume-dependent rate `constants' can be 
quite dangerous!


Rainer



Rainer Machne
Institute for Theoretical Chemistry
University of Vienna
Austria

Re: [SOSlib-devel] fixing the variable volume bug

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

> rather than have a separate array of variables for species concentration
> create a set of assignment rules that calculate the concentration

thanks, that would be a much better solution! so far the easiest from the 
implementation point of view.
also the jacobi matrix construction would automatically know about the 
necessary replacements.

And I just checked:
the semantic test suite folder contains models with variable compartments 
but they are not part of the testlist and there are no test or result 
files associated with them!

could you inquire what the reason was for not including them?


Rainer




On Tue, 20 Jun 2006, Andrew Finney wrote:

> Rainer
>
> rather than have a separate array of variables for species concentration
> create a set of assignment rules that calculate the concentration
>
> You can place identifiers for variables outside the sbml namespace by
> inserting an invalid character like '.' into the identifier string
>
> e.g. something thing like:
>
> create an assignment: s1 = s1.amount / compartment with ode as
>
> d s1.amount / dt = k * s1
>
> I'm sure this math is wrong but hopefully you get my point
>
> yours Andrew
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>

[SOSlib-devel] fixing the variable volume bug

[Andrew F. <af...@ph...>]

Rainer

rather than have a separate array of variables for species concentration
create a set of assignment rules that calculate the concentration

You can place identifiers for variables outside the sbml namespace by
inserting an invalid character like '.' into the identifier string=20

e.g. something thing like:

create an assignment: s1 =3D s1.amount / compartment with ode as

d s1.amount / dt =3D k * s1

I'm sure this math is wrong but hopefully you get my point

yours Andrew

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

[Andrew F. <af...@ph...>]

Rainer
=20
> do you know of any examples of models with a compartment with=20
> spatialDimension zero or hasOnlySubstanceUnits true?
>

Darren Wilkinson created a stoctastic test suite which consists of
models
of this form.  They should operate correctly in a deterministic
simulator as well.
=20
> can a compartment of dimension zero have a size?

No

>=20
> would we currently handle such structures correctly, at least=20
> for fixed compartment sizes?

I don't see why not

yours Andrew

> -----Original Message-----
> From: Rainer Machne [mailto:ra...@tb...]=20
> Sent: 20 June 2006 15:48
> To: Andrew Finney
> Cc: SOSlib Development
> Subject: Re: [SOSlib-devel] variable compartment bug: example[Scanned]
>=20
>=20
> >>> The units of the species are of the form substance/size units
>=20
> > "The units of the species are of the form substance/size=20
> units (i.e.,=20
> > concentration units, using a broad definition of=20
> concentration) if the=20
> > compartment's spatialDimensions is non-zero and=20
> hasOnlySubstanceUnits=20
> > has the value ``false''. The units of the species are of the form=20
> > substance if spatialDimensions is zero or hasOnlySubstanceUnits has=20
> > the value ``true''. The units of substance are those defined in the=20
> > substanceUnits, and the size units are those given in the=20
> > spatialSizeUnits field."
> > ...
>=20
> oops, i guess it was some freudian process that blinded me to=20
> the continuation of this sentence :)
>=20
> do you know of any examples of models with a compartment with=20
> spatialDimension zero or hasOnlySubstanceUnits true?
>=20
> can a compartment of dimension zero have a size?
>=20
> would we currently handle such structures correctly, at least=20
> for fixed compartment sizes?
>=20
> BTW, there are several models with variable compartments in=20
> the test suite, all in the group rulesForParametersAndCompartments.
>=20
> Could it be that SOSlib only succeeded because the program=20
> used to generate tested results also does it wrong?
>=20
> Rainer
>=20
>=20
> On Tue, 20 Jun 2006, Andrew Finney wrote:
>=20
> > Rainer
> >
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>
> >>> See section 4.6.4 here
> >>>
> >>=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level
> >> -
> >>> 2.html#SECTION00046000000000000000
> >>>
> >>> "
> >>> The units of the species are of the form substance/size units
> >>>
> >>> ...
> >>>
> >>> The units of the species are used in the following ways:
> >>>
> >>> ...
> >>>    * The species identifier has these units when it appears as a=20
> >>> numerical quantity in a mathematical formula expressed in MathML=20
> >>> (discussed in Section 3.6.3).
> >>> "
> >>>
> >>> So i guess, the units of species in math expressions are=20
> always in=20
> >>> concentrations!
> >>>
> >
> > WRONG
> >
> > They are either substance or concentration units depending on the=20
> > attributes of the species element.
> >
> > "The units of the species are of the form substance/size=20
> units (i.e.,=20
> > concentration units, using a broad definition of=20
> concentration) if the=20
> > compartment's spatialDimensions is non-zero and=20
> hasOnlySubstanceUnits=20
> > has the value ``false''. The units of the species are of the form=20
> > substance if spatialDimensions is zero or hasOnlySubstanceUnits has=20
> > the value ``true''. The units of substance are those defined in the=20
> > substanceUnits, and the size units are those given in the=20
> > spatialSizeUnits field."
> > ...
> >
> > "The units of the species are used in the following ways:
> >
> >
> > The species initialConcentration field has these units (see Section=20
> > 4.6.3).
> >
> > The species identifier has these units when it appears as a=20
> numerical=20
> > quantity in a mathematical formula expressed in MathML=20
> (discussed in=20
> > Section 3.6.3).
> >
> > The math field of an AssignmentRule structure determining=20
> the species'
> > quantity (see Section 4.8.2) has these units.
> >
> > In RateRule structures that set the rate of change of the species'
> > quantity (Section 4.8.3), the units on the rule's math=20
> field are the=20
> > units of the species divided by the built-in time units. "
> >
> > yours Andrew
> >
> >
> >
> >
> >> -----Original Message-----
> >> From: sbm...@li...
> >> [mailto:sbm...@li...] On=20
> Behalf Of=20
> >> Rainer Machne
> >> Sent: 20 June 2006 15:26
> >> To: Stefan Mueller
> >> Cc: sbm...@li...
> >> Subject: Re: [SOSlib-devel] variable compartment bug:=20
> >> example[Scanned]
> >>
> >>
> >>
> >> This
> >>
> >>> current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d=20
> S1 /dt =3D=20
> >>> f(S1, S2, ...)
> >>
> >> should of course be
> >>
> >> current ODE: d[S1]/dt =3D f([S1], [S2], ...) / V planned ODE: d
> >> S1 /dt =3D f([S1], [S2], ...)
> >>
> >>
> >> On Tue, 20 Jun 2006, Rainer Machne wrote:
> >>
> >>> stefan
> >>>
> >>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>
> >>>> yes.
> >>>> and that's exactly the same what you are planning, if you
> >> want to use
> >>>> conc ans subst in parallel.
> >>>> only the point of conversion differs.
> >>>
> >>> not really. i wasn't planning any conversion of the math,
> >> except for:
> >>>
> >>> current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d=20
> S1 /dt =3D=20
> >>> f(S1, S2, ...)
> >>>
> >>> then we have two data->value arrays, let's say
> >>>
> >>> data->odeValues =3D S1, S2, S3, etc.
> >>> and
> >>> data->concValues =3D [S1], [S2], [S3] etc.
> >>>
> >>>
> >>> The first array `odeValues' corresponds to the SUNDIALS
> >> values in NVector y.
> >>> The second array `concValues' will be used in evaluateAST, for=20
> >>> evaluation the math.
> >>>
> >>>
> >>> This would be equivalent of course, with just replacing=20
> all [S1] in=20
> >>> all math expressions with S1/volume, and only using one array
> >>> data->odeValues
> >>>
> >>> While the latter approach might be easier to implement it has a=20
> >>> probably the not so small drawback, that for each species
> >> in each math
> >>> expression the evaluation needs one more operation, i.e.
> >> the division,
> >>> while this division could be done only once for each
> >> species whenever
> >>> data->odeValues is updated (i.e. e.g. in the
> >> SOSlib/SUNDIALS funciton f).
> >>>
> >>> However, this might still only be correct, if compartments are=20
> >>> determined by assignment rules.
> >>> If they are determined by rate rules, then the replacement
> >> S1 -> S1/V
> >>> would also be required for constructing the jacobi matrix.
> >>>
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>
> >>> See section 4.6.4 here
> >>>
> >>=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level
> >> -
> >>> 2.html#SECTION00046000000000000000
> >>>
> >>> "
> >>> The units of the species are of the form substance/size units
> >>>
> >>> ...
> >>>
> >>> The units of the species are used in the following ways:
> >>>
> >>> ...
> >>>    * The species identifier has these units when it appears as a=20
> >>> numerical quantity in a mathematical formula expressed in MathML=20
> >>> (discussed in Section 3.6.3).
> >>> "
> >>>
> >>> So i guess, the units of species in math expressions are=20
> always in=20
> >>> concentrations!
> >>>
> >>> Rainer
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>>
> >>>> i mean, every species is in a well-defined compartment, isn't it?
> >>>> so you can use the formula [S]=3DS/V whenever you want...
> >>>>
> >>>> i would eliminate conc as early as possile, but that's a=20
> matter of
> >>> taste.
> >>>
> >>>
> >>>
> >>> On Tue, 20 Jun 2006, Stefan Mueller wrote:
> >>>
> >>>> hi rainer!
> >>>>
> >>>>> stefan
> >>>>>
> >>>>>>> Currently the ODE would be:
> >>>>>>>
> >>>>>>> d[S1]/dt =3D k1 * [S1] / V1
> >>>>>>
> >>>>>> this is not valid for variable compartments.
> >>>>>
> >>>>> i said that it's invalid for variable compartments some
> >> lines below:
> >>>>>>> which is only valid for constant compartments V1 and V2.
> >>>>>>
> >>>>>> please don't use it for comparison!
> >>>>>
> >>>>> it's not used for comparison, but it is the current
> >> implementation
> >>>>> in SOSlib!
> >>>>
> >>>> i said this just to avoid any additional misunderstanding.
> >>>> the whole topic seems to be quite confusing...
> >>>>
> >>>>>>> At least for the second ODE, depending in two different
> >> volumes,
> >>>>>>> the conversion to ODEs depending on substance instead of=20
> >>>>>>> concentration is not as simple as
> >>>>>>>
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>>
> >>>>>> why dou you refer to this law for a reaction of order 2?
> >>>>>
> >>>>> sorry, at first i wanted to make the second reaction in=20
> V2 of 2nd=20
> >>>>> order, e.g.,
> >>>>>
> >>>>> in V2:
> >>>>> S2 + S3 -> P: k2*[S2]*[S3]
> >>>>>
> >>>>>>> as your example before,
> >>>>>>> but for S2 I get:
> >>>>>>>
> >>>>>>> dS2/dt =3D k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
> >>>>>>>
> >>>>>>> Right?
> >>>>>>
> >>>>>> no!
> >>>>>> dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S2 it's as simple as =
that!
> >>>>>> (since dS2/dt =3D -dS1/dt as far as reaction S1->S2 is=20
> concerned.)=20
> >>>>>> that's the advantage of kinetic laws (and consequently=20
> odes) for=20
> >>>>>> substances!
> >>>>>
> >>>>> of course! sorry, i got quite confused obviously.
> >>>>>
> >>>>> and for the above second order reaction in V2 it would be:
> >>>>>
> >>>>>   dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
> >>>>>
> >>>>> Right this time?
> >>>>
> >>>> yes :)
> >>>>
> >>>>>>> Thus I would really have to check the type of reaction
> >> and edit a
> >>>>>>> lot of kinetic laws before constructing the ODEs.
> >>>>>>
> >>>>>> the above example shows that this is not necessary.
> >>>>>> maybe there are other examples...
> >>>>>
> >>>>> so you say, we could convert kinetic laws that don't explicitly=20
> >>>>> contain the volume by simply differentiating these=20
> three cases and
> >>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: k*[S1]
> >>>>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>>>
> >>>>> divide k by V(0) for 1st order reactions
> >>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *
> >> S1 * S2 *
> >>>>>>>> S3 /
> >>>>>
> >>>>> V^2
> >>>>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>
> >>>> yes.
> >>>> and that's exactly the same what you are planning, if you
> >> want to use
> >>>> conc ans subst in parallel.
> >>>> only the point of conversion differs.
> >>>>
> >>>> i mean, every species is in a well-defined compartment, isn't it?
> >>>> so you can use the formula [S]=3DS/V whenever you want...
> >>>>
> >>>> i would eliminate conc as early as possile, but that's a
> >> matter of taste.
> >>>>
> >>>>> but there are also enzymatic reactions, with in often
> >> quite complex
> >>>>> kinetic laws. we would really have to do unit checking,
> >> recognizing
> >>>>> the implemented reaction mechanisms etc.
> >>>>
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>>
> >>>>> also, we don't really know in which compartment=20
> reactions happen,=20
> >>>>> but could just try to get this information from the
> >> reactants' compartments.
> >>>>
> >>>> this does not matter, for 2 reasons:
> >>>>
> >>>> 1. if the kinetic law is specified correctly (including the=20
> >>>> compartment size), e.g. A+B->C: k*V*[A]*[B], the the
> >> compartment size
> >>>> is considered automatically.
> >>>>
> >>>> 2. the conversion [S]=3DS/V does not depend on reactions,
> >> but only on
> >>>> species (and their compartments).
> >>>>
> >>>>> so i guess, it would be much easier to
> >>>>>
> >>>>> a) use ODEs for substances but dependent on
> >> concentrations, here the
> >>>>> only thing we need to do is to leave out the compartment
> >> division at
> >>>>> the end of Species_odeFromReactions and have a double
> >> bookkeeping of
> >>>>> amounts and concentrations
> >>>>
> >>>> to my taste, double book keeping is dangerous and produces
> >> overhead.
> >>>>
> >>>>> and
> >>>>>
> >>>>> b) just ignore WRONG models that have a variable=20
> compartment but=20
> >>>>> don't have the compartment explicitly in their kinetic laws (or=20
> >>>>> check for such cases and issue an error);
> >>>>>
> >>>>> isn't it an issue of correct model writing rather then
> >> correct model
> >>>>> solving?
> >>>>
> >>>> i agree.
> >>>>
> >>>>> However, another problem might be that to calculate in amounts=20
> >>>>> rather then in concentrations also changes the use of absolute=20
> >>>>> errors and maybe also the general the performance of the
> >> integrator, i guess. ??
> >>>>
> >>>> abs errors might change, but that's not an issue.
> >>>> rel errors should stay the same.
> >>>> (in one-compartment models everything is just multiplied by the
> >>>> volume.)
> >>>>
> >>>>> maybe we will have to really implement both in parallel,
> >>>>
> >>>> i hope not :)
> >>>>
> >>>>> * ODEs for concentrations/time if no variable=20
> compartments are in=20
> >>>>> the model
> >>>>> * ODEs for ammounts/time if any compartment is variable
> >>>>>
> >>>>> ?
> >>>>>
> >>>>> Rainer
> >>>>
> >>>> maybe we discuss the next round via telephone :)
> >>>>
> >>>> cheers, s
> >>>>
> >>>>
> >>>>>>> Thus I think it just doesn't work with kinetic laws=20
> that don't=20
> >>>>>>> include the compartment explicitly!
> >>>>>>
> >>>>>> i don't know exactly what you mean.
> >>>>>> of course kinetic laws with explicit compartment are
> >> "nicer", but
> >>>>>> you can always transform any kinetic law to a "nice" one...
> >>>>>>
> >>>>>>> But if they include the compartment we can write ODEs for=20
> >>>>>>> substances but have concentrations in the right hand
> >> side of the
> >>>>>>> ODEs. We then just need to a second array for the
> >> concentrations, as discussed before.
> >>>>>>
> >>>>>> if we still use concentrations is a separate problem,
> >> independent
> >>>>>> of which we discussed until now.
> >>>>>> it's a matter of taste, if we eliminate the=20
> concentrations  from=20
> >>>>>> the kinetic laws and get a (possibly nonlinear)
> >> dependence on the
> >>>>>> compartments OR if we keep concentrations and always=20
> have linear=20
> >>>>>> dependence on the comp.
> >>>>>>
> >>>>>>> Rainer
> >>>>>>
> >>>>>> cheers,
> >>>>>> stefan.
> >>>>>>
> >>>>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
> >>>>>>>> let me summarize:
> >>>>>>>> kinetic laws have units substance/time, and they depend on=20
> >>>>>>>> concentrations or substances.
> >>>>>>>>
> >>>>>>>> we SHOULD write down odes for substances, and we COULD
> >> eliminate
> >>>>>>>> concs completely.
> >>>>>>>> (of course, we would have to convert initial concs to
> >>>>>>>> substances.)
> >>>>>>>>
> >>>>>>>> in any case, odes for substances depend on the
> >> compartment size
> >>>>>>>> (linearly or according to the order of the reaction):
> >>>>>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: k*[S1]
> >>>>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>>>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *
> >> S1 * S2 *
> >>>>>>>> S3 /
> >>>>>>>> V^2
> >>>>>>>>
> >>>>>>>> ...
> >>>>>>>>
> >>>>>>>> if the kinetic laws included the compartment
> >> explicitly (which is
> >>>>>>>> desirable), we would have the following situation:
> >>>>>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: r*[S1]*V
> >>>>>>>> dS1/dt =3D  r * [S1] * V=3D r * S1
> >>>>>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: r*[S1]*[S2]*V
> >>>>>>>> dS1/dt =3D r * [S1] * [S2] * V=3D r * S1 * S2 / V
> >>>>>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> >>>>>>>> dS1/dt =3D r * [S1] * [S2] * [S3] * V =3D r * S1 * S2 * S3 / =
V^2
> >>>>>>>>
> >>>>>>>> as stated earlier, odes for substances have got another=20
> >>>>>>>> advantage. for the compartment size we could use=20
> either a rate=20
> >>>>>>>> rule or an assignment
> >>>>>>>> rule: dV/dt =3D f(t, V, ...)
> >>>>>>>> or
> >>>>>>>> V=3D f(t, V, ...)
> >>>>>>>>
> >>>>>>>> cheers,
> >>>>>>>> stefan.
> >>>>>>>>
> >>>>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> >>>>>>>>> andrew:
> >>>>>>>>>
> >>>>>>>>> for the SBML system
> >>>>>>>>> 	S1 in C
> >>>>>>>>> 	S1->junk : k * [S1]
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>>
> >>>>>>>>> the right math is:
> >>>>>>>>> 	dS1/dt =3D k * [S1] / C(0) * C
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>> with initial cond:
> >>>>>>>>> 	S1(0) =3D [S1](0) * C(0)
> >>>>>>>>> 	C(0)
> >>>>>>>>>
> >>>>>>>>> notation:
> >>>>>>>>> 	S1 ... substance
> >>>>>>>>> 	[S1] ... conc
> >>>>>>>>>
> >>>>>>>>> the "proposed solution" (who proposed this?) has the right=20
> >>>>>>>>> units, but is wrong otherwise... :)
> >>>>>>>>>
> >>>>>>>>> as rainer and me already discussed, everything=20
> would be much=20
> >>>>>>>>> clearer, if the kinetic law included the compartment
> >> explicitly:
> >>>>>>>>> 	S1 in C
> >>>>>>>>> 	S1->junk : C * r * [S1]
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>>
> >>>>>>>>> with:
> >>>>>>>>> 	r ... the "real" chemical rate constant.
> >>>>>>>>>
> >>>>>>>>> rainer:
> >>>>>>>>>
> >>>>>>>>> chemists always knew that rate constants are not
> >> constant with
> >>>>>>>>> respect to pH, T, ... as a physicist i have to be fair! :)
> >>>>>>>>>
> >>>>>>>>> to conclude: the math is not really tricky.
> >>>>>>>>> odes for substances are better suited to handle=20
> both multiple=20
> >>>>>>>>> and variable compartments. of course, we have to=20
> take care of=20
> >>>>>>>>> many implementation details...
> >>>>>>>>>
> >>>>>>>>> cheers,
> >>>>>>>>> stefan.
> >>>>>>>>
> >>>>>>>> _______________________________________________
> >>>>>>>> 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
> >>>>
> >>>
> >>
> >>
> >> _______________________________________________
> >> 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
> >
>=20

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

[Andrew F. <af...@ph...>]

Rainer

> BTW, there are several models with variable compartments in=20
> the test suite, all in the group rulesForParametersAndCompartments.
>=20
> Could it be that SOSlib only succeeded because the program=20
> used to generate tested results also does it wrong?

Yes

I think you need to look at that model and the results
and make sure you're absolutely certain its wrong

If you are certain then I strongly urge you to make to email
SBML discuss to resolve it.

I created those tests and helped specify the behaviour of the=20
tool that created the results.  But it is possible I was mathematically
deluded.

Its obviously important that the community resolves this.

yours Andrew=20

> -----Original Message-----
> From: Rainer Machne [mailto:ra...@tb...]=20
> Sent: 20 June 2006 15:48
> To: Andrew Finney
> Cc: SOSlib Development
> Subject: Re: [SOSlib-devel] variable compartment bug: example[Scanned]
>=20
>=20
> >>> The units of the species are of the form substance/size units
>=20
> > "The units of the species are of the form substance/size=20
> units (i.e.,=20
> > concentration units, using a broad definition of=20
> concentration) if the=20
> > compartment's spatialDimensions is non-zero and=20
> hasOnlySubstanceUnits=20
> > has the value ``false''. The units of the species are of the form=20
> > substance if spatialDimensions is zero or hasOnlySubstanceUnits has=20
> > the value ``true''. The units of substance are those defined in the=20
> > substanceUnits, and the size units are those given in the=20
> > spatialSizeUnits field."
> > ...
>=20
> oops, i guess it was some freudian process that blinded me to=20
> the continuation of this sentence :)
>=20
> do you know of any examples of models with a compartment with=20
> spatialDimension zero or hasOnlySubstanceUnits true?
>=20
> can a compartment of dimension zero have a size?
>=20
> would we currently handle such structures correctly, at least=20
> for fixed compartment sizes?
>=20
> BTW, there are several models with variable compartments in=20
> the test suite, all in the group rulesForParametersAndCompartments.
>=20
> Could it be that SOSlib only succeeded because the program=20
> used to generate tested results also does it wrong?
>=20
> Rainer
>=20
>=20
> On Tue, 20 Jun 2006, Andrew Finney wrote:
>=20
> > Rainer
> >
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>
> >>> See section 4.6.4 here
> >>>
> >>=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level
> >> -
> >>> 2.html#SECTION00046000000000000000
> >>>
> >>> "
> >>> The units of the species are of the form substance/size units
> >>>
> >>> ...
> >>>
> >>> The units of the species are used in the following ways:
> >>>
> >>> ...
> >>>    * The species identifier has these units when it appears as a=20
> >>> numerical quantity in a mathematical formula expressed in MathML=20
> >>> (discussed in Section 3.6.3).
> >>> "
> >>>
> >>> So i guess, the units of species in math expressions are=20
> always in=20
> >>> concentrations!
> >>>
> >
> > WRONG
> >
> > They are either substance or concentration units depending on the=20
> > attributes of the species element.
> >
> > "The units of the species are of the form substance/size=20
> units (i.e.,=20
> > concentration units, using a broad definition of=20
> concentration) if the=20
> > compartment's spatialDimensions is non-zero and=20
> hasOnlySubstanceUnits=20
> > has the value ``false''. The units of the species are of the form=20
> > substance if spatialDimensions is zero or hasOnlySubstanceUnits has=20
> > the value ``true''. The units of substance are those defined in the=20
> > substanceUnits, and the size units are those given in the=20
> > spatialSizeUnits field."
> > ...
> >
> > "The units of the species are used in the following ways:
> >
> >
> > The species initialConcentration field has these units (see Section=20
> > 4.6.3).
> >
> > The species identifier has these units when it appears as a=20
> numerical=20
> > quantity in a mathematical formula expressed in MathML=20
> (discussed in=20
> > Section 3.6.3).
> >
> > The math field of an AssignmentRule structure determining=20
> the species'
> > quantity (see Section 4.8.2) has these units.
> >
> > In RateRule structures that set the rate of change of the species'
> > quantity (Section 4.8.3), the units on the rule's math=20
> field are the=20
> > units of the species divided by the built-in time units. "
> >
> > yours Andrew
> >
> >
> >
> >
> >> -----Original Message-----
> >> From: sbm...@li...
> >> [mailto:sbm...@li...] On=20
> Behalf Of=20
> >> Rainer Machne
> >> Sent: 20 June 2006 15:26
> >> To: Stefan Mueller
> >> Cc: sbm...@li...
> >> Subject: Re: [SOSlib-devel] variable compartment bug:=20
> >> example[Scanned]
> >>
> >>
> >>
> >> This
> >>
> >>> current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d=20
> S1 /dt =3D=20
> >>> f(S1, S2, ...)
> >>
> >> should of course be
> >>
> >> current ODE: d[S1]/dt =3D f([S1], [S2], ...) / V planned ODE: d
> >> S1 /dt =3D f([S1], [S2], ...)
> >>
> >>
> >> On Tue, 20 Jun 2006, Rainer Machne wrote:
> >>
> >>> stefan
> >>>
> >>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>
> >>>> yes.
> >>>> and that's exactly the same what you are planning, if you
> >> want to use
> >>>> conc ans subst in parallel.
> >>>> only the point of conversion differs.
> >>>
> >>> not really. i wasn't planning any conversion of the math,
> >> except for:
> >>>
> >>> current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d=20
> S1 /dt =3D=20
> >>> f(S1, S2, ...)
> >>>
> >>> then we have two data->value arrays, let's say
> >>>
> >>> data->odeValues =3D S1, S2, S3, etc.
> >>> and
> >>> data->concValues =3D [S1], [S2], [S3] etc.
> >>>
> >>>
> >>> The first array `odeValues' corresponds to the SUNDIALS
> >> values in NVector y.
> >>> The second array `concValues' will be used in evaluateAST, for=20
> >>> evaluation the math.
> >>>
> >>>
> >>> This would be equivalent of course, with just replacing=20
> all [S1] in=20
> >>> all math expressions with S1/volume, and only using one array
> >>> data->odeValues
> >>>
> >>> While the latter approach might be easier to implement it has a=20
> >>> probably the not so small drawback, that for each species
> >> in each math
> >>> expression the evaluation needs one more operation, i.e.
> >> the division,
> >>> while this division could be done only once for each
> >> species whenever
> >>> data->odeValues is updated (i.e. e.g. in the
> >> SOSlib/SUNDIALS funciton f).
> >>>
> >>> However, this might still only be correct, if compartments are=20
> >>> determined by assignment rules.
> >>> If they are determined by rate rules, then the replacement
> >> S1 -> S1/V
> >>> would also be required for constructing the jacobi matrix.
> >>>
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>
> >>> See section 4.6.4 here
> >>>
> >>=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level
> >> -
> >>> 2.html#SECTION00046000000000000000
> >>>
> >>> "
> >>> The units of the species are of the form substance/size units
> >>>
> >>> ...
> >>>
> >>> The units of the species are used in the following ways:
> >>>
> >>> ...
> >>>    * The species identifier has these units when it appears as a=20
> >>> numerical quantity in a mathematical formula expressed in MathML=20
> >>> (discussed in Section 3.6.3).
> >>> "
> >>>
> >>> So i guess, the units of species in math expressions are=20
> always in=20
> >>> concentrations!
> >>>
> >>> Rainer
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>>
> >>>> i mean, every species is in a well-defined compartment, isn't it?
> >>>> so you can use the formula [S]=3DS/V whenever you want...
> >>>>
> >>>> i would eliminate conc as early as possile, but that's a=20
> matter of
> >>> taste.
> >>>
> >>>
> >>>
> >>> On Tue, 20 Jun 2006, Stefan Mueller wrote:
> >>>
> >>>> hi rainer!
> >>>>
> >>>>> stefan
> >>>>>
> >>>>>>> Currently the ODE would be:
> >>>>>>>
> >>>>>>> d[S1]/dt =3D k1 * [S1] / V1
> >>>>>>
> >>>>>> this is not valid for variable compartments.
> >>>>>
> >>>>> i said that it's invalid for variable compartments some
> >> lines below:
> >>>>>>> which is only valid for constant compartments V1 and V2.
> >>>>>>
> >>>>>> please don't use it for comparison!
> >>>>>
> >>>>> it's not used for comparison, but it is the current
> >> implementation
> >>>>> in SOSlib!
> >>>>
> >>>> i said this just to avoid any additional misunderstanding.
> >>>> the whole topic seems to be quite confusing...
> >>>>
> >>>>>>> At least for the second ODE, depending in two different
> >> volumes,
> >>>>>>> the conversion to ODEs depending on substance instead of=20
> >>>>>>> concentration is not as simple as
> >>>>>>>
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>>
> >>>>>> why dou you refer to this law for a reaction of order 2?
> >>>>>
> >>>>> sorry, at first i wanted to make the second reaction in=20
> V2 of 2nd=20
> >>>>> order, e.g.,
> >>>>>
> >>>>> in V2:
> >>>>> S2 + S3 -> P: k2*[S2]*[S3]
> >>>>>
> >>>>>>> as your example before,
> >>>>>>> but for S2 I get:
> >>>>>>>
> >>>>>>> dS2/dt =3D k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
> >>>>>>>
> >>>>>>> Right?
> >>>>>>
> >>>>>> no!
> >>>>>> dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S2 it's as simple as =
that!
> >>>>>> (since dS2/dt =3D -dS1/dt as far as reaction S1->S2 is=20
> concerned.)=20
> >>>>>> that's the advantage of kinetic laws (and consequently=20
> odes) for=20
> >>>>>> substances!
> >>>>>
> >>>>> of course! sorry, i got quite confused obviously.
> >>>>>
> >>>>> and for the above second order reaction in V2 it would be:
> >>>>>
> >>>>>   dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
> >>>>>
> >>>>> Right this time?
> >>>>
> >>>> yes :)
> >>>>
> >>>>>>> Thus I would really have to check the type of reaction
> >> and edit a
> >>>>>>> lot of kinetic laws before constructing the ODEs.
> >>>>>>
> >>>>>> the above example shows that this is not necessary.
> >>>>>> maybe there are other examples...
> >>>>>
> >>>>> so you say, we could convert kinetic laws that don't explicitly=20
> >>>>> contain the volume by simply differentiating these=20
> three cases and
> >>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: k*[S1]
> >>>>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>>>
> >>>>> divide k by V(0) for 1st order reactions
> >>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *
> >> S1 * S2 *
> >>>>>>>> S3 /
> >>>>>
> >>>>> V^2
> >>>>>
> >>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>>>> order reactions
> >>>>
> >>>> yes.
> >>>> and that's exactly the same what you are planning, if you
> >> want to use
> >>>> conc ans subst in parallel.
> >>>> only the point of conversion differs.
> >>>>
> >>>> i mean, every species is in a well-defined compartment, isn't it?
> >>>> so you can use the formula [S]=3DS/V whenever you want...
> >>>>
> >>>> i would eliminate conc as early as possile, but that's a
> >> matter of taste.
> >>>>
> >>>>> but there are also enzymatic reactions, with in often
> >> quite complex
> >>>>> kinetic laws. we would really have to do unit checking,
> >> recognizing
> >>>>> the implemented reaction mechanisms etc.
> >>>>
> >>>> of course you have to check if the species in your kinetic
> >> law have
> >>>> dimension conc or substance!
> >>>> sbml allows the usage of both, doesn't it?
> >>>>
> >>>>> also, we don't really know in which compartment=20
> reactions happen,=20
> >>>>> but could just try to get this information from the
> >> reactants' compartments.
> >>>>
> >>>> this does not matter, for 2 reasons:
> >>>>
> >>>> 1. if the kinetic law is specified correctly (including the=20
> >>>> compartment size), e.g. A+B->C: k*V*[A]*[B], the the
> >> compartment size
> >>>> is considered automatically.
> >>>>
> >>>> 2. the conversion [S]=3DS/V does not depend on reactions,
> >> but only on
> >>>> species (and their compartments).
> >>>>
> >>>>> so i guess, it would be much easier to
> >>>>>
> >>>>> a) use ODEs for substances but dependent on
> >> concentrations, here the
> >>>>> only thing we need to do is to leave out the compartment
> >> division at
> >>>>> the end of Species_odeFromReactions and have a double
> >> bookkeeping of
> >>>>> amounts and concentrations
> >>>>
> >>>> to my taste, double book keeping is dangerous and produces
> >> overhead.
> >>>>
> >>>>> and
> >>>>>
> >>>>> b) just ignore WRONG models that have a variable=20
> compartment but=20
> >>>>> don't have the compartment explicitly in their kinetic laws (or=20
> >>>>> check for such cases and issue an error);
> >>>>>
> >>>>> isn't it an issue of correct model writing rather then
> >> correct model
> >>>>> solving?
> >>>>
> >>>> i agree.
> >>>>
> >>>>> However, another problem might be that to calculate in amounts=20
> >>>>> rather then in concentrations also changes the use of absolute=20
> >>>>> errors and maybe also the general the performance of the
> >> integrator, i guess. ??
> >>>>
> >>>> abs errors might change, but that's not an issue.
> >>>> rel errors should stay the same.
> >>>> (in one-compartment models everything is just multiplied by the
> >>>> volume.)
> >>>>
> >>>>> maybe we will have to really implement both in parallel,
> >>>>
> >>>> i hope not :)
> >>>>
> >>>>> * ODEs for concentrations/time if no variable=20
> compartments are in=20
> >>>>> the model
> >>>>> * ODEs for ammounts/time if any compartment is variable
> >>>>>
> >>>>> ?
> >>>>>
> >>>>> Rainer
> >>>>
> >>>> maybe we discuss the next round via telephone :)
> >>>>
> >>>> cheers, s
> >>>>
> >>>>
> >>>>>>> Thus I think it just doesn't work with kinetic laws=20
> that don't=20
> >>>>>>> include the compartment explicitly!
> >>>>>>
> >>>>>> i don't know exactly what you mean.
> >>>>>> of course kinetic laws with explicit compartment are
> >> "nicer", but
> >>>>>> you can always transform any kinetic law to a "nice" one...
> >>>>>>
> >>>>>>> But if they include the compartment we can write ODEs for=20
> >>>>>>> substances but have concentrations in the right hand
> >> side of the
> >>>>>>> ODEs. We then just need to a second array for the
> >> concentrations, as discussed before.
> >>>>>>
> >>>>>> if we still use concentrations is a separate problem,
> >> independent
> >>>>>> of which we discussed until now.
> >>>>>> it's a matter of taste, if we eliminate the=20
> concentrations  from=20
> >>>>>> the kinetic laws and get a (possibly nonlinear)
> >> dependence on the
> >>>>>> compartments OR if we keep concentrations and always=20
> have linear=20
> >>>>>> dependence on the comp.
> >>>>>>
> >>>>>>> Rainer
> >>>>>>
> >>>>>> cheers,
> >>>>>> stefan.
> >>>>>>
> >>>>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
> >>>>>>>> let me summarize:
> >>>>>>>> kinetic laws have units substance/time, and they depend on=20
> >>>>>>>> concentrations or substances.
> >>>>>>>>
> >>>>>>>> we SHOULD write down odes for substances, and we COULD
> >> eliminate
> >>>>>>>> concs completely.
> >>>>>>>> (of course, we would have to convert initial concs to
> >>>>>>>> substances.)
> >>>>>>>>
> >>>>>>>> in any case, odes for substances depend on the
> >> compartment size
> >>>>>>>> (linearly or according to the order of the reaction):
> >>>>>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: k*[S1]
> >>>>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>>>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / =
V
> >>>>>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *
> >> S1 * S2 *
> >>>>>>>> S3 /
> >>>>>>>> V^2
> >>>>>>>>
> >>>>>>>> ...
> >>>>>>>>
> >>>>>>>> if the kinetic laws included the compartment
> >> explicitly (which is
> >>>>>>>> desirable), we would have the following situation:
> >>>>>>>>
> >>>>>>>> order 1:
> >>>>>>>> S1->junk: r*[S1]*V
> >>>>>>>> dS1/dt =3D  r * [S1] * V=3D r * S1
> >>>>>>>>
> >>>>>>>> order 2:
> >>>>>>>> S1+S2->bla: r*[S1]*[S2]*V
> >>>>>>>> dS1/dt =3D r * [S1] * [S2] * V=3D r * S1 * S2 / V
> >>>>>>>>
> >>>>>>>> order 3:
> >>>>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> >>>>>>>> dS1/dt =3D r * [S1] * [S2] * [S3] * V =3D r * S1 * S2 * S3 / =
V^2
> >>>>>>>>
> >>>>>>>> as stated earlier, odes for substances have got another=20
> >>>>>>>> advantage. for the compartment size we could use=20
> either a rate=20
> >>>>>>>> rule or an assignment
> >>>>>>>> rule: dV/dt =3D f(t, V, ...)
> >>>>>>>> or
> >>>>>>>> V=3D f(t, V, ...)
> >>>>>>>>
> >>>>>>>> cheers,
> >>>>>>>> stefan.
> >>>>>>>>
> >>>>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> >>>>>>>>> andrew:
> >>>>>>>>>
> >>>>>>>>> for the SBML system
> >>>>>>>>> 	S1 in C
> >>>>>>>>> 	S1->junk : k * [S1]
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>>
> >>>>>>>>> the right math is:
> >>>>>>>>> 	dS1/dt =3D k * [S1] / C(0) * C
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>> with initial cond:
> >>>>>>>>> 	S1(0) =3D [S1](0) * C(0)
> >>>>>>>>> 	C(0)
> >>>>>>>>>
> >>>>>>>>> notation:
> >>>>>>>>> 	S1 ... substance
> >>>>>>>>> 	[S1] ... conc
> >>>>>>>>>
> >>>>>>>>> the "proposed solution" (who proposed this?) has the right=20
> >>>>>>>>> units, but is wrong otherwise... :)
> >>>>>>>>>
> >>>>>>>>> as rainer and me already discussed, everything=20
> would be much=20
> >>>>>>>>> clearer, if the kinetic law included the compartment
> >> explicitly:
> >>>>>>>>> 	S1 in C
> >>>>>>>>> 	S1->junk : C * r * [S1]
> >>>>>>>>> 	dC/dt =3D 1
> >>>>>>>>>
> >>>>>>>>> with:
> >>>>>>>>> 	r ... the "real" chemical rate constant.
> >>>>>>>>>
> >>>>>>>>> rainer:
> >>>>>>>>>
> >>>>>>>>> chemists always knew that rate constants are not
> >> constant with
> >>>>>>>>> respect to pH, T, ... as a physicist i have to be fair! :)
> >>>>>>>>>
> >>>>>>>>> to conclude: the math is not really tricky.
> >>>>>>>>> odes for substances are better suited to handle=20
> both multiple=20
> >>>>>>>>> and variable compartments. of course, we have to=20
> take care of=20
> >>>>>>>>> many implementation details...
> >>>>>>>>>
> >>>>>>>>> cheers,
> >>>>>>>>> stefan.
> >>>>>>>>
> >>>>>>>> _______________________________________________
> >>>>>>>> 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
> >>>>
> >>>
> >>
> >>
> >> _______________________________________________
> >> 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
> >
>=20

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

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

>>> The units of the species are of the form substance/size units

> "The units of the species are of the form substance/size units (i.e.,
> concentration units, using a broad definition of concentration) if the
> compartment's spatialDimensions is non-zero and hasOnlySubstanceUnits
> has the value ``false''. The units of the species are of the form
> substance if spatialDimensions is zero or hasOnlySubstanceUnits has the
> value ``true''. The units of substance are those defined in the
> substanceUnits, and the size units are those given in the
> spatialSizeUnits field."
> ...

oops, i guess it was some freudian process that blinded me to the 
continuation of this sentence :)

do you know of any examples of models with a compartment 
with spatialDimension zero or hasOnlySubstanceUnits true?

can a compartment of dimension zero have a size?

would we currently handle such structures correctly, at least for fixed 
compartment sizes?

BTW, there are several models with variable compartments in the test 
suite, all in the group rulesForParametersAndCompartments.

Could it be that SOSlib only succeeded because the program used to 
generate tested results also does it wrong?

Rainer


On Tue, 20 Jun 2006, Andrew Finney wrote:

> Rainer
>
>>>> of course you have to check if the species in your kinetic
>> law have
>>>> dimension conc or substance!
>>>> sbml allows the usage of both, doesn't it?
>>>
>>> See section 4.6.4 here
>>>
>> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-
>>> 2.html#SECTION00046000000000000000
>>>
>>> "
>>> The units of the species are of the form substance/size units
>>>
>>> ...
>>>
>>> The units of the species are used in the following ways:
>>>
>>> ...
>>>    * The species identifier has these units when it appears as a
>>> numerical quantity in a mathematical formula expressed in MathML
>>> (discussed in Section 3.6.3).
>>> "
>>>
>>> So i guess, the units of species in math expressions are always in
>>> concentrations!
>>>
>
> WRONG
>
> They are either substance or concentration units depending on the
> attributes of the species element.
>
> "The units of the species are of the form substance/size units (i.e.,
> concentration units, using a broad definition of concentration) if the
> compartment's spatialDimensions is non-zero and hasOnlySubstanceUnits
> has the value ``false''. The units of the species are of the form
> substance if spatialDimensions is zero or hasOnlySubstanceUnits has the
> value ``true''. The units of substance are those defined in the
> substanceUnits, and the size units are those given in the
> spatialSizeUnits field."
> ...
>
> "The units of the species are used in the following ways:
>
>
> The species initialConcentration field has these units (see Section
> 4.6.3).
>
> The species identifier has these units when it appears as a numerical
> quantity in a mathematical formula expressed in MathML (discussed in
> Section 3.6.3).
>
> The math field of an AssignmentRule structure determining the species'
> quantity (see Section 4.8.2) has these units.
>
> In RateRule structures that set the rate of change of the species'
> quantity (Section 4.8.3), the units on the rule's math field are the
> units of the species divided by the built-in time units. "
>
> yours Andrew
>
>
>
>
>> -----Original Message-----
>> From: sbm...@li...
>> [mailto:sbm...@li...] On
>> Behalf Of Rainer Machne
>> Sent: 20 June 2006 15:26
>> To: Stefan Mueller
>> Cc: sbm...@li...
>> Subject: Re: [SOSlib-devel] variable compartment bug: example[Scanned]
>>
>>
>>
>> This
>>
>>> current ODE: d[S1]/dt = f(S1, S2, ...) / V planned ODE: d S1 /dt =
>>> f(S1, S2, ...)
>>
>> should of course be
>>
>> current ODE: d[S1]/dt = f([S1], [S2], ...) / V planned ODE: d
>> S1 /dt = f([S1], [S2], ...)
>>
>>
>> On Tue, 20 Jun 2006, Rainer Machne wrote:
>>
>>> stefan
>>>
>>>
>>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd
>>>>> order reactions
>>>>
>>>> yes.
>>>> and that's exactly the same what you are planning, if you
>> want to use
>>>> conc ans subst in parallel.
>>>> only the point of conversion differs.
>>>
>>> not really. i wasn't planning any conversion of the math,
>> except for:
>>>
>>> current ODE: d[S1]/dt = f(S1, S2, ...) / V planned ODE: d S1 /dt =
>>> f(S1, S2, ...)
>>>
>>> then we have two data->value arrays, let's say
>>>
>>> data->odeValues = S1, S2, S3, etc.
>>> and
>>> data->concValues = [S1], [S2], [S3] etc.
>>>
>>>
>>> The first array `odeValues' corresponds to the SUNDIALS
>> values in NVector y.
>>> The second array `concValues' will be used in evaluateAST, for
>>> evaluation the math.
>>>
>>>
>>> This would be equivalent of course, with just replacing all [S1] in
>>> all math expressions with S1/volume, and only using one array
>>> data->odeValues
>>>
>>> While the latter approach might be easier to implement it has a
>>> probably the not so small drawback, that for each species
>> in each math
>>> expression the evaluation needs one more operation, i.e.
>> the division,
>>> while this division could be done only once for each
>> species whenever
>>> data->odeValues is updated (i.e. e.g. in the
>> SOSlib/SUNDIALS funciton f).
>>>
>>> However, this might still only be correct, if compartments are
>>> determined by assignment rules.
>>> If they are determined by rate rules, then the replacement
>> S1 -> S1/V
>>> would also be required for constructing the jacobi matrix.
>>>
>>>> of course you have to check if the species in your kinetic
>> law have
>>>> dimension conc or substance!
>>>> sbml allows the usage of both, doesn't it?
>>>
>>> See section 4.6.4 here
>>>
>> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-
>>> 2.html#SECTION00046000000000000000
>>>
>>> "
>>> The units of the species are of the form substance/size units
>>>
>>> ...
>>>
>>> The units of the species are used in the following ways:
>>>
>>> ...
>>>    * The species identifier has these units when it appears as a
>>> numerical quantity in a mathematical formula expressed in MathML
>>> (discussed in Section 3.6.3).
>>> "
>>>
>>> So i guess, the units of species in math expressions are always in
>>> concentrations!
>>>
>>> Rainer
>>>
>>>
>>>
>>>
>>>
>>>>
>>>> i mean, every species is in a well-defined compartment, isn't it?
>>>> so you can use the formula [S]=S/V whenever you want...
>>>>
>>>> i would eliminate conc as early as possile, but that's a matter of
>>> taste.
>>>
>>>
>>>
>>> On Tue, 20 Jun 2006, Stefan Mueller wrote:
>>>
>>>> hi rainer!
>>>>
>>>>> stefan
>>>>>
>>>>>>> Currently the ODE would be:
>>>>>>>
>>>>>>> d[S1]/dt = k1 * [S1] / V1
>>>>>>
>>>>>> this is not valid for variable compartments.
>>>>>
>>>>> i said that it's invalid for variable compartments some
>> lines below:
>>>>>>> which is only valid for constant compartments V1 and V2.
>>>>>>
>>>>>> please don't use it for comparison!
>>>>>
>>>>> it's not used for comparison, but it is the current
>> implementation
>>>>> in SOSlib!
>>>>
>>>> i said this just to avoid any additional misunderstanding.
>>>> the whole topic seems to be quite confusing...
>>>>
>>>>>>> At least for the second ODE, depending in two different
>> volumes,
>>>>>>> the conversion to ODEs depending on substance instead of
>>>>>>> concentration is not as simple as
>>>>>>>
>>>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>>
>>>>>> why dou you refer to this law for a reaction of order 2?
>>>>>
>>>>> sorry, at first i wanted to make the second reaction in V2 of 2nd
>>>>> order, e.g.,
>>>>>
>>>>> in V2:
>>>>> S2 + S3 -> P: k2*[S2]*[S3]
>>>>>
>>>>>>> as your example before,
>>>>>>> but for S2 I get:
>>>>>>>
>>>>>>> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>>>>>>>
>>>>>>> Right?
>>>>>>
>>>>>> no!
>>>>>> dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2 it's as simple as that!
>>>>>> (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
>>>>>> that's the advantage of kinetic laws (and consequently odes) for
>>>>>> substances!
>>>>>
>>>>> of course! sorry, i got quite confused obviously.
>>>>>
>>>>> and for the above second order reaction in V2 it would be:
>>>>>
>>>>>   dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
>>>>>
>>>>> Right this time?
>>>>
>>>> yes :)
>>>>
>>>>>>> Thus I would really have to check the type of reaction
>> and edit a
>>>>>>> lot of kinetic laws before constructing the ODEs.
>>>>>>
>>>>>> the above example shows that this is not necessary.
>>>>>> maybe there are other examples...
>>>>>
>>>>> so you say, we could convert kinetic laws that don't explicitly
>>>>> contain the volume by simply differentiating these three cases and
>>>>>
>>>>>>>> order 1:
>>>>>>>> S1->junk: k*[S1]
>>>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>>>
>>>>> divide k by V(0) for 1st order reactions
>>>>>
>>>>>>>> order 2:
>>>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>
>>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd
>>>>> order reactions
>>>>>
>>>>>>>> order 3:
>>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) *
>> S1 * S2 *
>>>>>>>> S3 /
>>>>>
>>>>> V^2
>>>>>
>>>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd
>>>>> order reactions
>>>>
>>>> yes.
>>>> and that's exactly the same what you are planning, if you
>> want to use
>>>> conc ans subst in parallel.
>>>> only the point of conversion differs.
>>>>
>>>> i mean, every species is in a well-defined compartment, isn't it?
>>>> so you can use the formula [S]=S/V whenever you want...
>>>>
>>>> i would eliminate conc as early as possile, but that's a
>> matter of taste.
>>>>
>>>>> but there are also enzymatic reactions, with in often
>> quite complex
>>>>> kinetic laws. we would really have to do unit checking,
>> recognizing
>>>>> the implemented reaction mechanisms etc.
>>>>
>>>> of course you have to check if the species in your kinetic
>> law have
>>>> dimension conc or substance!
>>>> sbml allows the usage of both, doesn't it?
>>>>
>>>>> also, we don't really know in which compartment reactions happen,
>>>>> but could just try to get this information from the
>> reactants' compartments.
>>>>
>>>> this does not matter, for 2 reasons:
>>>>
>>>> 1. if the kinetic law is specified correctly (including the
>>>> compartment size), e.g. A+B->C: k*V*[A]*[B], the the
>> compartment size
>>>> is considered automatically.
>>>>
>>>> 2. the conversion [S]=S/V does not depend on reactions,
>> but only on
>>>> species (and their compartments).
>>>>
>>>>> so i guess, it would be much easier to
>>>>>
>>>>> a) use ODEs for substances but dependent on
>> concentrations, here the
>>>>> only thing we need to do is to leave out the compartment
>> division at
>>>>> the end of Species_odeFromReactions and have a double
>> bookkeeping of
>>>>> amounts and concentrations
>>>>
>>>> to my taste, double book keeping is dangerous and produces
>> overhead.
>>>>
>>>>> and
>>>>>
>>>>> b) just ignore WRONG models that have a variable compartment but
>>>>> don't have the compartment explicitly in their kinetic laws (or
>>>>> check for such cases and issue an error);
>>>>>
>>>>> isn't it an issue of correct model writing rather then
>> correct model
>>>>> solving?
>>>>
>>>> i agree.
>>>>
>>>>> However, another problem might be that to calculate in amounts
>>>>> rather then in concentrations also changes the use of absolute
>>>>> errors and maybe also the general the performance of the
>> integrator, i guess. ??
>>>>
>>>> abs errors might change, but that's not an issue.
>>>> rel errors should stay the same.
>>>> (in one-compartment models everything is just multiplied by the
>>>> volume.)
>>>>
>>>>> maybe we will have to really implement both in parallel,
>>>>
>>>> i hope not :)
>>>>
>>>>> * ODEs for concentrations/time if no variable compartments are in
>>>>> the model
>>>>> * ODEs for ammounts/time if any compartment is variable
>>>>>
>>>>> ?
>>>>>
>>>>> Rainer
>>>>
>>>> maybe we discuss the next round via telephone :)
>>>>
>>>> cheers, s
>>>>
>>>>
>>>>>>> Thus I think it just doesn't work with kinetic laws that don't
>>>>>>> include the compartment explicitly!
>>>>>>
>>>>>> i don't know exactly what you mean.
>>>>>> of course kinetic laws with explicit compartment are
>> "nicer", but
>>>>>> you can always transform any kinetic law to a "nice" one...
>>>>>>
>>>>>>> But if they include the compartment we can write ODEs for
>>>>>>> substances but have concentrations in the right hand
>> side of the
>>>>>>> ODEs. We then just need to a second array for the
>> concentrations, as discussed before.
>>>>>>
>>>>>> if we still use concentrations is a separate problem,
>> independent
>>>>>> of which we discussed until now.
>>>>>> it's a matter of taste, if we eliminate the concentrations  from
>>>>>> the kinetic laws and get a (possibly nonlinear)
>> dependence on the
>>>>>> compartments OR if we keep concentrations and always have linear
>>>>>> dependence on the comp.
>>>>>>
>>>>>>> Rainer
>>>>>>
>>>>>> cheers,
>>>>>> stefan.
>>>>>>
>>>>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
>>>>>>>> let me summarize:
>>>>>>>> kinetic laws have units substance/time, and they depend on
>>>>>>>> concentrations or substances.
>>>>>>>>
>>>>>>>> we SHOULD write down odes for substances, and we COULD
>> eliminate
>>>>>>>> concs completely.
>>>>>>>> (of course, we would have to convert initial concs to
>>>>>>>> substances.)
>>>>>>>>
>>>>>>>> in any case, odes for substances depend on the
>> compartment size
>>>>>>>> (linearly or according to the order of the reaction):
>>>>>>>>
>>>>>>>> order 1:
>>>>>>>> S1->junk: k*[S1]
>>>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>>>>>>
>>>>>>>> order 2:
>>>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>>>>
>>>>>>>> order 3:
>>>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) *
>> S1 * S2 *
>>>>>>>> S3 /
>>>>>>>> V^2
>>>>>>>>
>>>>>>>> ...
>>>>>>>>
>>>>>>>> if the kinetic laws included the compartment
>> explicitly (which is
>>>>>>>> desirable), we would have the following situation:
>>>>>>>>
>>>>>>>> order 1:
>>>>>>>> S1->junk: r*[S1]*V
>>>>>>>> dS1/dt =  r * [S1] * V= r * S1
>>>>>>>>
>>>>>>>> order 2:
>>>>>>>> S1+S2->bla: r*[S1]*[S2]*V
>>>>>>>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>>>>>>>>
>>>>>>>> order 3:
>>>>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
>>>>>>>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>>>>>>>>
>>>>>>>> as stated earlier, odes for substances have got another
>>>>>>>> advantage. for the compartment size we could use either a rate
>>>>>>>> rule or an assignment
>>>>>>>> rule: dV/dt = f(t, V, ...)
>>>>>>>> or
>>>>>>>> V= f(t, V, ...)
>>>>>>>>
>>>>>>>> cheers,
>>>>>>>> stefan.
>>>>>>>>
>>>>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>>>>>>>>> andrew:
>>>>>>>>>
>>>>>>>>> for the SBML system
>>>>>>>>> 	S1 in C
>>>>>>>>> 	S1->junk : k * [S1]
>>>>>>>>> 	dC/dt = 1
>>>>>>>>>
>>>>>>>>> the right math is:
>>>>>>>>> 	dS1/dt = k * [S1] / C(0) * C
>>>>>>>>> 	dC/dt = 1
>>>>>>>>> with initial cond:
>>>>>>>>> 	S1(0) = [S1](0) * C(0)
>>>>>>>>> 	C(0)
>>>>>>>>>
>>>>>>>>> notation:
>>>>>>>>> 	S1 ... substance
>>>>>>>>> 	[S1] ... conc
>>>>>>>>>
>>>>>>>>> the "proposed solution" (who proposed this?) has the right
>>>>>>>>> units, but is wrong otherwise... :)
>>>>>>>>>
>>>>>>>>> as rainer and me already discussed, everything would be much
>>>>>>>>> clearer, if the kinetic law included the compartment
>> explicitly:
>>>>>>>>> 	S1 in C
>>>>>>>>> 	S1->junk : C * r * [S1]
>>>>>>>>> 	dC/dt = 1
>>>>>>>>>
>>>>>>>>> with:
>>>>>>>>> 	r ... the "real" chemical rate constant.
>>>>>>>>>
>>>>>>>>> rainer:
>>>>>>>>>
>>>>>>>>> chemists always knew that rate constants are not
>> constant with
>>>>>>>>> respect to pH, T, ... as a physicist i have to be fair! :)
>>>>>>>>>
>>>>>>>>> to conclude: the math is not really tricky.
>>>>>>>>> odes for substances are better suited to handle both multiple
>>>>>>>>> and variable compartments. of course, we have to take care of
>>>>>>>>> many implementation details...
>>>>>>>>>
>>>>>>>>> cheers,
>>>>>>>>> stefan.
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> 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
>>>>
>>>
>>
>>
>> _______________________________________________
>> 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
>

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

[Andrew F. <af...@ph...>]

Rainer

> >> of course you have to check if the species in your kinetic=20
> law have=20
> >> dimension conc or substance!
> >> sbml allows the usage of both, doesn't it?
> >
> > See section 4.6.4 here
> >=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-
> > 2.html#SECTION00046000000000000000
> >
> > "
> > The units of the species are of the form substance/size units
> >
> > ...
> >
> > The units of the species are used in the following ways:
> >
> > ...
> >    * The species identifier has these units when it appears as a=20
> > numerical quantity in a mathematical formula expressed in MathML=20
> > (discussed in Section 3.6.3).
> > "
> >
> > So i guess, the units of species in math expressions are always in=20
> > concentrations!
> >

WRONG

They are either substance or concentration units depending on the
attributes of the species element.

"The units of the species are of the form substance/size units (i.e.,
concentration units, using a broad definition of concentration) if the
compartment's spatialDimensions is non-zero and hasOnlySubstanceUnits
has the value ``false''. The units of the species are of the form
substance if spatialDimensions is zero or hasOnlySubstanceUnits has the
value ``true''. The units of substance are those defined in the
substanceUnits, and the size units are those given in the
spatialSizeUnits field."
...

"The units of the species are used in the following ways:=20


The species initialConcentration field has these units (see Section
4.6.3).=20

The species identifier has these units when it appears as a numerical
quantity in a mathematical formula expressed in MathML (discussed in
Section 3.6.3).=20

The math field of an AssignmentRule structure determining the species'
quantity (see Section 4.8.2) has these units.=20

In RateRule structures that set the rate of change of the species'
quantity (Section 4.8.3), the units on the rule's math field are the
units of the species divided by the built-in time units. "

yours Andrew


=20

> -----Original Message-----
> From: sbm...@li...=20
> [mailto:sbm...@li...] On=20
> Behalf Of Rainer Machne
> Sent: 20 June 2006 15:26
> To: Stefan Mueller
> Cc: sbm...@li...
> Subject: Re: [SOSlib-devel] variable compartment bug: example[Scanned]
>=20
>=20
>=20
> This
>=20
> > current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d S1 /dt =
=3D=20
> > f(S1, S2, ...)
>=20
> should of course be
>=20
> current ODE: d[S1]/dt =3D f([S1], [S2], ...) / V planned ODE: d=20
> S1 /dt =3D f([S1], [S2], ...)
>=20
>=20
> On Tue, 20 Jun 2006, Rainer Machne wrote:
>=20
> > stefan
> >
> >
> >>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>> order reactions
> >>=20
> >> yes.
> >> and that's exactly the same what you are planning, if you=20
> want to use=20
> >> conc ans subst in parallel.
> >> only the point of conversion differs.
> >
> > not really. i wasn't planning any conversion of the math,=20
> except for:
> >
> > current ODE: d[S1]/dt =3D f(S1, S2, ...) / V planned ODE: d S1 /dt =
=3D=20
> > f(S1, S2, ...)
> >
> > then we have two data->value arrays, let's say
> >
> > data->odeValues =3D S1, S2, S3, etc.
> > and
> > data->concValues =3D [S1], [S2], [S3] etc.
> >
> >
> > The first array `odeValues' corresponds to the SUNDIALS=20
> values in NVector y.
> > The second array `concValues' will be used in evaluateAST, for=20
> > evaluation the math.
> >
> >
> > This would be equivalent of course, with just replacing all [S1] in=20
> > all math expressions with S1/volume, and only using one array=20
> > data->odeValues
> >
> > While the latter approach might be easier to implement it has a=20
> > probably the not so small drawback, that for each species=20
> in each math=20
> > expression the evaluation needs one more operation, i.e.=20
> the division,=20
> > while this division could be done only once for each=20
> species whenever=20
> > data->odeValues is updated (i.e. e.g. in the=20
> SOSlib/SUNDIALS funciton f).
> >
> > However, this might still only be correct, if compartments are=20
> > determined by assignment rules.
> > If they are determined by rate rules, then the replacement=20
> S1 -> S1/V=20
> > would also be required for constructing the jacobi matrix.
> >
> >> of course you have to check if the species in your kinetic=20
> law have=20
> >> dimension conc or substance!
> >> sbml allows the usage of both, doesn't it?
> >
> > See section 4.6.4 here
> >=20
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-
> > 2.html#SECTION00046000000000000000
> >
> > "
> > The units of the species are of the form substance/size units
> >
> > ...
> >
> > The units of the species are used in the following ways:
> >
> > ...
> >    * The species identifier has these units when it appears as a=20
> > numerical quantity in a mathematical formula expressed in MathML=20
> > (discussed in Section 3.6.3).
> > "
> >
> > So i guess, the units of species in math expressions are always in=20
> > concentrations!
> >
> > Rainer
> >
> >
> >
> >
> >
> >>=20
> >> i mean, every species is in a well-defined compartment, isn't it?
> >> so you can use the formula [S]=3DS/V whenever you want...
> >>=20
> >> i would eliminate conc as early as possile, but that's a matter of
> > taste.
> >
> >
> >
> > On Tue, 20 Jun 2006, Stefan Mueller wrote:
> >
> >> hi rainer!
> >>=20
> >>> stefan
> >>>=20
> >>>>> Currently the ODE would be:
> >>>>>=20
> >>>>> d[S1]/dt =3D k1 * [S1] / V1
> >>>>=20
> >>>> this is not valid for variable compartments.
> >>>=20
> >>> i said that it's invalid for variable compartments some=20
> lines below:
> >>>>> which is only valid for constant compartments V1 and V2.
> >>>>=20
> >>>> please don't use it for comparison!
> >>>=20
> >>> it's not used for comparison, but it is the current=20
> implementation=20
> >>> in SOSlib!
> >>=20
> >> i said this just to avoid any additional misunderstanding.
> >> the whole topic seems to be quite confusing...
> >>=20
> >>>>> At least for the second ODE, depending in two different=20
> volumes,=20
> >>>>> the conversion to ODEs depending on substance instead of=20
> >>>>> concentration is not as simple as
> >>>>>=20
> >>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / V
> >>>>=20
> >>>> why dou you refer to this law for a reaction of order 2?
> >>>=20
> >>> sorry, at first i wanted to make the second reaction in V2 of 2nd=20
> >>> order, e.g.,
> >>>=20
> >>> in V2:
> >>> S2 + S3 -> P: k2*[S2]*[S3]
> >>>=20
> >>>>> as your example before,
> >>>>> but for S2 I get:
> >>>>>=20
> >>>>> dS2/dt =3D k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
> >>>>>=20
> >>>>> Right?
> >>>>=20
> >>>> no!
> >>>> dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S2 it's as simple as that!
> >>>> (since dS2/dt =3D -dS1/dt as far as reaction S1->S2 is =
concerned.)=20
> >>>> that's the advantage of kinetic laws (and consequently odes) for=20
> >>>> substances!
> >>>=20
> >>> of course! sorry, i got quite confused obviously.
> >>>=20
> >>> and for the above second order reaction in V2 it would be:
> >>>=20
> >>>   dS2/dt =3D k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
> >>>=20
> >>> Right this time?
> >>=20
> >> yes :)
> >>=20
> >>>>> Thus I would really have to check the type of reaction=20
> and edit a=20
> >>>>> lot of kinetic laws before constructing the ODEs.
> >>>>=20
> >>>> the above example shows that this is not necessary.
> >>>> maybe there are other examples...
> >>>=20
> >>> so you say, we could convert kinetic laws that don't explicitly=20
> >>> contain the volume by simply differentiating these three cases and
> >>>=20
> >>>>>> order 1:
> >>>>>> S1->junk: k*[S1]
> >>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>=20
> >>> divide k by V(0) for 1st order reactions
> >>>=20
> >>>>>> order 2:
> >>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / V
> >>>=20
> >>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>> order reactions
> >>>=20
> >>>>>> order 3:
> >>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *=20
> S1 * S2 *=20
> >>>>>> S3 /
> >>>=20
> >>> V^2
> >>>=20
> >>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd=20
> >>> order reactions
> >>=20
> >> yes.
> >> and that's exactly the same what you are planning, if you=20
> want to use=20
> >> conc ans subst in parallel.
> >> only the point of conversion differs.
> >>=20
> >> i mean, every species is in a well-defined compartment, isn't it?
> >> so you can use the formula [S]=3DS/V whenever you want...
> >>=20
> >> i would eliminate conc as early as possile, but that's a=20
> matter of taste.
> >>=20
> >>> but there are also enzymatic reactions, with in often=20
> quite complex=20
> >>> kinetic laws. we would really have to do unit checking,=20
> recognizing=20
> >>> the implemented reaction mechanisms etc.
> >>=20
> >> of course you have to check if the species in your kinetic=20
> law have=20
> >> dimension conc or substance!
> >> sbml allows the usage of both, doesn't it?
> >>=20
> >>> also, we don't really know in which compartment reactions happen,=20
> >>> but could just try to get this information from the=20
> reactants' compartments.
> >>=20
> >> this does not matter, for 2 reasons:
> >>=20
> >> 1. if the kinetic law is specified correctly (including the=20
> >> compartment size), e.g. A+B->C: k*V*[A]*[B], the the=20
> compartment size=20
> >> is considered automatically.
> >>=20
> >> 2. the conversion [S]=3DS/V does not depend on reactions,=20
> but only on=20
> >> species (and their compartments).
> >>=20
> >>> so i guess, it would be much easier to
> >>>=20
> >>> a) use ODEs for substances but dependent on=20
> concentrations, here the=20
> >>> only thing we need to do is to leave out the compartment=20
> division at=20
> >>> the end of Species_odeFromReactions and have a double=20
> bookkeeping of=20
> >>> amounts and concentrations
> >>=20
> >> to my taste, double book keeping is dangerous and produces=20
> overhead.
> >>=20
> >>> and
> >>>=20
> >>> b) just ignore WRONG models that have a variable compartment but=20
> >>> don't have the compartment explicitly in their kinetic laws (or=20
> >>> check for such cases and issue an error);
> >>>=20
> >>> isn't it an issue of correct model writing rather then=20
> correct model=20
> >>> solving?
> >>=20
> >> i agree.
> >>=20
> >>> However, another problem might be that to calculate in amounts=20
> >>> rather then in concentrations also changes the use of absolute=20
> >>> errors and maybe also the general the performance of the=20
> integrator, i guess. ??
> >>=20
> >> abs errors might change, but that's not an issue.
> >> rel errors should stay the same.
> >> (in one-compartment models everything is just multiplied by the=20
> >> volume.)
> >>=20
> >>> maybe we will have to really implement both in parallel,
> >>=20
> >> i hope not :)
> >>=20
> >>> * ODEs for concentrations/time if no variable compartments are in=20
> >>> the model
> >>> * ODEs for ammounts/time if any compartment is variable
> >>>=20
> >>> ?
> >>>=20
> >>> Rainer
> >>=20
> >> maybe we discuss the next round via telephone :)
> >>=20
> >> cheers, s
> >>=20
> >>=20
> >>>>> Thus I think it just doesn't work with kinetic laws that don't=20
> >>>>> include the compartment explicitly!
> >>>>=20
> >>>> i don't know exactly what you mean.
> >>>> of course kinetic laws with explicit compartment are=20
> "nicer", but=20
> >>>> you can always transform any kinetic law to a "nice" one...
> >>>>=20
> >>>>> But if they include the compartment we can write ODEs for=20
> >>>>> substances but have concentrations in the right hand=20
> side of the=20
> >>>>> ODEs. We then just need to a second array for the=20
> concentrations, as discussed before.
> >>>>=20
> >>>> if we still use concentrations is a separate problem,=20
> independent=20
> >>>> of which we discussed until now.
> >>>> it's a matter of taste, if we eliminate the concentrations  from=20
> >>>> the kinetic laws and get a (possibly nonlinear)=20
> dependence on the=20
> >>>> compartments OR if we keep concentrations and always have linear=20
> >>>> dependence on the comp.
> >>>>=20
> >>>>> Rainer
> >>>>=20
> >>>> cheers,
> >>>> stefan.
> >>>>=20
> >>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
> >>>>>> let me summarize:
> >>>>>> kinetic laws have units substance/time, and they depend on=20
> >>>>>> concentrations or substances.
> >>>>>>=20
> >>>>>> we SHOULD write down odes for substances, and we COULD=20
> eliminate=20
> >>>>>> concs completely.
> >>>>>> (of course, we would have to convert initial concs to=20
> >>>>>> substances.)
> >>>>>>=20
> >>>>>> in any case, odes for substances depend on the=20
> compartment size=20
> >>>>>> (linearly or according to the order of the reaction):
> >>>>>>=20
> >>>>>> order 1:
> >>>>>> S1->junk: k*[S1]
> >>>>>> dS1/dt =3D k * [S1] / V(0) * V =3D k/V(0) * S1
> >>>>>>=20
> >>>>>> order 2:
> >>>>>> S1+S2->bla: k*[S1]*[S2]
> >>>>>> dS1/dt =3D k * [S1] * [S2] / V(0) * V =3D k/V(0) * S1 * S2 / V
> >>>>>>=20
> >>>>>> order 3:
> >>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>>>>> dS1/dt =3D k * [S1] * [S2] * [S3] / V(0) * V =3D k/V(0) *=20
> S1 * S2 *=20
> >>>>>> S3 /
> >>>>>> V^2
> >>>>>>=20
> >>>>>> ...
> >>>>>>=20
> >>>>>> if the kinetic laws included the compartment=20
> explicitly (which is=20
> >>>>>> desirable), we would have the following situation:
> >>>>>>=20
> >>>>>> order 1:
> >>>>>> S1->junk: r*[S1]*V
> >>>>>> dS1/dt =3D  r * [S1] * V=3D r * S1
> >>>>>>=20
> >>>>>> order 2:
> >>>>>> S1+S2->bla: r*[S1]*[S2]*V
> >>>>>> dS1/dt =3D r * [S1] * [S2] * V=3D r * S1 * S2 / V
> >>>>>>=20
> >>>>>> order 3:
> >>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> >>>>>> dS1/dt =3D r * [S1] * [S2] * [S3] * V =3D r * S1 * S2 * S3 / =
V^2
> >>>>>>=20
> >>>>>> as stated earlier, odes for substances have got another=20
> >>>>>> advantage. for the compartment size we could use either a rate=20
> >>>>>> rule or an assignment
> >>>>>> rule: dV/dt =3D f(t, V, ...)
> >>>>>> or
> >>>>>> V=3D f(t, V, ...)
> >>>>>>=20
> >>>>>> cheers,
> >>>>>> stefan.
> >>>>>>=20
> >>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> >>>>>>> andrew:
> >>>>>>>=20
> >>>>>>> for the SBML system
> >>>>>>> 	S1 in C
> >>>>>>> 	S1->junk : k * [S1]
> >>>>>>> 	dC/dt =3D 1
> >>>>>>>=20
> >>>>>>> the right math is:
> >>>>>>> 	dS1/dt =3D k * [S1] / C(0) * C
> >>>>>>> 	dC/dt =3D 1
> >>>>>>> with initial cond:
> >>>>>>> 	S1(0) =3D [S1](0) * C(0)
> >>>>>>> 	C(0)
> >>>>>>>=20
> >>>>>>> notation:
> >>>>>>> 	S1 ... substance
> >>>>>>> 	[S1] ... conc
> >>>>>>>=20
> >>>>>>> the "proposed solution" (who proposed this?) has the right=20
> >>>>>>> units, but is wrong otherwise... :)
> >>>>>>>=20
> >>>>>>> as rainer and me already discussed, everything would be much=20
> >>>>>>> clearer, if the kinetic law included the compartment=20
> explicitly:
> >>>>>>> 	S1 in C
> >>>>>>> 	S1->junk : C * r * [S1]
> >>>>>>> 	dC/dt =3D 1
> >>>>>>>=20
> >>>>>>> with:
> >>>>>>> 	r ... the "real" chemical rate constant.
> >>>>>>>=20
> >>>>>>> rainer:
> >>>>>>>=20
> >>>>>>> chemists always knew that rate constants are not=20
> constant with=20
> >>>>>>> respect to pH, T, ... as a physicist i have to be fair! :)
> >>>>>>>=20
> >>>>>>> to conclude: the math is not really tricky.
> >>>>>>> odes for substances are better suited to handle both multiple=20
> >>>>>>> and variable compartments. of course, we have to take care of=20
> >>>>>>> many implementation details...
> >>>>>>>=20
> >>>>>>> cheers,
> >>>>>>> stefan.
> >>>>>>=20
> >>>>>> _______________________________________________
> >>>>>> sbmlsolver-devel mailing list
> >>>>>> sbm...@li...
> >>>>>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
> >>>>>=20
> >>>>> _______________________________________________
> >>>>> sbmlsolver-devel mailing list
> >>>>> sbm...@li...
> >>>>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
> >>>>=20
> >>>> _______________________________________________
> >>>> sbmlsolver-devel mailing list
> >>>> sbm...@li...
> >>>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
> >>>=20
> >>> _______________________________________________
> >>> sbmlsolver-devel mailing list
> >>> sbm...@li...
> >>> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
> >>=20
> >>=20
> >> _______________________________________________
> >> sbmlsolver-devel mailing list
> >> sbm...@li...
> >> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
> >>=20
> >
>=20
>=20
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>=20

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

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

> However, this might still only be correct, if compartments are 
determined by
> assignment rules.
> If they are determined by rate rules, then the replacement S1 -> S1/V 
would
> also be required for constructing the jacobi matrix.

I just realized that it needs generally be done for constructing the 
jacobi because the jacoby should of course also be e.g. (dS1/dt)/dS2 and 
not (dS1/dt)/d[S2].

r


On Tue, 20 Jun 2006, Rainer Machne wrote:

> stefan
>
>
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>> 
>> yes.
>> and that's exactly the same what you are planning,
>> if you want to use conc ans subst in parallel.
>> only the point of conversion differs.
>
> not really. i wasn't planning any conversion of the math,
> except for:
>
> current ODE: d[S1]/dt = f(S1, S2, ...) / V
> planned ODE: d S1 /dt = f(S1, S2, ...)
>
> then we have two data->value arrays, let's say
>
> data->odeValues = S1, S2, S3, etc.
> and
> data->concValues = [S1], [S2], [S3] etc.
>
>
> The first array `odeValues' corresponds to the SUNDIALS values in NVector y.
> The second array `concValues' will be used in evaluateAST, for evaluation the 
> math.
>
>
> This would be equivalent of course, with just replacing all [S1] in all math 
> expressions with S1/volume, and only using one array data->odeValues
>
> While the latter approach might be easier to implement it has a probably the 
> not so small drawback, that for each species in each math expression the 
> evaluation needs one more operation, i.e. the division, while this division 
> could be done only once for each species whenever data->odeValues is updated 
> (i.e. e.g. in the SOSlib/SUNDIALS funciton f).
>
> However, this might still only be correct, if compartments are determined by 
> assignment rules.
> If they are determined by rate rules, then the replacement S1 -> S1/V would 
> also be required for constructing the jacobi matrix.
>
>> of course you have to check if the species in your kinetic law have 
>> dimension
>> conc or substance!
>> sbml allows the usage of both, doesn't it?
>
> See section 4.6.4 here
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-2.html#SECTION00046000000000000000
>
> "
> The units of the species are of the form substance/size units
>
> ...
>
> The units of the species are used in the following ways:
>
> ...
>    * The species identifier has these units when it appears as a numerical 
> quantity in a mathematical formula expressed in MathML (discussed in Section 
> 3.6.3).
> "
>
> So i guess, the units of species in math expressions are always in 
> concentrations!
>
> Rainer
>
>
>
>
>
>> 
>> i mean, every species is in a well-defined compartment, isn't it?
>> so you can use the formula [S]=S/V whenever you want...
>> 
>> i would eliminate conc as early as possile, but that's a matter of 
> taste.
>
>
>
> On Tue, 20 Jun 2006, Stefan Mueller wrote:
>
>> hi rainer!
>> 
>>> stefan
>>> 
>>>>> Currently the ODE would be:
>>>>> 
>>>>> d[S1]/dt = k1 * [S1] / V1
>>>> 
>>>> this is not valid for variable compartments.
>>> 
>>> i said that it's invalid for variable compartments some lines below:
>>>>> which is only valid for constant compartments V1 and V2.
>>>> 
>>>> please don't use it for comparison!
>>> 
>>> it's not used for comparison, but it is the current implementation in
>>> SOSlib!
>> 
>> i said this just to avoid any additional misunderstanding.
>> the whole topic seems to be quite confusing...
>> 
>>>>> At least for the second ODE, depending in two different volumes,
>>>>> the conversion to ODEs depending on substance instead of concentration
>>>>> is not as simple as
>>>>> 
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>> 
>>>> why dou you refer to this law for a reaction of order 2?
>>> 
>>> sorry, at first i wanted to make the second reaction in V2 of 2nd order,
>>> e.g.,
>>> 
>>> in V2:
>>> S2 + S3 -> P: k2*[S2]*[S3]
>>> 
>>>>> as your example before,
>>>>> but for S2 I get:
>>>>> 
>>>>> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>>>>> 
>>>>> Right?
>>>> 
>>>> no!
>>>> dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2
>>>> it's as simple as that!
>>>> (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
>>>> that's the advantage of kinetic laws (and consequently odes) for
>>>> substances!
>>> 
>>> of course! sorry, i got quite confused obviously.
>>> 
>>> and for the above second order reaction in V2 it would be:
>>> 
>>>   dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
>>> 
>>> Right this time?
>> 
>> yes :)
>> 
>>>>> Thus I would really have to check the type of reaction and edit a lot of
>>>>> kinetic laws before constructing the ODEs.
>>>> 
>>>> the above example shows that this is not necessary.
>>>> maybe there are other examples...
>>> 
>>> so you say, we could convert kinetic laws that don't explicitly contain
>>> the volume by simply differentiating these three cases and
>>> 
>>>>>> order 1:
>>>>>> S1->junk: k*[S1]
>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>> 
>>> divide k by V(0) for 1st order reactions
>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>> 
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>> 
>>> V^2
>>> 
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>> 
>> yes.
>> and that's exactly the same what you are planning,
>> if you want to use conc ans subst in parallel.
>> only the point of conversion differs.
>> 
>> i mean, every species is in a well-defined compartment, isn't it?
>> so you can use the formula [S]=S/V whenever you want...
>> 
>> i would eliminate conc as early as possile, but that's a matter of taste.
>> 
>>> but there are also enzymatic reactions, with in often quite complex
>>> kinetic laws. we would really have to do unit checking, recognizing the
>>> implemented reaction mechanisms etc.
>> 
>> of course you have to check if the species in your kinetic law have 
>> dimension
>> conc or substance!
>> sbml allows the usage of both, doesn't it?
>> 
>>> also, we don't really know in which compartment reactions happen, but
>>> could just try to get this information from the reactants' compartments.
>> 
>> this does not matter, for 2 reasons:
>> 
>> 1. if the kinetic law is specified correctly (including the compartment 
>> size),
>> e.g. A+B->C: k*V*[A]*[B],
>> the the compartment size is considered automatically.
>> 
>> 2. the conversion [S]=S/V does not depend on reactions, but only on species
>> (and their compartments).
>> 
>>> so i guess, it would be much easier to
>>> 
>>> a) use ODEs for substances but dependent on concentrations, here the only
>>> thing we need to do is to leave out the compartment division at the end of
>>> Species_odeFromReactions and have a double bookkeeping of amounts and
>>> concentrations
>> 
>> to my taste, double book keeping is dangerous and produces overhead.
>> 
>>> and
>>> 
>>> b) just ignore WRONG models that have a variable compartment but don't
>>> have the compartment explicitly in their kinetic laws
>>> (or check for such cases and issue an error);
>>> 
>>> isn't it an issue of correct model writing rather then correct model
>>> solving?
>> 
>> i agree.
>> 
>>> However, another problem might be that to calculate in amounts rather then
>>> in concentrations also changes the use of absolute errors and maybe
>>> also the general the performance of the integrator, i guess. ??
>> 
>> abs errors might change, but that's not an issue.
>> rel errors should stay the same.
>> (in one-compartment models everything is just multiplied by the volume.)
>> 
>>> maybe we will have to really implement both in parallel,
>> 
>> i hope not :)
>> 
>>> * ODEs for concentrations/time if no variable compartments are in the
>>> model
>>> * ODEs for ammounts/time if any compartment is variable
>>> 
>>> ?
>>> 
>>> Rainer
>> 
>> maybe we discuss the next round via telephone :)
>> 
>> cheers, s
>> 
>> 
>>>>> Thus I think it just doesn't work with kinetic laws that don't
>>>>> include the compartment explicitly!
>>>> 
>>>> i don't know exactly what you mean.
>>>> of course kinetic laws with explicit compartment are "nicer",
>>>> but you can always transform any kinetic law to a "nice" one...
>>>> 
>>>>> But if they include the compartment we can write ODEs for substances but
>>>>> have concentrations in the right hand side of the ODEs. We then just
>>>>> need to a second array for the concentrations, as discussed before.
>>>> 
>>>> if we still use concentrations is a separate problem,
>>>> independent of which we discussed until now.
>>>> it's a matter of taste, if we eliminate the concentrations  from the
>>>> kinetic laws and get a (possibly nonlinear) dependence on the
>>>> compartments OR if we keep concentrations and always have linear
>>>> dependence on the comp.
>>>> 
>>>>> Rainer
>>>> 
>>>> cheers,
>>>> stefan.
>>>> 
>>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
>>>>>> let me summarize:
>>>>>> kinetic laws have units substance/time,
>>>>>> and they depend on concentrations or substances.
>>>>>> 
>>>>>> we SHOULD write down odes for substances,
>>>>>> and we COULD eliminate concs completely.
>>>>>> (of course, we would have to convert initial concs to substances.)
>>>>>> 
>>>>>> in any case, odes for substances depend on the compartment size
>>>>>> (linearly or according to the order of the reaction):
>>>>>> 
>>>>>> order 1:
>>>>>> S1->junk: k*[S1]
>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>>>>> V^2
>>>>>> 
>>>>>> ...
>>>>>> 
>>>>>> if the kinetic laws included the compartment explicitly (which is
>>>>>> desirable), we would have the following situation:
>>>>>> 
>>>>>> order 1:
>>>>>> S1->junk: r*[S1]*V
>>>>>> dS1/dt =  r * [S1] * V= r * S1
>>>>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: r*[S1]*[S2]*V
>>>>>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>>>>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
>>>>>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>>>>>> 
>>>>>> as stated earlier, odes for substances have got another advantage. for
>>>>>> the compartment size we could use either a rate rule or an assignment
>>>>>> rule: dV/dt = f(t, V, ...)
>>>>>> or
>>>>>> V= f(t, V, ...)
>>>>>> 
>>>>>> cheers,
>>>>>> stefan.
>>>>>> 
>>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>>>>>>> andrew:
>>>>>>> 
>>>>>>> for the SBML system
>>>>>>> 	S1 in C
>>>>>>> 	S1->junk : k * [S1]
>>>>>>> 	dC/dt = 1
>>>>>>> 
>>>>>>> the right math is:
>>>>>>> 	dS1/dt = k * [S1] / C(0) * C
>>>>>>> 	dC/dt = 1
>>>>>>> with initial cond:
>>>>>>> 	S1(0) = [S1](0) * C(0)
>>>>>>> 	C(0)
>>>>>>> 
>>>>>>> notation:
>>>>>>> 	S1 ... substance
>>>>>>> 	[S1] ... conc
>>>>>>> 
>>>>>>> the "proposed solution" (who proposed this?) has the right units,
>>>>>>> but is wrong otherwise... :)
>>>>>>> 
>>>>>>> as rainer and me already discussed, everything would be much clearer,
>>>>>>> if the kinetic law included the compartment explicitly:
>>>>>>> 	S1 in C
>>>>>>> 	S1->junk : C * r * [S1]
>>>>>>> 	dC/dt = 1
>>>>>>> 
>>>>>>> with:
>>>>>>> 	r ... the "real" chemical rate constant.
>>>>>>> 
>>>>>>> rainer:
>>>>>>> 
>>>>>>> chemists always knew that rate constants are not constant with respect
>>>>>>> to pH, T, ... as a physicist i have to be fair! :)
>>>>>>> 
>>>>>>> to conclude: the math is not really tricky.
>>>>>>> odes for substances are better suited to handle both multiple and
>>>>>>> variable compartments. of course, we have to take care of many
>>>>>>> implementation details...
>>>>>>> 
>>>>>>> cheers,
>>>>>>> stefan.
>>>>>> 
>>>>>> _______________________________________________
>>>>>> 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
>> 
>

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

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

This

> current ODE: d[S1]/dt = f(S1, S2, ...) / V
> planned ODE: d S1 /dt = f(S1, S2, ...)

should of course be

current ODE: d[S1]/dt = f([S1], [S2], ...) / V
planned ODE: d S1 /dt = f([S1], [S2], ...)


On Tue, 20 Jun 2006, Rainer Machne wrote:

> stefan
>
>
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>> 
>> yes.
>> and that's exactly the same what you are planning,
>> if you want to use conc ans subst in parallel.
>> only the point of conversion differs.
>
> not really. i wasn't planning any conversion of the math,
> except for:
>
> current ODE: d[S1]/dt = f(S1, S2, ...) / V
> planned ODE: d S1 /dt = f(S1, S2, ...)
>
> then we have two data->value arrays, let's say
>
> data->odeValues = S1, S2, S3, etc.
> and
> data->concValues = [S1], [S2], [S3] etc.
>
>
> The first array `odeValues' corresponds to the SUNDIALS values in NVector y.
> The second array `concValues' will be used in evaluateAST, for evaluation the 
> math.
>
>
> This would be equivalent of course, with just replacing all [S1] in all math 
> expressions with S1/volume, and only using one array data->odeValues
>
> While the latter approach might be easier to implement it has a probably the 
> not so small drawback, that for each species in each math expression the 
> evaluation needs one more operation, i.e. the division, while this division 
> could be done only once for each species whenever data->odeValues is updated 
> (i.e. e.g. in the SOSlib/SUNDIALS funciton f).
>
> However, this might still only be correct, if compartments are determined by 
> assignment rules.
> If they are determined by rate rules, then the replacement S1 -> S1/V would 
> also be required for constructing the jacobi matrix.
>
>> of course you have to check if the species in your kinetic law have 
>> dimension
>> conc or substance!
>> sbml allows the usage of both, doesn't it?
>
> See section 4.6.4 here
> http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-2.html#SECTION00046000000000000000
>
> "
> The units of the species are of the form substance/size units
>
> ...
>
> The units of the species are used in the following ways:
>
> ...
>    * The species identifier has these units when it appears as a numerical 
> quantity in a mathematical formula expressed in MathML (discussed in Section 
> 3.6.3).
> "
>
> So i guess, the units of species in math expressions are always in 
> concentrations!
>
> Rainer
>
>
>
>
>
>> 
>> i mean, every species is in a well-defined compartment, isn't it?
>> so you can use the formula [S]=S/V whenever you want...
>> 
>> i would eliminate conc as early as possile, but that's a matter of 
> taste.
>
>
>
> On Tue, 20 Jun 2006, Stefan Mueller wrote:
>
>> hi rainer!
>> 
>>> stefan
>>> 
>>>>> Currently the ODE would be:
>>>>> 
>>>>> d[S1]/dt = k1 * [S1] / V1
>>>> 
>>>> this is not valid for variable compartments.
>>> 
>>> i said that it's invalid for variable compartments some lines below:
>>>>> which is only valid for constant compartments V1 and V2.
>>>> 
>>>> please don't use it for comparison!
>>> 
>>> it's not used for comparison, but it is the current implementation in
>>> SOSlib!
>> 
>> i said this just to avoid any additional misunderstanding.
>> the whole topic seems to be quite confusing...
>> 
>>>>> At least for the second ODE, depending in two different volumes,
>>>>> the conversion to ODEs depending on substance instead of concentration
>>>>> is not as simple as
>>>>> 
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>> 
>>>> why dou you refer to this law for a reaction of order 2?
>>> 
>>> sorry, at first i wanted to make the second reaction in V2 of 2nd order,
>>> e.g.,
>>> 
>>> in V2:
>>> S2 + S3 -> P: k2*[S2]*[S3]
>>> 
>>>>> as your example before,
>>>>> but for S2 I get:
>>>>> 
>>>>> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>>>>> 
>>>>> Right?
>>>> 
>>>> no!
>>>> dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2
>>>> it's as simple as that!
>>>> (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
>>>> that's the advantage of kinetic laws (and consequently odes) for
>>>> substances!
>>> 
>>> of course! sorry, i got quite confused obviously.
>>> 
>>> and for the above second order reaction in V2 it would be:
>>> 
>>>   dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
>>> 
>>> Right this time?
>> 
>> yes :)
>> 
>>>>> Thus I would really have to check the type of reaction and edit a lot of
>>>>> kinetic laws before constructing the ODEs.
>>>> 
>>>> the above example shows that this is not necessary.
>>>> maybe there are other examples...
>>> 
>>> so you say, we could convert kinetic laws that don't explicitly contain
>>> the volume by simply differentiating these three cases and
>>> 
>>>>>> order 1:
>>>>>> S1->junk: k*[S1]
>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>> 
>>> divide k by V(0) for 1st order reactions
>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>> 
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>> 
>>> V^2
>>> 
>>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>>> reactions
>> 
>> yes.
>> and that's exactly the same what you are planning,
>> if you want to use conc ans subst in parallel.
>> only the point of conversion differs.
>> 
>> i mean, every species is in a well-defined compartment, isn't it?
>> so you can use the formula [S]=S/V whenever you want...
>> 
>> i would eliminate conc as early as possile, but that's a matter of taste.
>> 
>>> but there are also enzymatic reactions, with in often quite complex
>>> kinetic laws. we would really have to do unit checking, recognizing the
>>> implemented reaction mechanisms etc.
>> 
>> of course you have to check if the species in your kinetic law have 
>> dimension
>> conc or substance!
>> sbml allows the usage of both, doesn't it?
>> 
>>> also, we don't really know in which compartment reactions happen, but
>>> could just try to get this information from the reactants' compartments.
>> 
>> this does not matter, for 2 reasons:
>> 
>> 1. if the kinetic law is specified correctly (including the compartment 
>> size),
>> e.g. A+B->C: k*V*[A]*[B],
>> the the compartment size is considered automatically.
>> 
>> 2. the conversion [S]=S/V does not depend on reactions, but only on species
>> (and their compartments).
>> 
>>> so i guess, it would be much easier to
>>> 
>>> a) use ODEs for substances but dependent on concentrations, here the only
>>> thing we need to do is to leave out the compartment division at the end of
>>> Species_odeFromReactions and have a double bookkeeping of amounts and
>>> concentrations
>> 
>> to my taste, double book keeping is dangerous and produces overhead.
>> 
>>> and
>>> 
>>> b) just ignore WRONG models that have a variable compartment but don't
>>> have the compartment explicitly in their kinetic laws
>>> (or check for such cases and issue an error);
>>> 
>>> isn't it an issue of correct model writing rather then correct model
>>> solving?
>> 
>> i agree.
>> 
>>> However, another problem might be that to calculate in amounts rather then
>>> in concentrations also changes the use of absolute errors and maybe
>>> also the general the performance of the integrator, i guess. ??
>> 
>> abs errors might change, but that's not an issue.
>> rel errors should stay the same.
>> (in one-compartment models everything is just multiplied by the volume.)
>> 
>>> maybe we will have to really implement both in parallel,
>> 
>> i hope not :)
>> 
>>> * ODEs for concentrations/time if no variable compartments are in the
>>> model
>>> * ODEs for ammounts/time if any compartment is variable
>>> 
>>> ?
>>> 
>>> Rainer
>> 
>> maybe we discuss the next round via telephone :)
>> 
>> cheers, s
>> 
>> 
>>>>> Thus I think it just doesn't work with kinetic laws that don't
>>>>> include the compartment explicitly!
>>>> 
>>>> i don't know exactly what you mean.
>>>> of course kinetic laws with explicit compartment are "nicer",
>>>> but you can always transform any kinetic law to a "nice" one...
>>>> 
>>>>> But if they include the compartment we can write ODEs for substances but
>>>>> have concentrations in the right hand side of the ODEs. We then just
>>>>> need to a second array for the concentrations, as discussed before.
>>>> 
>>>> if we still use concentrations is a separate problem,
>>>> independent of which we discussed until now.
>>>> it's a matter of taste, if we eliminate the concentrations  from the
>>>> kinetic laws and get a (possibly nonlinear) dependence on the
>>>> compartments OR if we keep concentrations and always have linear
>>>> dependence on the comp.
>>>> 
>>>>> Rainer
>>>> 
>>>> cheers,
>>>> stefan.
>>>> 
>>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
>>>>>> let me summarize:
>>>>>> kinetic laws have units substance/time,
>>>>>> and they depend on concentrations or substances.
>>>>>> 
>>>>>> we SHOULD write down odes for substances,
>>>>>> and we COULD eliminate concs completely.
>>>>>> (of course, we would have to convert initial concs to substances.)
>>>>>> 
>>>>>> in any case, odes for substances depend on the compartment size
>>>>>> (linearly or according to the order of the reaction):
>>>>>> 
>>>>>> order 1:
>>>>>> S1->junk: k*[S1]
>>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>>>>> V^2
>>>>>> 
>>>>>> ...
>>>>>> 
>>>>>> if the kinetic laws included the compartment explicitly (which is
>>>>>> desirable), we would have the following situation:
>>>>>> 
>>>>>> order 1:
>>>>>> S1->junk: r*[S1]*V
>>>>>> dS1/dt =  r * [S1] * V= r * S1
>>>>>> 
>>>>>> order 2:
>>>>>> S1+S2->bla: r*[S1]*[S2]*V
>>>>>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>>>>>> 
>>>>>> order 3:
>>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
>>>>>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>>>>>> 
>>>>>> as stated earlier, odes for substances have got another advantage. for
>>>>>> the compartment size we could use either a rate rule or an assignment
>>>>>> rule: dV/dt = f(t, V, ...)
>>>>>> or
>>>>>> V= f(t, V, ...)
>>>>>> 
>>>>>> cheers,
>>>>>> stefan.
>>>>>> 
>>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>>>>>>> andrew:
>>>>>>> 
>>>>>>> for the SBML system
>>>>>>> 	S1 in C
>>>>>>> 	S1->junk : k * [S1]
>>>>>>> 	dC/dt = 1
>>>>>>> 
>>>>>>> the right math is:
>>>>>>> 	dS1/dt = k * [S1] / C(0) * C
>>>>>>> 	dC/dt = 1
>>>>>>> with initial cond:
>>>>>>> 	S1(0) = [S1](0) * C(0)
>>>>>>> 	C(0)
>>>>>>> 
>>>>>>> notation:
>>>>>>> 	S1 ... substance
>>>>>>> 	[S1] ... conc
>>>>>>> 
>>>>>>> the "proposed solution" (who proposed this?) has the right units,
>>>>>>> but is wrong otherwise... :)
>>>>>>> 
>>>>>>> as rainer and me already discussed, everything would be much clearer,
>>>>>>> if the kinetic law included the compartment explicitly:
>>>>>>> 	S1 in C
>>>>>>> 	S1->junk : C * r * [S1]
>>>>>>> 	dC/dt = 1
>>>>>>> 
>>>>>>> with:
>>>>>>> 	r ... the "real" chemical rate constant.
>>>>>>> 
>>>>>>> rainer:
>>>>>>> 
>>>>>>> chemists always knew that rate constants are not constant with respect
>>>>>>> to pH, T, ... as a physicist i have to be fair! :)
>>>>>>> 
>>>>>>> to conclude: the math is not really tricky.
>>>>>>> odes for substances are better suited to handle both multiple and
>>>>>>> variable compartments. of course, we have to take care of many
>>>>>>> implementation details...
>>>>>>> 
>>>>>>> cheers,
>>>>>>> stefan.
>>>>>> 
>>>>>> _______________________________________________
>>>>>> 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
>> 
>

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

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

stefan


>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>> reactions
>
> yes.
> and that's exactly the same what you are planning,
> if you want to use conc ans subst in parallel.
> only the point of conversion differs.

not really. i wasn't planning any conversion of the math,
  except for:

current ODE: d[S1]/dt = f(S1, S2, ...) / V
planned ODE: d S1 /dt = f(S1, S2, ...)

then we have two data->value arrays, let's say

data->odeValues = S1, S2, S3, etc.
and
data->concValues = [S1], [S2], [S3] etc.


The first array `odeValues' corresponds to the SUNDIALS values in NVector 
y.
The second array `concValues' will be used in evaluateAST, for evaluation 
the math.


This would be equivalent of course, with just replacing all [S1] in all 
math expressions with S1/volume, and only using one array data->odeValues

While the latter approach might be easier to implement it has a probably 
the not so small drawback, that for each species in each math expression 
the evaluation needs one more operation, i.e. the division, while this 
division could be done only once for each species whenever data->odeValues 
is updated (i.e. e.g. in the SOSlib/SUNDIALS funciton f).

However, this might still only be correct, if compartments are determined 
by assignment rules.
If they are determined by rate rules, then the replacement S1 -> S1/V 
would also be required for constructing the jacobi matrix.

> of course you have to check if the species in your kinetic law have dimension
> conc or substance!
> sbml allows the usage of both, doesn't it?

See section 4.6.4 here
http://sbml.org/specifications/sbml-level-2/version-1/html/sbml-level-2.html#SECTION00046000000000000000

"
The units of the species are of the form substance/size units

...

The units of the species are used in the following ways:

...
     * The species identifier has these units when it appears as a 
numerical quantity in a mathematical formula expressed in MathML 
(discussed in Section 3.6.3).
"

So i guess, the units of species in math expressions are always in 
concentrations!

Rainer





>
> i mean, every species is in a well-defined compartment, isn't it?
> so you can use the formula [S]=S/V whenever you want...
>
> i would eliminate conc as early as possile, but that's a matter of 
taste.



On Tue, 20 Jun 2006, Stefan Mueller wrote:

> hi rainer!
>
>> stefan
>>
>>>> Currently the ODE would be:
>>>>
>>>> d[S1]/dt = k1 * [S1] / V1
>>>
>>> this is not valid for variable compartments.
>>
>> i said that it's invalid for variable compartments some lines below:
>>>> which is only valid for constant compartments V1 and V2.
>>>
>>> please don't use it for comparison!
>>
>> it's not used for comparison, but it is the current implementation in
>> SOSlib!
>
> i said this just to avoid any additional misunderstanding.
> the whole topic seems to be quite confusing...
>
>>>> At least for the second ODE, depending in two different volumes,
>>>> the conversion to ODEs depending on substance instead of concentration
>>>> is not as simple as
>>>>
>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>
>>> why dou you refer to this law for a reaction of order 2?
>>
>> sorry, at first i wanted to make the second reaction in V2 of 2nd order,
>> e.g.,
>>
>> in V2:
>> S2 + S3 -> P: k2*[S2]*[S3]
>>
>>>> as your example before,
>>>> but for S2 I get:
>>>>
>>>> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>>>>
>>>> Right?
>>>
>>> no!
>>> dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2
>>> it's as simple as that!
>>> (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
>>> that's the advantage of kinetic laws (and consequently odes) for
>>> substances!
>>
>> of course! sorry, i got quite confused obviously.
>>
>> and for the above second order reaction in V2 it would be:
>>
>>   dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
>>
>> Right this time?
>
> yes :)
>
>>>> Thus I would really have to check the type of reaction and edit a lot of
>>>> kinetic laws before constructing the ODEs.
>>>
>>> the above example shows that this is not necessary.
>>> maybe there are other examples...
>>
>> so you say, we could convert kinetic laws that don't explicitly contain
>> the volume by simply differentiating these three cases and
>>
>>>>> order 1:
>>>>> S1->junk: k*[S1]
>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>
>> divide k by V(0) for 1st order reactions
>>
>>>>> order 2:
>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>
>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>> reactions
>>
>>>>> order 3:
>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>
>> V^2
>>
>> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
>> reactions
>
> yes.
> and that's exactly the same what you are planning,
> if you want to use conc ans subst in parallel.
> only the point of conversion differs.
>
> i mean, every species is in a well-defined compartment, isn't it?
> so you can use the formula [S]=S/V whenever you want...
>
> i would eliminate conc as early as possile, but that's a matter of taste.
>
>> but there are also enzymatic reactions, with in often quite complex
>> kinetic laws. we would really have to do unit checking, recognizing the
>> implemented reaction mechanisms etc.
>
> of course you have to check if the species in your kinetic law have dimension
> conc or substance!
> sbml allows the usage of both, doesn't it?
>
>> also, we don't really know in which compartment reactions happen, but
>> could just try to get this information from the reactants' compartments.
>
> this does not matter, for 2 reasons:
>
> 1. if the kinetic law is specified correctly (including the compartment size),
> e.g. A+B->C: k*V*[A]*[B],
> the the compartment size is considered automatically.
>
> 2. the conversion [S]=S/V does not depend on reactions, but only on species
> (and their compartments).
>
>> so i guess, it would be much easier to
>>
>> a) use ODEs for substances but dependent on concentrations, here the only
>> thing we need to do is to leave out the compartment division at the end of
>> Species_odeFromReactions and have a double bookkeeping of amounts and
>> concentrations
>
> to my taste, double book keeping is dangerous and produces overhead.
>
>> and
>>
>> b) just ignore WRONG models that have a variable compartment but don't
>> have the compartment explicitly in their kinetic laws
>> (or check for such cases and issue an error);
>>
>> isn't it an issue of correct model writing rather then correct model
>> solving?
>
> i agree.
>
>> However, another problem might be that to calculate in amounts rather then
>> in concentrations also changes the use of absolute errors and maybe
>> also the general the performance of the integrator, i guess. ??
>
> abs errors might change, but that's not an issue.
> rel errors should stay the same.
> (in one-compartment models everything is just multiplied by the volume.)
>
>> maybe we will have to really implement both in parallel,
>
> i hope not :)
>
>> * ODEs for concentrations/time if no variable compartments are in the
>> model
>> * ODEs for ammounts/time if any compartment is variable
>>
>> ?
>>
>> Rainer
>
> maybe we discuss the next round via telephone :)
>
> cheers, s
>
>
>>>> Thus I think it just doesn't work with kinetic laws that don't
>>>> include the compartment explicitly!
>>>
>>> i don't know exactly what you mean.
>>> of course kinetic laws with explicit compartment are "nicer",
>>> but you can always transform any kinetic law to a "nice" one...
>>>
>>>> But if they include the compartment we can write ODEs for substances but
>>>> have concentrations in the right hand side of the ODEs. We then just
>>>> need to a second array for the concentrations, as discussed before.
>>>
>>> if we still use concentrations is a separate problem,
>>> independent of which we discussed until now.
>>> it's a matter of taste, if we eliminate the concentrations  from the
>>> kinetic laws and get a (possibly nonlinear) dependence on the
>>> compartments OR if we keep concentrations and always have linear
>>> dependence on the comp.
>>>
>>>> Rainer
>>>
>>> cheers,
>>> stefan.
>>>
>>>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
>>>>> let me summarize:
>>>>> kinetic laws have units substance/time,
>>>>> and they depend on concentrations or substances.
>>>>>
>>>>> we SHOULD write down odes for substances,
>>>>> and we COULD eliminate concs completely.
>>>>> (of course, we would have to convert initial concs to substances.)
>>>>>
>>>>> in any case, odes for substances depend on the compartment size
>>>>> (linearly or according to the order of the reaction):
>>>>>
>>>>> order 1:
>>>>> S1->junk: k*[S1]
>>>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>>>
>>>>> order 2:
>>>>> S1+S2->bla: k*[S1]*[S2]
>>>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>>>
>>>>> order 3:
>>>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>>>>> V^2
>>>>>
>>>>> ...
>>>>>
>>>>> if the kinetic laws included the compartment explicitly (which is
>>>>> desirable), we would have the following situation:
>>>>>
>>>>> order 1:
>>>>> S1->junk: r*[S1]*V
>>>>> dS1/dt =  r * [S1] * V= r * S1
>>>>>
>>>>> order 2:
>>>>> S1+S2->bla: r*[S1]*[S2]*V
>>>>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>>>>>
>>>>> order 3:
>>>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
>>>>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>>>>>
>>>>> as stated earlier, odes for substances have got another advantage. for
>>>>> the compartment size we could use either a rate rule or an assignment
>>>>> rule: dV/dt = f(t, V, ...)
>>>>> or
>>>>> V= f(t, V, ...)
>>>>>
>>>>> cheers,
>>>>> stefan.
>>>>>
>>>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>>>>>> andrew:
>>>>>>
>>>>>> for the SBML system
>>>>>> 	S1 in C
>>>>>> 	S1->junk : k * [S1]
>>>>>> 	dC/dt = 1
>>>>>>
>>>>>> the right math is:
>>>>>> 	dS1/dt = k * [S1] / C(0) * C
>>>>>> 	dC/dt = 1
>>>>>> with initial cond:
>>>>>> 	S1(0) = [S1](0) * C(0)
>>>>>> 	C(0)
>>>>>>
>>>>>> notation:
>>>>>> 	S1 ... substance
>>>>>> 	[S1] ... conc
>>>>>>
>>>>>> the "proposed solution" (who proposed this?) has the right units,
>>>>>> but is wrong otherwise... :)
>>>>>>
>>>>>> as rainer and me already discussed, everything would be much clearer,
>>>>>> if the kinetic law included the compartment explicitly:
>>>>>> 	S1 in C
>>>>>> 	S1->junk : C * r * [S1]
>>>>>> 	dC/dt = 1
>>>>>>
>>>>>> with:
>>>>>> 	r ... the "real" chemical rate constant.
>>>>>>
>>>>>> rainer:
>>>>>>
>>>>>> chemists always knew that rate constants are not constant with respect
>>>>>> to pH, T, ... as a physicist i have to be fair! :)
>>>>>>
>>>>>> to conclude: the math is not really tricky.
>>>>>> odes for substances are better suited to handle both multiple and
>>>>>> variable compartments. of course, we have to take care of many
>>>>>> implementation details...
>>>>>>
>>>>>> cheers,
>>>>>> stefan.
>>>>>
>>>>> _______________________________________________
>>>>> 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
>

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

[Stefan M. <ste...@oe...>]

hi rainer!

> stefan
>
> >> Currently the ODE would be:
> >>
> >> d[S1]/dt = k1 * [S1] / V1
> >
> > this is not valid for variable compartments.
>
> i said that it's invalid for variable compartments some lines below:
> >> which is only valid for constant compartments V1 and V2.
> >
> > please don't use it for comparison!
>
> it's not used for comparison, but it is the current implementation in
> SOSlib!

i said this just to avoid any additional misunderstanding.
the whole topic seems to be quite confusing...

> >> At least for the second ODE, depending in two different volumes,
> >> the conversion to ODEs depending on substance instead of concentration
> >> is not as simple as
> >>
> >>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
> >
> > why dou you refer to this law for a reaction of order 2?
>
> sorry, at first i wanted to make the second reaction in V2 of 2nd order,
> e.g.,
>
> in V2:
> S2 + S3 -> P: k2*[S2]*[S3]
>
> >> as your example before,
> >> but for S2 I get:
> >>
> >> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
> >>
> >> Right?
> >
> > no!
> > dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2
> > it's as simple as that!
> > (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
> > that's the advantage of kinetic laws (and consequently odes) for
> > substances!
>
> of course! sorry, i got quite confused obviously.
>
> and for the above second order reaction in V2 it would be:
>
>   dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2
>
> Right this time?

yes :)

> >> Thus I would really have to check the type of reaction and edit a lot of
> >> kinetic laws before constructing the ODEs.
> >
> > the above example shows that this is not necessary.
> > maybe there are other examples...
>
> so you say, we could convert kinetic laws that don't explicitly contain
> the volume by simply differentiating these three cases and
>
> >>> order 1:
> >>> S1->junk: k*[S1]
> >>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>
> divide k by V(0) for 1st order reactions
>
> >>> order 2:
> >>> S1+S2->bla: k*[S1]*[S2]
> >>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>
> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
> reactions
>
> >>> order 3:
> >>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
>
> V^2
>
> divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
> reactions

yes.
and that's exactly the same what you are planning,
if you want to use conc ans subst in parallel.
only the point of conversion differs.

i mean, every species is in a well-defined compartment, isn't it?
so you can use the formula [S]=S/V whenever you want...

i would eliminate conc as early as possile, but that's a matter of taste.

> but there are also enzymatic reactions, with in often quite complex
> kinetic laws. we would really have to do unit checking, recognizing the
> implemented reaction mechanisms etc.

of course you have to check if the species in your kinetic law have dimension 
conc or substance!
sbml allows the usage of both, doesn't it?

> also, we don't really know in which compartment reactions happen, but
> could just try to get this information from the reactants' compartments.

this does not matter, for 2 reasons:

1. if the kinetic law is specified correctly (including the compartment size),
e.g. A+B->C: k*V*[A]*[B],
the the compartment size is considered automatically.

2. the conversion [S]=S/V does not depend on reactions, but only on species
(and their compartments).

> so i guess, it would be much easier to
>
> a) use ODEs for substances but dependent on concentrations, here the only
> thing we need to do is to leave out the compartment division at the end of
> Species_odeFromReactions and have a double bookkeeping of amounts and
> concentrations

to my taste, double book keeping is dangerous and produces overhead.

> and
>
> b) just ignore WRONG models that have a variable compartment but don't
> have the compartment explicitly in their kinetic laws
> (or check for such cases and issue an error);
>
> isn't it an issue of correct model writing rather then correct model
> solving?

i agree.

> However, another problem might be that to calculate in amounts rather then
> in concentrations also changes the use of absolute errors and maybe
> also the general the performance of the integrator, i guess. ??

abs errors might change, but that's not an issue.
rel errors should stay the same.
(in one-compartment models everything is just multiplied by the volume.)

> maybe we will have to really implement both in parallel,

i hope not :)

> * ODEs for concentrations/time if no variable compartments are in the
> model
> * ODEs for ammounts/time if any compartment is variable
>
> ?
>
> Rainer

maybe we discuss the next round via telephone :)

cheers, s


> >> Thus I think it just doesn't work with kinetic laws that don't
> >> include the compartment explicitly!
> >
> > i don't know exactly what you mean.
> > of course kinetic laws with explicit compartment are "nicer",
> > but you can always transform any kinetic law to a "nice" one...
> >
> >> But if they include the compartment we can write ODEs for substances but
> >> have concentrations in the right hand side of the ODEs. We then just
> >> need to a second array for the concentrations, as discussed before.
> >
> > if we still use concentrations is a separate problem,
> > independent of which we discussed until now.
> > it's a matter of taste, if we eliminate the concentrations  from the
> > kinetic laws and get a (possibly nonlinear) dependence on the
> > compartments OR if we keep concentrations and always have linear
> > dependence on the comp.
> >
> >> Rainer
> >
> > cheers,
> > stefan.
> >
> >> On Mon, 19 Jun 2006, Stefan Mueller wrote:
> >>> let me summarize:
> >>> kinetic laws have units substance/time,
> >>> and they depend on concentrations or substances.
> >>>
> >>> we SHOULD write down odes for substances,
> >>> and we COULD eliminate concs completely.
> >>> (of course, we would have to convert initial concs to substances.)
> >>>
> >>> in any case, odes for substances depend on the compartment size
> >>> (linearly or according to the order of the reaction):
> >>>
> >>> order 1:
> >>> S1->junk: k*[S1]
> >>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
> >>>
> >>> order 2:
> >>> S1+S2->bla: k*[S1]*[S2]
> >>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
> >>>
> >>> order 3:
> >>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> >>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 /
> >>> V^2
> >>>
> >>> ...
> >>>
> >>> if the kinetic laws included the compartment explicitly (which is
> >>> desirable), we would have the following situation:
> >>>
> >>> order 1:
> >>> S1->junk: r*[S1]*V
> >>> dS1/dt =  r * [S1] * V= r * S1
> >>>
> >>> order 2:
> >>> S1+S2->bla: r*[S1]*[S2]*V
> >>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
> >>>
> >>> order 3:
> >>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> >>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
> >>>
> >>> as stated earlier, odes for substances have got another advantage. for
> >>> the compartment size we could use either a rate rule or an assignment
> >>> rule: dV/dt = f(t, V, ...)
> >>> or
> >>> V= f(t, V, ...)
> >>>
> >>> cheers,
> >>> stefan.
> >>>
> >>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> >>>> andrew:
> >>>>
> >>>> for the SBML system
> >>>> 	S1 in C
> >>>> 	S1->junk : k * [S1]
> >>>> 	dC/dt = 1
> >>>>
> >>>> the right math is:
> >>>> 	dS1/dt = k * [S1] / C(0) * C
> >>>> 	dC/dt = 1
> >>>> with initial cond:
> >>>> 	S1(0) = [S1](0) * C(0)
> >>>> 	C(0)
> >>>>
> >>>> notation:
> >>>> 	S1 ... substance
> >>>> 	[S1] ... conc
> >>>>
> >>>> the "proposed solution" (who proposed this?) has the right units,
> >>>> but is wrong otherwise... :)
> >>>>
> >>>> as rainer and me already discussed, everything would be much clearer,
> >>>> if the kinetic law included the compartment explicitly:
> >>>> 	S1 in C
> >>>> 	S1->junk : C * r * [S1]
> >>>> 	dC/dt = 1
> >>>>
> >>>> with:
> >>>> 	r ... the "real" chemical rate constant.
> >>>>
> >>>> rainer:
> >>>>
> >>>> chemists always knew that rate constants are not constant with respect
> >>>> to pH, T, ... as a physicist i have to be fair! :)
> >>>>
> >>>> to conclude: the math is not really tricky.
> >>>> odes for substances are better suited to handle both multiple and
> >>>> variable compartments. of course, we have to take care of many
> >>>> implementation details...
> >>>>
> >>>> cheers,
> >>>> stefan.
> >>>
> >>> _______________________________________________
> >>> 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

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

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

stefan




>> Currently the ODE would be:
>>
>> d[S1]/dt = k1 * [S1] / V1
>
> this is not valid for variable compartments.

i said that it's invalid for variable compartments some lines below:

>> which is only valid for constant compartments V1 and V2.

> please don't use it for comparison!

it's not used for comparison, but it is the current implementation in 
SOSlib!


>> At least for the second ODE, depending in two different volumes,
>> the conversion to ODEs depending on substance instead of concentration is
>> not as simple as
>>
>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>
> why dou you refer to this law for a reaction of order 2?
>

sorry, at first i wanted to make the second reaction in V2 of 2nd order, 
e.g.,

in V2:
S2 + S3 -> P: k2*[S2]*[S3]


>> as your example before,
>> but for S2 I get:
>>
>> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>>
>> Right?
>
> no!
> dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2
> it's as simple as that!
> (since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
> that's the advantage of kinetic laws (and consequently odes) for substances!

of course! sorry, i got quite confused obviously.

and for the above second order reaction in V2 it would be:

  dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S1 * S3 / V2

Right this time?



>
>> Thus I would really have to check the type of reaction and edit a lot of
>> kinetic laws before constructing the ODEs.
>
> the above example shows that this is not necessary.
> maybe there are other examples...


so you say, we could convert kinetic laws that don't explicitly contain 
the volume by simply differentiating these three cases and


>>> order 1:
>>> S1->junk: k*[S1]
>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>

divide k by V(0) for 1st order reactions


>>> order 2:
>>> S1+S2->bla: k*[S1]*[S2]
>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V

divide k by V(0) and the whole kinetic law again by V^2 for 2nd order 
reactions


>>>
>>> order 3:
>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 / 
V^2

divide k by V(0) and the whole kinetic law again by V^2 for 2nd order
reactions


but there are also enzymatic reactions, with in often quite complex 
kinetic laws. we would really have to do unit checking, recognizing the 
implemented reaction mechanisms etc.

also, we don't really know in which compartment reactions happen, but 
could just try to get this information from the reactants' compartments.

so i guess, it would be much easier to

a) use ODEs for substances but dependent on concentrations, here the only 
thing we need to do is to leave out the compartment division at the end of 
Species_odeFromReactions and have a double bookkeeping of amounts and 
concentrations

and

b) just ignore WRONG models that have a variable compartment but don't 
have the compartment explicitly in their kinetic laws
(or check for such cases and issue an error);

isn't it an issue of correct model writing rather then correct model 
solving?


However, another problem might be that to calculate in amounts rather then 
in concentrations also changes the use of absolute errors and maybe 
also the general the performance of the integrator, i guess. ??


maybe we will have to really implement both in parallel,
* ODEs for concentrations/time if no variable compartments are in the 
model 
* ODEs for ammounts/time if any compartment is variable

?

Rainer



>> Thus I think it just doesn't work with kinetic laws that don't
>> include the compartment explicitly!
>
> i don't know exactly what you mean.
> of course kinetic laws with explicit compartment are "nicer",
> but you can always transform any kinetic law to a "nice" one...
>
>> But if they include the compartment we can write ODEs for substances but
>> have concentrations in the right hand side of the ODEs. We then just need
>> to a second array for the concentrations, as discussed before.
>
> if we still use concentrations is a separate problem,
> independent of which we discussed until now.
> it's a matter of taste, if we eliminate the concentrations  from the kinetic
> laws and get a (possibly nonlinear) dependence on the compartments OR
> if we keep concentrations and always have linear dependence on the comp.
>
>> Rainer
>
> cheers,
> stefan.
>
>
>> On Mon, 19 Jun 2006, Stefan Mueller wrote:
>>> let me summarize:
>>> kinetic laws have units substance/time,
>>> and they depend on concentrations or substances.
>>>
>>> we SHOULD write down odes for substances,
>>> and we COULD eliminate concs completely.
>>> (of course, we would have to convert initial concs to substances.)
>>>
>>> in any case, odes for substances depend on the compartment size
>>> (linearly or according to the order of the reaction):
>>>
>>> order 1:
>>> S1->junk: k*[S1]
>>> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>>>
>>> order 2:
>>> S1+S2->bla: k*[S1]*[S2]
>>> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>>>
>>> order 3:
>>> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
>>> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 / V^2
>>>
>>> ...
>>>
>>> if the kinetic laws included the compartment explicitly (which is
>>> desirable), we would have the following situation:
>>>
>>> order 1:
>>> S1->junk: r*[S1]*V
>>> dS1/dt =  r * [S1] * V= r * S1
>>>
>>> order 2:
>>> S1+S2->bla: r*[S1]*[S2]*V
>>> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>>>
>>> order 3:
>>> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
>>> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>>>
>>> as stated earlier, odes for substances have got another advantage. for
>>> the compartment size we could use either a rate rule or an assignment
>>> rule: dV/dt = f(t, V, ...)
>>> or
>>> V= f(t, V, ...)
>>>
>>> cheers,
>>> stefan.
>>>
>>> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>>>> andrew:
>>>>
>>>> for the SBML system
>>>> 	S1 in C
>>>> 	S1->junk : k * [S1]
>>>> 	dC/dt = 1
>>>>
>>>> the right math is:
>>>> 	dS1/dt = k * [S1] / C(0) * C
>>>> 	dC/dt = 1
>>>> with initial cond:
>>>> 	S1(0) = [S1](0) * C(0)
>>>> 	C(0)
>>>>
>>>> notation:
>>>> 	S1 ... substance
>>>> 	[S1] ... conc
>>>>
>>>> the "proposed solution" (who proposed this?) has the right units,
>>>> but is wrong otherwise... :)
>>>>
>>>> as rainer and me already discussed, everything would be much clearer,
>>>> if the kinetic law included the compartment explicitly:
>>>> 	S1 in C
>>>> 	S1->junk : C * r * [S1]
>>>> 	dC/dt = 1
>>>>
>>>> with:
>>>> 	r ... the "real" chemical rate constant.
>>>>
>>>> rainer:
>>>>
>>>> chemists always knew that rate constants are not constant with respect
>>>> to pH, T, ... as a physicist i have to be fair! :)
>>>>
>>>> to conclude: the math is not really tricky.
>>>> odes for substances are better suited to handle both multiple and
>>>> variable compartments. of course, we have to take care of many
>>>> implementation details...
>>>>
>>>> cheers,
>>>> stefan.
>>>
>>> _______________________________________________
>>> 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
>

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

[Stefan M. <ste...@oe...>]

rainer

> Stefan
>
> i guess it might be different if we have to consider multiple
> compartments
>
> S1 is in V1;
> S2 and S3 are in V2;
>
> S1->S2: k1*[S1]

well, this is not a chemical reaction, but a transport process,
but let's see...

> S2->S3: k2*[S2]
>
> [S1] = S1 / V1
> [S2] = S2 / V2
> [S3] = S3 / V3
>
> Thus the two parameters include two different initial volumes, i.e.:
> k1 [volume/sec] = r1 * V1(0)
> k2 [volume/sec] = r2 * V2(0)
>
> Currently the ODE would be:
>
> d[S1]/dt = k1 * [S1] / V1

this is not valid for variable compartments.
please don't use it for comparison!
btw, this ode is missing a minus :)

> d[S2]/dt = ( k1*[S1] - k2*[S2] ) / V2
> d[S3]/dt  k2*[S2]*[S3] / V3
>
> which is only valid for constant compartments V1 and V2.
>
> At least for the second ODE, depending in two different volumes,
> the conversion to ODEs depending on substance instead of concentration is
> not as simple as
>
> > dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V

why dou you refer to this law for a reaction of order 2?

> (all the ODEs in your email are missing a minus, BTW!)
>
> It should probably be something like:
>
> dS1/dt = - k1 * [S1] / V1(0) * V1 = - k1/V1(0) * S1

yes, this is the right one!
(of order 1)

> as your example before,
> but for S2 I get:
>
> dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2
>
> Right?

no!
dS2/dt = k1/V1(0) * S1 - k2/V2(0) * S2

it's as simple as that!
(since dS2/dt = -dS1/dt as far as reaction S1->S2 is concerned.)
that's the advantage of kinetic laws (and consequently odes) for substances!

> Thus I would really have to check the type of reaction and edit a lot of
> kinetic laws before constructing the ODEs.

the above example shows that this is not necessary.
maybe there are other examples...

> Thus I think it just doesn't work with kinetic laws that don't
> include the compartment explicitly!

i don't know exactly what you mean.
of course kinetic laws with explicit compartment are "nicer",
but you can always transform any kinetic law to a "nice" one...

> But if they include the compartment we can write ODEs for substances but
> have concentrations in the right hand side of the ODEs. We then just need
> to a second array for the concentrations, as discussed before.

if we still use concentrations is a separate problem,
independent of which we discussed until now.
it's a matter of taste, if we eliminate the concentrations  from the kinetic 
laws and get a (possibly nonlinear) dependence on the compartments OR
if we keep concentrations and always have linear dependence on the comp. 

> Rainer

cheers,
stefan.


> On Mon, 19 Jun 2006, Stefan Mueller wrote:
> > let me summarize:
> > kinetic laws have units substance/time,
> > and they depend on concentrations or substances.
> >
> > we SHOULD write down odes for substances,
> > and we COULD eliminate concs completely.
> > (of course, we would have to convert initial concs to substances.)
> >
> > in any case, odes for substances depend on the compartment size
> > (linearly or according to the order of the reaction):
> >
> > order 1:
> > S1->junk: k*[S1]
> > dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
> >
> > order 2:
> > S1+S2->bla: k*[S1]*[S2]
> > dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
> >
> > order 3:
> > S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> > dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 / V^2
> >
> > ...
> >
> > if the kinetic laws included the compartment explicitly (which is
> > desirable), we would have the following situation:
> >
> > order 1:
> > S1->junk: r*[S1]*V
> > dS1/dt =  r * [S1] * V= r * S1
> >
> > order 2:
> > S1+S2->bla: r*[S1]*[S2]*V
> > dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
> >
> > order 3:
> > S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> > dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
> >
> > as stated earlier, odes for substances have got another advantage. for
> > the compartment size we could use either a rate rule or an assignment
> > rule: dV/dt = f(t, V, ...)
> > or
> > V= f(t, V, ...)
> >
> > cheers,
> > stefan.
> >
> > Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> >> andrew:
> >>
> >> for the SBML system
> >> 	S1 in C
> >> 	S1->junk : k * [S1]
> >> 	dC/dt = 1
> >>
> >> the right math is:
> >> 	dS1/dt = k * [S1] / C(0) * C
> >> 	dC/dt = 1
> >> with initial cond:
> >> 	S1(0) = [S1](0) * C(0)
> >> 	C(0)
> >>
> >> notation:
> >> 	S1 ... substance
> >> 	[S1] ... conc
> >>
> >> the "proposed solution" (who proposed this?) has the right units,
> >> but is wrong otherwise... :)
> >>
> >> as rainer and me already discussed, everything would be much clearer,
> >> if the kinetic law included the compartment explicitly:
> >> 	S1 in C
> >> 	S1->junk : C * r * [S1]
> >> 	dC/dt = 1
> >>
> >> with:
> >> 	r ... the "real" chemical rate constant.
> >>
> >> rainer:
> >>
> >> chemists always knew that rate constants are not constant with respect
> >> to pH, T, ... as a physicist i have to be fair! :)
> >>
> >> to conclude: the math is not really tricky.
> >> odes for substances are better suited to handle both multiple and
> >> variable compartments. of course, we have to take care of many
> >> implementation details...
> >>
> >> cheers,
> >> stefan.
> >
> > _______________________________________________
> > 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

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

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

Stefan

i guess it might be different if we have to consider multiple 
compartments

S1 is in V1;
S2 and S3 are in V2;

S1->S2: k1*[S1]
S2->S3: k2*[S2]

[S1] = S1 / V1
[S2] = S2 / V2
[S3] = S3 / V3

Thus the two parameters include two different initial volumes, i.e.:
k1 [volume/sec] = r1 * V1(0)
k2 [volume/sec] = r2 * V2(0)


Currently the ODE would be:

d[S1]/dt = k1 * [S1] / V1
d[S2]/dt = ( k1*[S1] - k2*[S2] ) / V2
d[S3]/dt  k2*[S2]*[S3] / V3

which is only valid for constant compartments V1 and V2.

At least for the second ODE, depending in two different volumes, 
the conversion to ODEs depending on substance instead of concentration is 
not as simple as

> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
(all the ODEs in your email are missing a minus, BTW!)

It should probably be something like:

dS1/dt = - k1 * [S1] / V1(0) * V1 = - k1/V1(0) * S1

as your example before,
but for S2 I get:

dS2/dt = k1/V1(0) * S1/V1 * V2  -  k2/V2(0) * S2

Right?

Thus I would really have to check the type of reaction and edit a lot of 
kinetic laws before constructing the ODEs.

Thus I think it just doesn't work with kinetic laws that don't
include the compartment explicitly!

But if they include the compartment we can write ODEs for substances but 
have concentrations in the right hand side of the ODEs. We then just need 
to a second array for the concentrations, as discussed before.

Rainer



On Mon, 19 Jun 2006, Stefan Mueller wrote:

> let me summarize:
> kinetic laws have units substance/time,
> and they depend on concentrations or substances.
>
> we SHOULD write down odes for substances,
> and we COULD eliminate concs completely.
> (of course, we would have to convert initial concs to substances.)
>
> in any case, odes for substances depend on the compartment size
> (linearly or according to the order of the reaction):
>
> order 1:
> S1->junk: k*[S1]
> dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1
>
> order 2:
> S1+S2->bla: k*[S1]*[S2]
> dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V
>
> order 3:
> S1+S2+S3->wow: k*[S1]*[S2]*[S3]
> dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 / V^2
>
> ...
>
> if the kinetic laws included the compartment explicitly (which is desirable),
> we would have the following situation:
>
> order 1:
> S1->junk: r*[S1]*V
> dS1/dt =  r * [S1] * V= r * S1
>
> order 2:
> S1+S2->bla: r*[S1]*[S2]*V
> dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V
>
> order 3:
> S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
> dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2
>
> as stated earlier, odes for substances have got another advantage. for the
> compartment size we could use either a rate rule or an assignment rule:
> dV/dt = f(t, V, ...)
> or
> V= f(t, V, ...)
>
> cheers,
> stefan.
>
>
> Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
>> andrew:
>>
>> for the SBML system
>> 	S1 in C
>> 	S1->junk : k * [S1]
>> 	dC/dt = 1
>>
>> the right math is:
>> 	dS1/dt = k * [S1] / C(0) * C
>> 	dC/dt = 1
>> with initial cond:
>> 	S1(0) = [S1](0) * C(0)
>> 	C(0)
>>
>> notation:
>> 	S1 ... substance
>> 	[S1] ... conc
>>
>> the "proposed solution" (who proposed this?) has the right units,
>> but is wrong otherwise... :)
>>
>> as rainer and me already discussed, everything would be much clearer,
>> if the kinetic law included the compartment explicitly:
>> 	S1 in C
>> 	S1->junk : C * r * [S1]
>> 	dC/dt = 1
>>
>> with:
>> 	r ... the "real" chemical rate constant.
>>
>> rainer:
>>
>> chemists always knew that rate constants are not constant with respect to
>> pH, T, ... as a physicist i have to be fair! :)
>>
>> to conclude: the math is not really tricky.
>> odes for substances are better suited to handle both multiple and variable
>> compartments. of course, we have to take care of many implementation
>> details...
>>
>> cheers,
>> stefan.
>
>
> _______________________________________________
> sbmlsolver-devel mailing list
> sbm...@li...
> https://lists.sourceforge.net/lists/listinfo/sbmlsolver-devel
>

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

[Stefan M. <ste...@oe...>]

let me summarize:
kinetic laws have units substance/time,
and they depend on concentrations or substances.

we SHOULD write down odes for substances,
and we COULD eliminate concs completely.
(of course, we would have to convert initial concs to substances.)

in any case, odes for substances depend on the compartment size
(linearly or according to the order of the reaction):

order 1:
S1->junk: k*[S1]
dS1/dt = k * [S1] / V(0) * V = k/V(0) * S1

order 2:
S1+S2->bla: k*[S1]*[S2]
dS1/dt = k * [S1] * [S2] / V(0) * V = k/V(0) * S1 * S2 / V

order 3:
S1+S2+S3->wow: k*[S1]*[S2]*[S3]
dS1/dt = k * [S1] * [S2] * [S3] / V(0) * V = k/V(0) * S1 * S2 * S3 / V^2

...

if the kinetic laws included the compartment explicitly (which is desirable),
we would have the following situation:

order 1:
S1->junk: r*[S1]*V
dS1/dt =  r * [S1] * V= r * S1

order 2:
S1+S2->bla: r*[S1]*[S2]*V
dS1/dt = r * [S1] * [S2] * V= r * S1 * S2 / V

order 3:
S1+S2+S3->wow: r*[S1]*[S2]*[S3]*V
dS1/dt = r * [S1] * [S2] * [S3] * V = r * S1 * S2 * S3 / V^2

as stated earlier, odes for substances have got another advantage. for the 
compartment size we could use either a rate rule or an assignment rule:
dV/dt = f(t, V, ...)
or
V= f(t, V, ...)

cheers,
stefan.


Am Montag 19 Juni 2006 11:17 schrieb Stefan Mueller:
> andrew:
>
> for the SBML system
> 	S1 in C
> 	S1->junk : k * [S1]
> 	dC/dt = 1
>
> the right math is:
> 	dS1/dt = k * [S1] / C(0) * C
> 	dC/dt = 1
> with initial cond:
> 	S1(0) = [S1](0) * C(0)
> 	C(0)
>
> notation:
> 	S1 ... substance
> 	[S1] ... conc
>
> the "proposed solution" (who proposed this?) has the right units,
> but is wrong otherwise... :)
>
> as rainer and me already discussed, everything would be much clearer,
> if the kinetic law included the compartment explicitly:
> 	S1 in C
> 	S1->junk : C * r * [S1]
> 	dC/dt = 1
>
> with:
> 	r ... the "real" chemical rate constant.
>
> rainer:
>
> chemists always knew that rate constants are not constant with respect to
> pH, T, ... as a physicist i have to be fair! :)
>
> to conclude: the math is not really tricky.
> odes for substances are better suited to handle both multiple and variable
> compartments. of course, we have to take care of many implementation
> details...
>
> cheers,
> stefan.

[SOSlib-devel] variable compartment bug: example

[Stefan M. <ste...@oe...>]

andrew:

for the SBML system
	S1 in C
	S1->junk : k * [S1]
	dC/dt = 1

the right math is:
	dS1/dt = k * [S1] / C(0) * C
	dC/dt = 1
with initial cond:
	S1(0) = [S1](0) * C(0)
	C(0)

notation:
	S1 ... substance
	[S1] ... conc

the "proposed solution" (who proposed this?) has the right units,
but is wrong otherwise... :)

as rainer and me already discussed, everything would be much clearer,
if the kinetic law included the compartment explicitly:
	S1 in C
	S1->junk : C * r * [S1]
	dC/dt = 1

with:
	r ... the "real" chemical rate constant.

rainer:

chemists always knew that rate constants are not constant with respect to pH, 
T, ... as a physicist i have to be fair! :)

to conclude: the math is not really tricky.
odes for substances are better suited to handle both multiple and variable 
compartments. of course, we have to take care of many implementation 
details...

cheers,
stefan.

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
>