Kinetic Models

Introduction to Kinetic Models

A kinetic model can be description of metabolism, which is constrained by mass conservation where the flows are limited by rate laws. As an example, let’s observe a hypothetical metabolite \(A\) with the metabolite flows defined in the figure below. The rate of change of \(A\) is given by the following equation:

\[\frac{dN_A}{dt} = V \times v_1(A,\theta_1) - V \times v_2(A, \theta_2)\]

If we know the functions for \(v_1\) and \(v_2\) we can determine the concentration of \(A\) over time, given an initial state \(A_0\), the compartment volume \(V\), and kinetic parameters \(\theta_i\).

\[A_t = \int_0^t v_1(A,\theta_1) - v_2(A, \theta_2) dt\]

A metabolic network is formed by connecting metabolites through reactions, thus creating a tightly coupled network. For example, if \(v_2\) leads to a metabolite \(B\), the concentration of \(B\) is dependent on \(A\). The resulting system of ODEs must be solved simultaneously.

The network can be represented with a stoichiometric matrix (S) in which each row represents a metabolite, each column a reaction and each entry represents the corresponding stoichiometric coefficient for the metabolite in the reaction. The following figure is a three reaction network which is represented in a stoichiometric matrix


Using the stoichiometric matrix, we can express the system of ODEs in concise matrix form:

\[\begin{split}\frac{dC_i}{dt} &= \Sigma_j S_{ij} v_j \\ \frac{dC}{dt} &= Sv\end{split}\]

where \(S_{ij}\) is the stoichiometric coefficient of metabolite \(i\) in reaction \(j\) and \(v_j\) is the flux for reaction \(j\).

Steady State Assumption

Maud evaluates the metabolism to a psuedo-steady-state such that:

\[\frac{dC_i}{dt} = 0\]

or in other terms,

\[Sv = 0\]

The following assumptions are made about the system:

  1. Enzyme concentrations are constant over time, and,

  2. The flows In and Out of a metabolite due to reactions are much higher than the dilution rate due to growth.

Assumption 2. can be relaxed by accounting for the growth rate, however, we leave this to the user to implement.

Flux vector representation

In a kinetic model, each flux is expressed as a rate law:

\[v_j = f(E,C,k, \Delta G^0)\]


  • \(j\) is the reaction \(j\);

  • \(E\) is the enzyme concentration;

  • \(C\) are the metabolite concentrations;

  • \(k\) are the kinetic parameters;

  • \(\Delta G^0\) is the standard Gibbs energy of reaction

The rate laws can incorporate allosteric effectors, competitive inhibition, reaction mechanisms, and varied protein concentrations. Common rate laws include: mass action kinetics where the rate is a scaled-product of the reactant concentrations, or Michaelis-Menten kinetics which reproduce enzyme saturation behavior.

The current supported mechanisms are reversible_modular_rate_law and irreversible_modular_rate_law, which is augmented using the Generalised Monod-Wyman-Changeux formalism for allostery (see the documentation page about enzyme kinetics for details). Additional rate laws can be easily implemented by the user and assigning them flags in the corresponding “” and “model.stan” files.

Out of Scope

Maud’s current purpose is to fit steady state data and currently isn’t structured to accept anything other than a single point. As such, the following applications are not implemented.

Dynamic model fitting

Whilst theoretically possible and useful for circumstances such as substrate pulses, fitting concentrations at multiple timepoints would require a significant restructuring of the current implementation. Therefore, we limit measurement evalutions to steady state evaluations.

Limit Cycles

Parameterisations which result in limit cycles are currently unsupported. A limit cycle is a stable oscillatory solution that doesn’t achieve a steady state, it has been observed in Yeast glycolysis and can be approximated by taking the average of the ossciations as the measurement of the system, as seen in How Yeast Cells Synchronize their Glycolytic Oscillations: A Perturbation Analytic Treatment. However, the current method does not permit for averaging and will indicate that the solution is unstable, regardless of simulation time.

Solving Systems of ODEs

To determine \(Sv = 0\) we evaluate the system of ODEs as an initial value problem by integrating the system over an arbitrary period of time that is defined by the user. Tolerances for steady state evaluation can be altered in the source code under - src/maud/model.stan. This acts as an indicator which will prompt the user if the timepoint selected is insufficient, or may be a result of Limit Cycles. The principal assumption of this method is that there is a unique solution for every initial point, assuming every \(x_0 \gt 0\). By initialising the concentrations as close to the measured values as possible: simulation time is minimised; and of the possible non-unique solutions, the solution that is closest to the initial value and hence measurement will ideally be selected. There is no other implementation in Maud to account for multiple solutions.