different simulation results from SimBiology code and Matlab code about covid19-qsp-model HOT 10 CLOSED

Thank you for your feedback.
My preliminary analysis shows there are also the difference between DBSolve and Sibiology simulations.
It looks like some problems in Simbiology version.
We are working on it.

from covid19-qsp-model.

To investigate the reproducibility I simulate in different formats and software comparing with DBSolve format as a standard:

• format: SLV, simulation in DBSolve - standard
• format: Simbio, simulate in Simbiology - other results
• format: SBML, simulations in Simbiology - other results
• format: SBML, simulation in DBSolve - same results
• format: SBML, simulation in COPASY - same results

Both Simbiology and SBML simulated in Simbiology (version R2018a) show the other results.

My next step is to try pure Matlab for simulations.

from covid19-qsp-model.

thanks for the update!
I guess probably something is wrong with setups in SimBiology.
Looking forward to your new results!

from covid19-qsp-model.

I've simulated the Matlab version of the model. The results are the same as in DBSolve. I still do not understand what's wrong with Simbio version of code.

The comparison of simulations can be found in dist branch
https://github.com/insysbio/covid19-qsp-model/tree/dist

from covid19-qsp-model.

that's also what I don't understand.
I also did the simulation in Julia (with some modifications for sure to make it run without SimSolver) and the result is the same as in Matlab.
That's why I also think that something is wrong with SimBiology!

from covid19-qsp-model.

After looking at your model, it looks like V_tran_pc_ipc and V_tran_ipc_vpc are dimensionally inconsistent.

For instance, in the following equation:
dydt(3) = V_mat_pc -V_tran_pc_ipc -V_apo_pc;
dydt(3), V_mat_pc and V_apo_pc have units of kcell/hour (amount/time) whereas V_tran_pc_ipc has units of kcell/liter/hour (concentration/time).
The same happens for dydt(4) and dydt(5).

SimBiology honors the specified units and in this case, since a conversion from concentration to amount is missing, it automatically performs a volume correction and does not report any inconsistencies in dimensionality.
If you modify your MATLAB code to multiply V_tran_pc_ipc and V_tran_ipc_vpc by Vol_alv to convert them in dimensions of to amount/time, you will get the same results as in SimBiology.

So, it looks like SimBiology is providing the correct answer given the units specified in the model.

Now, I don’t know whether this inconsistency was due to an omission of the volume scaling or whether the reaction rates are correct, and units are incorrect for some parameters. However, the parameters used in these reaction rates are also used in other reaction rates. So, one can’t change their units without making other reaction rates dimensionally inconsistent.

from covid19-qsp-model.

Thank you, @ChezJe .
You are right. There are some inconsistencies between units in the reactions.
Simbiology solves correctly.
I guess we will fix it soon.
The general problem in finding the bug was that Simbiology tries to guess the correct ODEs based on calculated units.
It looks like our fault in reaction units was fixed by Simbiology internally by multiplying by a volume without any warning.
As a result, the ODE in Simbiology was not the same as in other tools.

from covid19-qsp-model.

We agree that adding volume scaling behind the scenes can lead to confusion.

The bottom line is that users are not writing ODEs in SimBiology but they instead specify reaction rates from which ODEs are then derived. Users can specify reaction rates in concentration or amount per unit time and SimBiology allows it. We don’t consider this an error because many users do specify rates in concentration/time as done in your model.

When building the ODEs, SimBiology converts all reaction rates to dimensionally consistent reaction fluxes adding volume scaling where necessary. Release R2020b includes color coding in the equations view of the modelling app to differentiate species from parameters and from volumes. This will help spot volume scaling.

Given your last change to the units of some parameters to fix this issue, we see that your intent was to specify all reaction rates in amount/time, but there is no way for SimBiology to infer that the units did not match this intent. However, SimBiology is able to detect dimensional inconsistencies in the revised model and reports an error.

from covid19-qsp-model.

don’t consider this an error because many users do specify rates in concentration/time as done in your model.

Of course it is not an error. I understand this is an approach or a feature. But I guess it might act as a pitfall in many cases. And it does.
It is not evident how to scale a reaction in a proper way if user forget about it. What if we have transport between compartments, which is the right compartment to multiply?

One of the possible solution in this situation would be the optional Unit property in Reaction object. In that case a user could set the unit directly and find inconsistencies in reactions if they exist.

from covid19-qsp-model.

Thank you @ChezJe, @ZZP12 for your help.
There is no difference between dbsolve, simbio and matlab.
The issue have been solved.
See dev branch and the code in master.

from covid19-qsp-model.