Add multi-substrate Michaelis–Menten kinetics

doc/modelling/reaction/michaelis_menten_kinetics.rst

@@ -16,27 +16,31 @@ where :math:`S` is the stoichiometric matrix and the fluxes are given by
 .. math::
 
     \begin{aligned}
-        \nu_j = \frac{\mu_{\mathrm{max},j} \, c_S}{k_{\mathrm{MM},j} + c_S},
+        \nu_j = \mu_{\mathrm{max},j} \prod_{i \in \mathcal{S}_j} \frac{{c_i}}{k_{\mathrm{MM},j} + c_i},
     \end{aligned}
 
 where
 
 - :math:`j` is the reaction index,
-- :math:`c_S` is the substrate component,
+- :math:`\mathcal{S}_j` is the set of substrate indices for reaction :math:`j`,
 - :math:`\mu_{\mathrm{max},j}`, is the limiting rate approached by the system at saturation,
-- :math:`k_{\mathrm{MM},j}` is the Michaelis constant, which is defined as the concentration of substrate at which the reaction rate is half ov :math:`\mu_{\mathrm{max},j}`.
+- :math:`k_{\mathrm{MM},j}` is the Michaelis constant, which is defined as the concentration of substrate at which the reaction rate is half of :math:`\mu_{\mathrm{max},j}`.
 
 
-In addition, the reaction might be inhibited by other components.
+In addition, the reaction might be inhibited by other components due to noncompetitive inhibition.
 In this case, the flux has the form
 
 .. math::
 
     \begin{aligned}
-        \nu_j = \frac{\mu_{\mathrm{max},j} c_S}{k_{\mathrm{MM},j} + c_S} \cdot \frac{1}{1 + \sum_i c_{Ii}/k_{I,j,i}}.
+        \nu_j = \mu_{\mathrm{max},j} \prod_{i \in \mathcal{S}_j}\frac{c_i}{k_{\mathrm{MM},j} + c_i} \cdot \frac{1}{1 + \sum_{k = 1}^{n} \frac{c_{k}}{/k_{I,j,k}}}.
     \end{aligned}
 
-Here, :math:k_{\mathrm{I},j,i} is the inhibition constant with respect to component :math:i in reaction :math:j.
-If :math:k_{\mathrm{I},j,i} \leq 0, component :math:i does not act as an inhibitor.
+Here, :math: `k_{\mathrm{I},j,k}` is the inhibition constant with respect to component :math:`k` in reaction :math:`j`.
+If :math: `k_{\mathrm{I},j,k} \leq 0`, component :math:`k` does not act as an inhibitor.
 
 Note: Currently, the model does not allow substrates to function as inhibitors.
+
+References
+^^^^^^^^^^
+.. [1] Irwin H. Segel (1975). Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady State Enzyme Systems.

src/libcadet/model/reaction/MichaelisMentenReaction.cpp

@@ -142,7 +142,10 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 			const unsigned int nReactions = numElements / nComp;
 
 			_stoichiometryBulk.resize(nComp, nReactions);
-			_idxSubstrate = std::vector<int>(nReactions, -1);
+			_idxInhibitor.resize(nReactions);
+			_idxSubstrate.resize(nReactions);
+			_oldInterface = false;
+
 		}
 
 		return true;
@@ -157,7 +160,9 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 	ParamHandler_t _paramHandler; //!< Handles parameters and their dependence on external functions
 
 	linalg::ActiveDenseMatrix _stoichiometryBulk;
