[SOSlib-devel] Request for help with a multicompartment problem

[Michael H. <mh...@ca...>]

Rainer and the rest,

Thanks for finding those errors in the examples in the SBML
spec; we'll try to get them out.

But since you have been thinking about this topic, I'd like
to ask you how you think the following class of problems can
be solved.  It has been confounding Nicolas and myself.

Suppose we have a reaction system consisting of 3
compartments, with species A in Va, B in Vb, and C in Vc
(I'm using the volumes as the names of the compartments),
and the reaction is

   A + B -> C

The compartments have constant volume; no varying volumes
here.  Suppose further that the reaction is said to obey the
simplest possible kind of kinetics, a kind of mass action
(but a weird one because multiple compartments are
involved).  So let's say the rate is

  k * [A] * [B]

Here, k must have units of a 2nd order rate constant, which
is 1/(Molar * time) = (volume)/(substance*time).  As initial
values, we might be given the values of concentrations [A],
[B], [C], or we might be given the substance amounts; it
won't matter as we'll see below.  The key point is that the
reactants are in different compartments.

Now we want to construct the ODEs.

First let's start naively just to make sure all the grounds
are covered.  The first attempt we might make is to say that
the rates of change are like in traditional kinetics:

  d[A]/dt = -k * [A] * [B]
  d[B]/dt = -k * [A] * [B]
  d[C]/dt =  k * [A] * [B]

But this can easily be shown to be wrong.  Borrowing an
example that Nicolas has written for the l2v2 spec, let's
say that volume Vc = 2*Va (i.e., Vc is twice the volume of
Va).  Then if we calculate the change of substance amounts
(moles or items) in compartments Va and Vc (and here I'm
using "A" to stand for the amount of substance of species A,
"C" to stand for the amount of substance of species C):

  dA/dt = Va*d[A]/dt = Va * (-k * [A] * [B]^2)

  dC/dt = Vc*d[C]/dt = Vc * (k * [A] * [B]^2)
                     = 2 * Va * (k * [A] * [B]^2)

So we just created matter out of nothing.  This simply
proves what we already knew, which is that we have to
express the rate in terms of substance changes, not
concentration changes.

OK, fine, so let's try to do that.  Since we know the
compartment volumes and we can know the concentrations of
each species, we have the relationships

  [A] = A/Va
  [B] = B/Vb
  [C] = C/Vc

We then substitute this into the ODEs:

  d[A]/dt = -k * (A/Va) * (B/Vb)
  d[B]/dt = -k * (A/Va) * (B/Vb)
  d[C]/dt =  k * (A/Va) * (B/Vb)

So far this is perfectly valid, but it didn't help anything;
the overall expressions are still in terms of
concentrations.  We need to make them to be about species
amount, not concentrations.

But how?

Attempt #1: Multiply each ODE by the volume in which that
species is located.  Let's take for example d[C]/dt, and
multiply it by the volume in which [C] is located:

  dC/dt = Vc*d[C]/dt = k * (A/Va) * (B/Vb) * Vc

This seems reasonable; after all, C is located in Vc, and
Vc*[C] must give me number of molecules of C.  But it can't
be right.  Suppose I change the model by doubling the size
of Vc, without changing anything else in the system.
According to the above, the rate of appearance of molecules
of C will double if I do that, which is impossible because
actually the number of molecules must be solely determined
by the reactants.

Attempt #2: The conclusion from the above must be that
multiplying by Vc is the wrong way around.  The core rate
expression must be reformulated in terms of substance per
time, not concentration per time, and *then* concentrations
can be made out of it by dividing by the volumes as needed.
So let's try to do that.

  new rate expr = k * ([A]*Va) * ([B]*Vb)

This gets rid of the concentration terms, but now the units
are wrong, because they are (volume * substance)/time rather
than substance/time.  Hmmm.  Could this be a clue?  Maybe we
just need to divide by a volume?  That seems promising.  So:

  dA/dt = ( -k * ([A]*Va) * ([B]*Vb) )/Va
  dB/dt = ( -k * ([A]*Va) * ([B]*Vb) )/Vb
  dC/dt = (  k * ([A]*Va) * ([B]*Vb) )/Vc

The units look right.  But wait, this has the same problem
as before: if Vc changes and nothing else in the system
changes, the rate of appearance of C will change.

Attempt #3: Maybe we simply divided by the wrong volume in
the equations above; maybe it should be always Va or Vb.
Oh, but, there are two choices: Va and Vb.  If we can only
use one, how do we pick which one we use?

Attempt #4: Ok, so the volume division idea isn't working.
What if the issue is not about dividing by volumes, but
instead, the problem is that the units of k need to be
changed?  Maybe what we really need is

  dA/dt = -k' * ([A]*Va) * ([B]*Vb)
  dB/dt = -k' * ([A]*Va) * ([B]*Vb)
  dC/dt =  k' * ([A]*Va) * ([B]*Vb)

Then k' can be used to sort out the units.  Except.... what
*should* be the correction applied to k?  How can we create
k' systematically?  What if the rate law is something
entirely different -- without a k', or with other constants?
What will the correction be then?


So at this point, we're stuck.  The conclusion must be one
of the following:

  Option (a): Conclude that attempt #4 is the right way, and
    we need to work out how to construct a proper correction
    term for the general case.

  Option (b): Conclude that there is a fundamental error in
    one or more of the attempts #1-4, or there is another
    way out.

  Option (c): Conclude that we must declare the situation
    impossible and require that all reactants must be in the
    same compartment.

Unfortunately, Nicolas points out that (c) is violated by
some existing models -- they already use multiple
compartments for reactants -- so people must believe it's
doable.  Option (c) is not a way out.

What is your opinion?  What is the solution?  How should
this multicompartment model be solved?  

MH