CADET-Core
CADET-Core copied to clipboard
cadet
Simple Michaelis-Menten reaction kinetics
Implements liquid phase Michaelis-Menten reaction kinetics of the form
S \nu,
where $S$ is the stoichiometric matrix and the fluxes are given by
\nu_i = \frac{\mu_{\mathrm{max},i} c_S}{k_{\mathrm{MM},i} + c_S}.
The substrate component $c_S$ is identified by the index of the first negative entry in the stoichiometry of this reaction.
In addition, the reaction might be inhibited by other components. In this case, the flux has the form
\nu_i = \frac{\mu_{\mathrm{max},i} c_S}{k_{\mathrm{MM},i} + c_S} \prod_j \frac{k_{\mathrm{I},i,j}}{k_{\mathrm{I},i,j} + c_{\mathrm{I},j}}.
The value of $k_{\mathrm{I},i,j}$ decides whether component $j$ inhibits reaction $i$. If $k_{\mathrm{I},i,j} < 0$, the component does not inhibit the reaction.
Tests pass, Jacobian appears to be correct. Please perform some tests to see if the model works or what should be changed.
What if $k_{I,i,j}=0$? Wouldn't it make sense to say that compoent $j$ then also doesn't inhibit reaction $i$?
What if kI,i,j=0? Wouldn't it make sense to say that compoent j then also doesn't inhibit reaction i?
Changed it to use <= 0 instad of < 0.
I will provide some examples and Python code from my lecture on computational biotechnology for the testing.
The tests can also be based on the existing implementation which is based on micro-kinetics. e.g., with $k_r=0$ and $E_0=1$ such that $\mu_{max}=k_{cat} and $K_{MM}=\frac{1}{k_f}$.
Moreover, we could easily cover Hill kinetics by including an exponent.
mm_reaction.zip So I am currently testing this with an example from Erics lecture, and I want to make sure I set it up correctly.
The example given is:
with
I calculated B to be 100.65 in steady state where B_dot = u-v = 0.
I set up my CSTR (V = 1000) with a reaction model following Sams entry post as
vmax_A = 400
vmax_B = 100
KA = 2.5
KB = 5.0
IB = 35
A = 30.0
B = 100.65
init_c = [A, B]
#Reaction i = rows ?
#Component j = columns ?
model.root.input.model.unit_000.reaction_model = 'MICHAELIS_MENTEN'
model.root.input.model.unit_000.reaction_bulk.mm_kmm = [KA, KB]
model.root.input.model.unit_000.reaction_bulk.mm_ki = [-1.0, IB, #Component B Inhibits Reaction 1
-1.0, -1.0]
model.root.input.model.unit_000.reaction_bulk.mm_vmax = [vmax_A, vmax_B]
model.root.input.model.unit_000.reaction_bulk.mm_stoichiometry_bulk = [-1.0, 0.0, #Substrate in Reaction 1 is A, so -1.0?
0.0, 1.0] #Product of Reaction 2 is B, so 1.0?
from my understanding I would now expect to see Components A (blue) and B (orange) in steady state, but I get
which looks like I set up the stoichiometry wrong. I played around a little and found in the documentation of the other reaction models, that there the components are noted in the rows and the reactions in the columns, which confused me as here reactions are noted as i and components as j and I could not find a documentation for the mm reaction yet. Playing around I managed to get different profiles, but not what I expected.
I don't understand the sentence "The substrate component is identified by the index of the first negative entry in the stoichiometry of this reaction." because for the S-Matrix it seems like the opposite, negative entry gets produced and positive entry is consumed.
So first I want to make sure I understand the parameters right and set up the matrices correctly.
@sleweke-bayer
We now broke it down to a MRE, that shows that although A is marked as a substrate, A is synthesized while B is consumed. Maybe a sign error in the implementation of the reaction in CADET? The notebook is attached. MichalelisMenten_MRE.zip
# Variables
vmax_A = 400
KA = 2.5
A = 100.0
B = 0.0
init_c = [A, B]
# Reaction
model.root.input.model.unit_000.reaction_model = 'MICHAELIS_MENTEN'
model.root.input.model.unit_000.reaction_bulk.mm_kmm = [KA]
model.root.input.model.unit_000.reaction_bulk.mm_ki = [-1.0]
model.root.input.model.unit_000.reaction_bulk.mm_vmax = [vmax_A]
model.root.input.model.unit_000.reaction_bulk.mm_stoichiometry_bulk = [-1.0,
1.0]
Sam and I could not reproduce the error based solely on the h5 file. Will check my files again and come back to this.
I cannot reproduce the problem. The simple Michaelis-Menten test case appears to work for me.
Found the issue, I was on another (outdated) feature/michaelis_menten branch... Thanks for taking the time today Sam, at least we know now that it seems to be working! Will continue testing and report.
I created the basics for a test that checks the mm kinetics against mass action law microcinetics. It seems to work just fine!
´mm_gt.h5´ uses the new michaelis menten reaction model, ´mm_microkinetics_test.h5´ uses microkinetics and the old mass action law reaction model to get the same results within a very tight error tolerance (<1e-5) (component 1 and 2 of the michaelis menten solution equal component 2 and 4 of the microcinetics solution ).
Now this has to be integrated into a test where the solutions are compared. As I am not super familiar with C++ I would either need help (and some time) to do this or give this last task to someone else.
We added a Numerical Reference Test comparing a case from Erics Lecture calculated with FD against a simulation.
From our side the branch is now ready to merge
I just realized that this PR didn't contain any documentation or interface specification. In future, we should make sure to add that before merging. Or at least create a separate issue/PR.
This could also be a good idea for a template @ronald-jaepel
