Hasan,
I believe there is a mismatch in the units of the ODE. According to the SBML website this problem is described as:
The model contains one compartment called "compartment". There are two species named S1 and S2 and two parameters named k1 and k2. The model contains two reactions defined as:
| Reaction |
Rate |
| S1 → 2 S2 |
k1 · S1 · compartment |
| 2 S2 → S1 |
k2 · S2 · S2 · compartment |
The initial conditions are as follows:
| |
Value |
Units |
| Initial amount of S1 |
1.5 × 10−4 |
mole |
| Initial amount of S2 |
0 × |
mole |
| Value of parameter k1 |
0.35 |
second-1 |
| Value of parameter k2 |
180 |
litre mole-1 second-1 |
| Volume of compartment "compartment" |
0.3 |
litre |
The species values are given as amounts of substance to make it easier to use the model in a discrete stochastic simulator, but (as per usual SBML principles) their symbols represent their values in concentration units where they appear in expressions.
I think you are using a0 to represent s1, and a1 to represent s2 in your ODE RHS equations, but to be consistent lets keep s1 and s2.
Your RHS side equations are therefore:
ds1/dt = -(compartment*k1*s1)+(compartment*k2*s2^2)
ds2/dt = -2*(compartment*k2*s2^2)+2*(compartment*k1*s1)
Your initial values for the s1 and s2 are 0.00015 and 0, just as the problem description states. This implies that the units of S1 and s2 are moles.
The graph of s1,s2 VS time from the SBML website also supports this, as the s1 trace begins at 0.00015 and the s2 trace begins at 0.
Now lets look at the units that derive from the first RHS equation:
ds1/dt = -(compartment*k1*s1)+(compartment*k2*s2^2)
units of compartment = litre
units of k1=1/second
units of s1=mole
ds1/dt is the rate of change of s1 over time, so this should be mole/second.
So the first term (compartment*k1*s1) has units of (litre)(1/second)(mole) or (litre*mole/second). This does not match the units of ds1/dt.
The SBML description has a note at the bottom that seems to describe the issue
The species values are given as amounts of substance to make it easier to use the model in a discrete stochastic simulator, but (as per usual SBML principles) their symbols represent their values in concentration units where they appear in expressions.
The symbol s1 in the reaction rate equation you used to derive the ODE has units of concentration, which I think should be Molarity (mole/litre). However, the variable s1 and s2 in the ODE have units of mole, and so should be be scaled by volume. s1->s1/compartment, s2->s2/compartment
This makes the RHS equations:
ds1/dt = -(compartment*k1*(s1/compartment))+(compartment*k2*(s2/compartment)^2)
ds2/dt = -2*(compartment*k2*(s2/compartment)^2)+2*(compartment*k1*(s1/compartment))
simplifying we get:
ds1/dt = -(k1)*s1+(k2/compartment)*s2^2
ds2/dt = 2*(k1)*s1 - 2*(k2/compartment)*s2^2
Which are the equations I found empirically yesterday.
Could this be the problem?
-Jim
