Skip to content

Commit

Permalink
Add arbitrary inhibition to Michaelis-Menten reaction
Browse files Browse the repository at this point in the history 
Inhibition of a component is determined by the inhibition coefficient,
that is, a negative value marks a component as non-inhibiting.
  • Loading branch information
@sleweke-bayer
sleweke-bayer committed Mar 15, 2024
1 parent 7820f14 commit fe80b47
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 20 deletions.
74 changes: 57 additions & 17 deletions src/libcadet/model/reaction/MichaelisMentenReaction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,19 +34,13 @@
[
{ "type": "ScalarReactionDependentParameter", "varName": "vMax", "confName": "MM_VMAX"},
{ "type": "ScalarReactionDependentParameter", "varName": "kMM", "confName": "MM_KMM"},
{ "type": "ScalarReactionDependentParameter", "varName": "kInhibit", "confName": "MM_KI"}
{ "type": "ComponentDependentReactionDependentParameter", "varName": "kInhibit", "confName": "MM_KI"}
]
}
</codegen>*/

/* Parameter description
------------------------
kFwdBulk = Forward rate for reactions in bulk volume
kBwdBulk = Backward rate for reactions in bulk volume
kFwdLiquid = Forward rate for reactions in particle liquid phase
kBwdLiquid = Backward rate for reactions in particle liquid phase
kFwdSolid = Forward rate for reactions in particle solid phase
kBwdSolid = Backward rate for reactions in particle solid phase
*/


Expand Down Expand Up @@ -95,14 +89,27 @@ namespace


/**
* @brief Defines the multi component Langmuir binding model
* @details Implements the Langmuir adsorption model: \f[ \begin{align}
* \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} q_i
* @brief Defines a Michaelis-Menten reaction kinetic with simple inhibition
* @details Implements the Michaelis-Menten kinetics: \f[ \begin{align}
* S \nu,
* \end{align} \f]
* Multiple bound states are not supported.
* Components without bound state (i.e., non-binding components) are supported.
* where \f$ S \f$ is the stoichiometric matrix and the fluxes are given by
* \f[ \begin{align}
* \nu_i = \frac{\mu_{\mathrm{max},i} c_S}{k_{\mathrm{MM},i} + c_S}.
* \end{align} \f]
* The substrate component \f$ c_S \f$ 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
* \f[ \begin{align}
* \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}}.
* \end{align} \f]
* The value of \f$ k_{\mathrm{I},i,j} \f$ decides whether component \f$ j \f$
* inhibits reaction \f$ i \f$. If \f$ k_{\mathrm{I},i,j} < 0 \f$, the component
* does not inhibit the reaction.
*
* See @cite Langmuir1916.
* Only reactions in liquid phase are supported (no solid phase or cross-phase reactions).
* @tparam ParamHandler_t Type that can add support for external function dependence
*/
template <class ParamHandler_t>
Expand Down Expand Up @@ -209,7 +216,16 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
continue;
}

fluxes[r] = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->vMax[r]) * y[idxSubs] / (static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kMM[r]) + y[idxSubs] * (1.0 + y[idxSubs] / static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kInhibit[r])));
fluxes[r] = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->vMax[r]) * y[idxSubs] / (static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kMM[r]) + y[idxSubs]);

for (int comp = 0; comp < _nComp; ++comp)
{
const flux_t kI = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kInhibit[_nComp * r + comp]);
if (kI < 0.0)
continue;

fluxes[r] *= kI / (kI + y[comp]);
}
}

// Add reaction terms to residual
Expand Down Expand Up @@ -243,11 +259,21 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
if (idxSubs == -1)
continue;

double inhibit = 1.0;
for (int comp = 0; comp < _nComp; ++comp)
{
const double kI = static_cast<double>(p->kInhibit[_nComp * r + comp]);
if (kI < 0.0)
continue;

inhibit *= kI / (kI + y[comp]);
}

const double vMax = static_cast<double>(p->vMax[r]);
const double kIn = static_cast<double>(p->kInhibit[r]);
const double kMM = static_cast<double>(p->kMM[r]);
const double kk = kIn * kMM;
const double fluxGrad = kIn * vMax * (kk - sqr(y[idxSubs])) / sqr(y[idxSubs] * (kIn + y[idxSubs]) + kk);
const double denom = kMM + y[idxSubs];
const double fluxGrad = vMax / denom * (1.0 - y[idxSubs] / denom) * inhibit;
const double flux = vMax * y[idxSubs] / denom * inhibit;

// Add gradients to Jacobian
RowIterator curJac = jac;
Expand All @@ -256,6 +282,20 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
const double colFactor = static_cast<double>(_stoichiometryBulk.native(row, r)) * factor;
curJac[idxSubs - static_cast<int>(row)] += colFactor * fluxGrad;
}

curJac = jac;
for (unsigned int row = 0; row < _nComp; ++row, ++curJac)
{
const double colFactor = static_cast<double>(_stoichiometryBulk.native(row, r)) * factor;
for (int comp = 0; comp < _nComp; ++comp)
{
const double kI = static_cast<double>(p->kInhibit[_nComp * r + comp]);
if (kI < 0.0)
continue;

curJac[comp - static_cast<int>(row)] -= colFactor * flux / (kI + y[comp]);
}
}
}
}

Expand Down
23 changes: 20 additions & 3 deletions test/ReactionModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,19 +77,36 @@ TEST_CASE("MassActionLaw kinetic analytic Jacobian vs AD", "[MassActionLaw],[Rea
);
}

TEST_CASE("MichaelisMenten kinetic analytic Jacobian vs AD", "[MichaelisMenten],[ReactionModel],[Jacobian],[AD]")
TEST_CASE("MichaelisMenten kinetic analytic Jacobian vs AD without inhibition", "[MichaelisMenten],[ReactionModel],[Jacobian],[AD]")
{
const unsigned int nBound[] = {1, 2, 1};
const double point[] = {1.0, 2.0, 1.4, 2.1, 0.2, 1.1, 1.8};
cadet::test::reaction::testDynamicJacobianAD("MICHAELIS_MENTEN", 3, nBound,
R"json({
"MM_KMM": [1.0, 2.0, 0.4],
"MM_KI": [1.0, 0.2, 1.5],
"MM_KI": [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0],
"MM_VMAX": [1.0, 0.2, 1.5],
"MM_STOICHIOMETRY_BULK": [ 1.0, -2.0, 3.0,
-1.0, 0.0, -2.0,
0.0, 1.0, 1.0]
})json",
point, 1e-15, 1e-15
);
}
}

TEST_CASE("MichaelisMenten kinetic analytic Jacobian vs AD with inhibition", "[MichaelisMenten],[ReactionModel],[Jacobian],[AD]")
{
const unsigned int nBound[] = {1, 2, 1};
const double point[] = {1.0, 2.0, 1.4, 2.1, 0.2, 1.1, 1.8};
cadet::test::reaction::testDynamicJacobianAD("MICHAELIS_MENTEN", 3, nBound,
R"json({
"MM_KMM": [1.0, 2.0, 0.4],
"MM_KI": [-1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 3.0, 2.0, -1.0],
"MM_VMAX": [1.0, 0.2, 1.5],
"MM_STOICHIOMETRY_BULK": [ 1.0, -2.0, 3.0,
-1.0, 0.0, -2.0,
0.0, 1.0, 1.0]
})json",
point, 1e-15, 1e-15
);
}

0 comments on commit fe80b47

Please sign in to comment.