-	std::vector<int> _idxSubstrate;
+	std::vector<std::vector<int>> _idxSubstrate;
+	std::vector<std::vector<int>> _idxInhibitor;
+	bool _oldInterface;
 
 	virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx)
 	{
@@ -167,6 +172,11 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 		if ((_stoichiometryBulk.columns() > 0) && ((_paramHandler.vMax().size() < _stoichiometryBulk.columns()) || (_paramHandler.kMM().size() < _stoichiometryBulk.columns())))
 			throw InvalidParameterException("MM_VMAX and MM_KMM have to have the same size (number of reactions)");
 
+		if (_stoichiometryBulk.rows() > 1 && (_paramHandler.kMM().size() == _stoichiometryBulk.columns() * _stoichiometryBulk.rows()))
+			_oldInterface = false;
+		else
+			_oldInterface = true;
+
 		if ((_stoichiometryBulk.columns() > 0) && (_paramHandler.kInhibit().size() < _stoichiometryBulk.columns() * _nComp))
 			throw InvalidParameterException("MM_KI have to have the size (number of reactions) x (number of components)");
 
@@ -180,25 +190,43 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 			std::copy(s.begin(), s.end(), _stoichiometryBulk.data());
 
 			// Find substrate index (first educt)
-			for (int i = 0; i < _stoichiometryBulk.columns(); ++i)
+			_idxSubstrate.clear();
+			_idxInhibitor.clear();
+			const std::vector<double> ki = paramProvider.getDoubleArray("MM_KI");
+			for (int r = 0; r < _stoichiometryBulk.columns(); ++r)
 			{
-				_idxSubstrate[i] = -1;
+				std::vector<int> idxSubstrateReaction_r;
+				std::vector<int> idxInhibitorReaction_r;
+				for (int c = 0; c < _nComp; ++c)
+				{
+					if (_stoichiometryBulk.native(c, r) < 0.0)
+						idxSubstrateReaction_r.push_back(c);
+
+					double kI = ki[_nComp * r + c];
+					if (kI > 0.0)
+						idxInhibitorReaction_r.push_back(c);
 
-				for (int j = 0; j < _nComp; ++j)
+				}
+				if(idxSubstrateReaction_r.empty())
+					throw InvalidParameterException("No substrate found in reaction " + std::to_string(r));
+
+				//throw error if inhibitor is substrate
+				if(!idxInhibitorReaction_r.empty())
 				{
-					if (_stoichiometryBulk.native(j, i) < 0.0)
+					// check if inhibitor is substrate
+					for (int idxSubstrateReaction_r : idxSubstrateReaction_r)
 					{
-						_idxSubstrate[i] = j;
-						break;
+						if(std::find(idxInhibitorReaction_r.begin(), idxInhibitorReaction_r.end(), idxSubstrateReaction_r) != idxInhibitorReaction_r.end())
+							throw InvalidParameterException("Michaelis-Menten reaction: Inhibitor is also substrate in reaction " + std::to_string(r) + "this is not supported yet");
 					}
 				}
-				if (_idxSubstrate[i] == -1)
-					throw InvalidParameterException("Michaelis Menten: No substrate found in reaction " + std::to_string(i));
+
+				_idxSubstrate.push_back(idxSubstrateReaction_r);
+				_idxInhibitor.push_back(idxInhibitorReaction_r);
 			}
 
 		}
 		registerCompRowMatrix(_parameters, unitOpIdx, parTypeIdx, "MM_STOICHIOMETRY_BULK", _stoichiometryBulk);
-
 		return true;
 	}
 
@@ -213,30 +241,31 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 		BufferedArray<flux_t> fluxes = workSpace.array<flux_t>(_stoichiometryBulk.columns());
 		for (int r = 0; r < _stoichiometryBulk.columns(); ++r)
 		{
-			const int idxSubs = _idxSubstrate[r];
-			if (idxSubs == -1)
+			int nSub = _idxSubstrate[r].size();
+			int nInh = _idxInhibitor[r].size();
+			flux_t vProd = 1.0;
+			for (int sIdx = 0; sIdx < nSub; sIdx++)
 			{
-				fluxes[r] = 0.0;
-				continue;
-			}
-
-			flux_t inhSum = 0.0;
-			for (int comp = 0; comp < _nComp; ++comp)
-			{
-				const flux_t kI = static_cast<typename DoubleActiveDemoter<flux_t, ParamType>::type>(p->kInhibit[_nComp * r + comp]);
-				if (kI <= 0.0)
-					continue;
+				int s = _idxSubstrate[r][sIdx];
+				flux_t inhSum = 0;
+				for (int iIdx = 0; iIdx < nInh; iIdx++)
+				{
+					int i = _idxInhibitor[r][iIdx];
+					flux_t ki_rsi = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kInhibit[_nComp * r + i]);
+					inhSum += y[i] / ki_rsi;
+				}
 
-				if(comp == idxSubs)
-					throw InvalidParameterException("Michaelis Menten: Inhibition of substrate is not supported yet");
-				inhSum += y[comp]/kI;
-			}
+				flux_t kMM_rs = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kMM[r]);
+				if (!_oldInterface)
+					// kMM_rs = kMM[r][s]
+					 kMM_rs = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kMM[_nComp * r + s]);
 
