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

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

Rainer
=20
I think I understand your English argument

my last email was really just a comment on
how you expressed your argument :-)

I don't exactly understand why the math doesn't work.
CVODE solves any correctly defined ODE system
if CVODE produces the wrong results then we've got the formulae wrong.  =
(I'm just restating the issue slightly differently)

given=20

d[s]/dt =3D f([S])/C         [1]

when=20

dC/dt =3D g([S])             [2]

Either [1] is the correct formulae or it isn't
when
C is varying volume and
f is a function returning substance/time units and
g is a function returning volume/time units and
[S] is the set of all concentrations. =20

Intuitively I can't see why this formulae is wrong
but hey what do I know!

My calculus really isn't up to figuring this out!
sorry I'm a software engineer.  I guess I should know for SBML and all =
that but I don't.

Have we got it wrong?
Perhaps we have.

I don't have problem in changing SOS to using amount units for the =
species variables.  However you will still have a problem when a model =
contains a rate rule in concentration units.  What happens then?
(in that case as na=EFve softy I would multiply by the volume!)

yours Andrew



> -----Original Message-----
> From: Rainer Machne [mailto:ra...@tb...]=20
> Sent: 13 June 2006 15:09
> To: Andrew Finney
> Cc: SOSlib Development
> Subject: RE: [SOSlib-devel] variable compartment bug[Scanned]
>=20
> Andrew
>=20
>=20
> I have added a corrected and more clear version of the example file.
>=20
>=20
> > All you need to do is ensure that the initialAmount values=20
> are used to=20
> > calculate the initialConcentration values again by dividing by the=20
> > volume.  Have we implemented that?
>=20
> Yes, the internal SBML ode model only contains=20
> initialConcentrations; initialAmounts from the input model=20
> have all been converted by dividing by compartment size.
>=20
>=20
> > if S1 is in a compartment with volume C then
> >
> > d[S1]/dt =3D (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.
>=20
> Yes, that is what we have.
>=20
>=20
> > Why do you need two sets of values?
>=20
> Ok, I'll try to explain in more detail.
>=20
> I change the notation to
> S ... amount of molecule S
> s ... concentration, [S] =3D S/v
> v ... volume
>=20
>=20
> Imagine we start with initialConcentration s(0) =3D 1,
> and integrate from 0 to time 1.
>=20
> What cvode does is to add the concentration increase ds to=20
> the initial=20
> concentration:
>=20
> s(t) =3D s(0) + ds
>=20
> if the volume is doubled during that time e.g. by a rate=20
> rule, ds should
> correctly include the changing volume because its part of the ODE.
>=20
> However, s(0)=3D1 has never been corrected for the changing volume!
>=20
> If the ODE ds/dt was 0 during this time, e.g. because of a missing=20
> activator of the reaction, s(t) would still change to 0.5,=20
> because the=20
> volume has doubled.
>=20
> This is not reflected in SOSlib!
>=20
> To the example XML:
>=20
> > Note the initial volume of the compartment is not specified=20
> in the xml.
>=20
> Sorry, I tried the model with old version which still had the=20
> `missing=20
> parameter bug' from last week.
> See the attached corrected file.
>=20
> The compartment volume increases to 11 because of its rate rule
> d(volume)/dt =3D + 0.1
>=20
> Substance s1 starts with concentration 1 and is totally=20
> consumed by the=20
> irreversible reaction
>=20
> S1 -> S2; kineticLaw =3D k+S1
>=20
> which moves all molecules to S2.
>=20
> Substance S2 would have a final concentration of 1, if the=20
> compartment=20
> still had volume 1, but not for a final volume of 11.
>=20
> The final concentration should be 1/11.
>=20
> The totalConcentration in the example file should also=20
> decrease to 1/11.
>=20
>=20
> SOSlib and the ODE division by volume only corrects the change of=20
> concentration but not the initial concentration, or generally=20
> spoken, the=20
> concentration of the last time step, which also has decreased.
>=20
> The only solution i can think of, is to define ODEs in=20
> substance/time and=20
> do the division independently of the ODE system.
>=20
> Rainer
>=20
>=20
>=20
>=20
>=20
>=20
>=20
>=20
>=20
>=20
> On Tue, 13 Jun 2006, Andrew Finney wrote:
>=20
> > 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 =3D (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=20
> 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=20
> 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=20
> values again by
> > dividing by the volume.  Have
> > we implemented that?
> >
> > Note the initial volume of the compartment is not specified=20
> 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
> >>
> >
>=20