Often biological systems are modeled using Ordinary Differential Equations (ODEs). Biological ODE models are continuous deterministic models in which variables represent concentrations of species (particles). Models are based on average reaction rates, measured from experiments which relate to the change of concentrations over time. With the reaction equation, initial amount of each species and the rate of every reaction a model can be constructed.
There are several methods to model biological systems, some of them are:
Each possible state of the system is described by a linear autonomous ODE. Solution of the kth equation at time t is a real number giving the probability of the model being in that particular state at time t. Generally, CME method has a high dimension that the model cannot handle analytically or computationally. A description with motivating example to derive the CME is presented in “D. J. Higham, “Modeling and Simulating Chemical Reactions,” ”.
A set of non-linear, autonomous Stochastic Differential Equations (SDEs). One SDE for each chemical species. Solution of the jth equation at time t is a real-valued random variable representing the amount of species j at time t. A description and example to derive the CLE is presented in “D. J. Higham, “Modeling and Simulating Chemical Reactions,” ”.
A set of non-linear, autonomous ODEs. One ODE for each chemical species. Solution of the jth equation at time t is a real number representing the concentration of species j at time t
An example of modeling a simple substrate-enzyme reaction using Reaction Rate Equations is presented below.
An illustration of the deterministic reaction rate method is a simple substrate-enzyme reaction. This reaction is also known as the Michaelis-Menten enzyme kinetic system. The reaction scheme is as follows:
This scheme describes how enzyme E catylizes a reaction wherein substrate S is converted into product P. This reaction takes place in a fixed volume which is assumed to be well mixed and in thermal equilibrium. This reaction scheme can be separated in the following three reactions:
When substrate S and enzyme E collide complex C is formed which consists of substrate S bounded to enzyme E. Substrate-enzyme complex C can undergo two different reactions. Complex C can disintegrate into enzyme E and substrate S again. Also a reaction can occur wherein bounded substrate S in complex C is reconfigured into product P, where after complex C falls apart in product P and enzyme E. The first reaction is characterized by mass-action rate constant k1 [min-1]. Parameter k1 denotes the reaction speed of the search of the substrate and the enzyme for each other, k2 denotes the failure to reconfigure the substrate after forming complex C, and k3 denotes successful reconfiguration of substrate. The dynamics of the system can be written as the following set of ODEs:
[S], [E], [C] and [P] denote the concentration of substrate S, enzyme E, complex C and product P respectively. Note that enzymes are either unbound represented by symbol E or bound with substrate S to form complex C.
This set of ODEs can be numerically solved, given the initial concentrations and the reaction rates.is an initial value problem where, given the initial concentrations and reaction rates, an ODE simulation. The results for these ODEs are presented in Figure 1. The initial concentrations and reaction rates are used from Chapter 7 of Stochastic Modelling for Systems Biology.
can be conducted.
This simulation can be generated with different tools (MATLAB, Hybrid Chi or COPASI), resulting in the same outcome. The files and code that are needed to simulate the simple substrate-enzyme reaction can be found __here__.
In this example is presented that the response of a deterministic model of the chemical reaction is univocally determined ,given the initial concentrations and reaction rates. Deterministic models are so widespread because usually the biological system under consideration involve a large number of molecules, therefore the average effect is well described through deterministic equations.
The system considered in this example has been described by Leonor Michaelis and Maud Menten for the first time in 1913 (Die Kinetik der Invertinwirkung.). They assumed that the reconfiguration step k3, from C to P and E, was the bottleneck and much slower than k2. Therefore they neglected k3 in some parts of the model. In A note on the kinetics of enzyme action, Briggs and Haldane proposed that the total concentration of the enzyme is much smaller than the substrate concentration. Since this is typically true for biological reactions, this assumption is usually valid. Total enzyme concentration [E]t is equal to the sum of bound and unbound enzyme concentrations:
If [S] » [E]t a steady state is reached in which complex concentration [C] does not change in time:
The enzyme concentration can be substituted by the total enzyme concentration substracted with the complex concentration.
Removing the brackets and rewriting this equation gives:
Dividing with constant k1:
This relation between the constants is the Michaelis-Menten constant KM:
Implementing this constant in the equation gives:
Isolating the complex concentration:
Implementing this equation in the ODE for the product P:
The maximal velocity of the reaction, Vmax, is denoted by the mass-action rate constant of the product forming step, k3, multiplied by the total enzyme concentration [E]t.
If unbound enzyme [E] and complex concentration [C] do not change in time, the decrease in concentration of substrate [S] is equal to the increase in concentration of product [P]. The ODEs for the substrate-enzyme reaction are:
This set of ODEs with the same initial conditions as above can be simulated by using the same tools. The result is presented in Figure 2. The results of the simulation with Michaelis-Menten constant are almost the same to the results of the complete reaction.
Figure 2: Simulation results of the Michaelis Menten reaction
In contrast to a deterministic model, a stochastic model involves random variables. Therefore the substrate-enzyme reaction behavior cannot be predicted a priori, although it can be statistically characterized. Certainly biological reactions fall in the category of stochastic systems, since the very basic steps of every molecular reaction can be described only in terms of its probability of occurrence. Moreover, the diffusion of molecules is a realization of a random walk process (Brownian motion). When the system to be described is based on the interaction of few molecules, or we want to simulate the functioning of a little pool of cells it is necessary to use stochastic models.
There are a few different methods to solve the stochastic system and compute the system evolution with differential equations. We use the stochastic simulation algorithm (SSA), also known as Gillespie's algorithm for the simple substrate-enzyme reaction. In [Higham 2008] the SSA algorithm is explained and a result of the reaction described in the deterministic case is presented in Figure 3.
Figure 3: Simulation of the Stochastic substrate-enzyme model
It can be seen in this example that a stochastic simulation can be performed using differential equations. With this approach only exponential distributions can modeled and simulated. There is also no distinction between the distributions for the search and process times.
The modeling of biological systems with ODEs is also known as conventional mass-action-law-based-theory.
Two fundamental stochastic properties of the law of mass-action are:
When dealing with systems containing non-Markovian elements, results from the ODEs are not satisfying. There is also no distinction between two stochastic independent processes in the ODE model. These processes are the time needed before a substrate and enzyme collide (search time) and the time needed for reconfiguration of the substrate (reconfiguration time). For non-Markovian or stochastic independent biological systems it might be more appropriate to use other approaches. Hereafter we present a modeling approach as is common in modeling manufacturing lines. With these models the througput and flowtime of products can be calculated.