-			const flux_t vMax = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->vMax[r]);
-			const flux_t kMM = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->kMM[r]);
+				vProd *= y[s] / ((kMM_rs + y[s]) * (1+inhSum));
 
-			fluxes[r] = vMax * y[idxSubs] / ((kMM + y[idxSubs]) * (1 + inhSum));
+			}
 
+			flux_t vmax_r = static_cast<typename DoubleActiveDemoter<flux_t, active>::type>(p->vMax[r]);
+			fluxes[r] = vmax_r * vProd;
 		}
 
 		// Add reaction terms to residual
@@ -266,52 +295,133 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase
 
 		for (int r = 0; r < _stoichiometryBulk.columns(); ++r)
 		{
-			const int idxSubs = _idxSubstrate[r];
-			if (idxSubs == -1)
-				continue;
+			// Calculate flux and inhibitor sum
+			int nSub = _idxSubstrate[r].size();
+			int nInh = _idxInhibitor[r].size();
 
-			double inhSum = 0.0;
-			for (int comp = 0; comp < _nComp; ++comp)
+			double flux = 1.0;
+			std::vector<double> substratFlux(nSub, 0.0);
+			std::vector<double> inhSumOfSub(nSub, 0.0);
+
+			for (int sIdx = 0; sIdx < nSub; sIdx++)
 			{
-				const double kI = static_cast<double>(p->kInhibit[_nComp * r + comp]);
-				if (kI <= 0.0)
-					continue;
+				int s = _idxSubstrate[r][sIdx];
+				double inhSum = 0.0;
+				for (int iIdx = 0; iIdx < nInh; iIdx++)
+				{
+					int i = _idxInhibitor[r][iIdx];
+					double kI_rsi = static_cast<double>(p->kInhibit[_nComp * r + i]);
 
-				if (comp == idxSubs)
-					throw InvalidParameterException("Michaelis Menten: Substrate-inhibition is not supported yet");
-				inhSum += y[comp] / kI;
-			}
+					inhSum += y[i] / kI_rsi;
+				}
 
+				double kMM_rs = static_cast<double>(p->kMM[r]);
+				if (!_oldInterface)
+					// kMM_rs = kMM[r][s]
+					 kMM_rs = static_cast<double>(p->kMM[_nComp * r + s]);
 
-			const double vMax = static_cast<double>(p->vMax[r]);
-			const double kMM = static_cast<double>(p->kMM[r]);
+				double sFlux = y[s] / ((kMM_rs + y[s]) * (1 + inhSum));
+				flux *= sFlux;
+				// save flux and inhibitor sum for each substrate
+				substratFlux[sIdx] = sFlux;
+				inhSumOfSub[sIdx] = inhSum;
+			}
 
-			const double nom = vMax * y[idxSubs];
-			const double denom = (kMM + y[idxSubs]);
+			double vmax_r = static_cast<double>(p->vMax[r]);
+			flux *= vmax_r;
 
-			// Add gradients to Jacobian
-			// for each substrate component
-			const double dvds = vMax / denom * (1.0 - y[idxSubs] / denom) * 1 / (1 + inhSum);
-			RowIterator curJac = jac;
-			for (int row = 0; row < _nComp; ++row, ++curJac)
+			// calculate derivative for every component jac_flux = dv/dy
+			for (int comp = 0; comp < _nComp; ++comp)
 			{
-				const double colFactor = static_cast<double>(_stoichiometryBulk.native(row, r)) * factor;
-				curJac[idxSubs - static_cast<int>(row)] += colFactor * dvds;
-			}
+				// 1. Check if cIdx is substrate or inhibitor
+				bool isSubstrate = false;
+				int substrateIdx = -1;
+				for (int sIdx = 0; sIdx < nSub; sIdx++)
+				{
+					if (comp == _idxSubstrate[r][sIdx])
+					{
+						isSubstrate = true;
+						substrateIdx = sIdx;
+						break;
+					}
+				}
 
-			// for each inhibitor component
-			curJac = jac;
-			for (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)
+				bool isInhibitor = false;
+				int inhibitorIdx = -1;
+				for (int iIdx = 0; iIdx < nInh; ++iIdx)
+				{
+					if (comp == _idxInhibitor[r][iIdx])
+					{
+						isInhibitor = true;
+						inhibitorIdx = iIdx;
+						break;
+					}
+				}
+
+				double dvdy = 0.0;
+
+				// case 1: comp is a substrate and not an inhibitor
+				// dvds = vmax * prod_i!=s (y[i]/(kMM_rs + y[i])*(1+inhSum)) * df/ds
+				// df/ds = ((kMM_rs + y[i])*(1+inhSum)) - y[i](1+inhSum) / ((kMM_rs + y[i])*(1+inhSum))^2
+				//       = (kMM_rs) / ((kMM_rs + y[i])*(1+inhSum)^2)
+				if (isSubstrate && !isInhibitor)
+				{
+					int s = comp; // comp is a substrate
+					double kMM_rs = static_cast<double>(p->kMM[r]);
+					if (!_oldInterface)
+						kMM_rs = static_cast<double>(p->kMM[_nComp * r + s]);
+
+					double inhSum = inhSumOfSub[substrateIdx];
+					double dnom = (kMM_rs + y[s])* (kMM_rs + y[s])*(1+inhSum);
+
+					double dsubFluxds = kMM_rs/dnom ;
+					double factor = flux / substratFlux[substrateIdx];
+
+					dvdy += factor * dsubFluxds;
+
+				}
+
+				//case 2: comp is a inhibitor and not a substrate
+				// dvdI = vmax * sum_L(prod_!i=l (y[i]/(kMM_rs + y[i])*(1+inhSum)) * df/dI
+				// L is the set of all substrates which are inhibited by I
+				// df/dI =  - ((kMM_rs+y[s])/kI_ri)/((kMM_rs + y[s])*(1+inhSum))^2
+				//		 =  - (y[s])/((kMM_rs + y[s])*(1+inhSum)^2 kI_ri)
+				if (isInhibitor && !isSubstrate)
 				{
-					const double kI = static_cast<double>(p->kInhibit[_nComp * r + comp]);
-					if (kI <= 0.0)
-						continue;
+					int i = comp;
+					double kI_ri = static_cast<double>(p->kInhibit[_nComp * r + i]);
+
+					for (int sIdx = 0; sIdx < nSub; sIdx++)
+					{
+						int s = _idxSubstrate[r][sIdx];
+						double kMM_rs = static_cast<double>(p->kMM[r]);
+						if (!_oldInterface)
+							kMM_rs = static_cast<double>(p->kMM[_nComp * r + s]);
 
-					double dvdi = - (vMax * y[idxSubs] * (kMM + y[idxSubs])) / (denom * denom * kI);
-					curJac[comp - static_cast<int>(row)] += colFactor * dvdi;
+						double inhSum = inhSumOfSub[sIdx];
+
+						double denom = (kMM_rs + y[s])*(1+inhSum)* (1 + inhSum);
+						double dinhFluxdi = - (y[s])/ (denom * kI_ri);
+						double factor = flux / substratFlux[sIdx];
+
+						dvdy += factor * dinhFluxdi;
+					}
+				}
+
+				if (isInhibitor && isSubstrate)
+				{
+					throw(InvalidParameterException("Michaelis-Menten reaction: Inhibitor and substrate in the same reaction is not supported yet"));
+				}
+
+				if (std::abs(dvdy) < 1e-18)
+					continue;
+
+				// Add gradients to Jacobian
+				RowIterator curJac = jac;
+				for (int row = 0; row < _nComp; ++row, ++curJac)
+				{
+					const double colFactor = static_cast<double>(_stoichiometryBulk.native(row, r)) * factor;
+					curJac[comp - static_cast<int>(row)] += colFactor * dvdy;
 				}
 			}
 		}