Team:TU Delft-Leiden/Modeling

Introduction to deterministic modeling

The three main systems that form a part of our iGEM project, the landmine promoter, the formation of (conductive) curli, and the assembly of the Extracellular Electron transport (EET) pathway, all involve various biological mechanisms that are not fully understood yet. To gain more insight in these mechanisms, we will apply the technique of deterministic modeling. Deterministic modeling, as opposed to stochastic modeling, does not involve randomness, and therefore yields “exact” solutions. In this section, a general outline of the strategy used to set up and analyze a deterministic model will be provided.

The first step in deterministic modeling should be to get a clear view of what you want to model. In all our cases, we want to obtain a model that predicts the amount of a certain protein or protein complex formed at a certain time. Once this modeling goal is established, you should get a very good and detailed understanding of the system at hand. What compounds are involved? How do the different compounds react with each other? Is localization and transport an important part of the system? Answering these questions (and many more) will mostly be achieved by doing extensive literature research.

Once a clear overview of the model is established, you should decide which processes are most important. Although it is tempting to describe your system in as much detail as possible, this will make your model cluttered and difficult to work with and analyze. Besides that, every reaction you include introduces at least one extra parameter, such as a reaction rate. To find exact values for those parameters is nigh-impossible, and even making an educated guess is not as easy as it seems. To decide whether a process is important or not, a good strategy is to see if certain steps are described as rate-limiting. Rate-limiting steps are steps that take a lot of time and therefore have a big influence on the time behavior of your system. Another class of important processes are processes in which new compounds are formed. However, not every reaction step in such a process should be described as a separate reaction. For example if you have A turns to B, and B turns to C, and C turns to D (A -> B -> C -> D), this can be summarized as A -> D immediately. The reaction rate of this summarized reaction can be estimated to be the lowest rate in the process (the rate-limiting step).

Once you have decided which processes are important enough to include in your model, you should write down reaction equations for the processes. An generic example of a reaction equation looks like this:

$A\; +\; B\; ->k->\; C$ (1)

This reaction equation tells us that compound A and B react at a reaction rate k to form compound C. Writing down such equations from the information you found during research is the most important step of modeling. This is where you write down what you think is an accurate description of how the system works in real life.

After establishing a system of reaction equations, the next step is to convert this system to a system of coupled Ordinary Differential Equations (ODEs). ODEs are a broad class of differential equations, which have in common that they contain a function of one dependent unknown variable and its derivatives. In the case of the deterministic modeling of biological systems, this means that for every compound (the dependent variable), you write down how its concentration changes in time (the time derivative of the dependent variable). To make this more clear, we will write down the system of ODEs describing reaction (1) step-by-step.

We will first consider the concentration of the compound A ([A]) as the dependent variable. A reacts with B to form C; this means A will be removed from the system due to this reaction. To have this reaction occurring, at least one A and one B is needed. The rate at which this reaction occurs is therefore proportional to both the concentration of A and B, and of course the reaction rate k. Written down as an ODE, you have this:

$d/dt[A]=\; -k[A][B]$

For compound B, exactly the same happens, so the differential equation is similar:

$d/dt[B]=\; -k[A][B]$

For the increase of [C], it can easily be seen that this equals the decrease of A (or B), and therefore the ODE will read:

$d/dt[C]=\; k[A][B]$

The system of ODEs we have arrived at is of course extremely simple. In real deterministic modeling of biological systems, usually a lot more compounds are involved which might react in different ways, yielding more and longer ODEs. It is important to keep track of all compounds and to make sure that you have a closed system.

When you have written down a system of ODEs, there are a couple of different strategies you can pursue. Perhaps the easiest is to make a Matlab script containing your system of ODEs and using the function ode45. This solves the system in an iterative way, and will give you the concentration of all compounds as a function of time. Although this is easy, it can be quite computationally intensive. A more elegant way, which can be applied to small systems, is to solve the system by hand, to obtain a closed form solution. Depending on your system and your math skill, this can be easy, hard, or impossible. A third option would be to search for steady state solutions. Steady state means that the system does not change anymore, i.e. all time derivatives are zero. Finding steady state solutions is rather easy; however, it is not possible to find a meaningful steady-state solution for a lot of systems. The generic system described above has a steady state solution, but the only thing it will tell you is that nothing happens when [A] or [B] is zero.

Once you have found a solution to your system, you would like it to show (approximately) the same results as found during the lab work or in literature. To realize this, you need to choose your unknown parameters in such a way that your modeling results match the data. This is called fitting. When you have found an analytical solution to your system, you can use a range of Matlab functions to do this, such as nlinfit. If your system is more complex, or if you want to fit data depending on something else than the independent variable (time in most cases), this will usually not work. Although this sounds quite non-scientific, the best approach in such a case is to guess your parameters and adjust them in a iterative fashion until you have found a fit that matches your data. Following this approach, you will most probably not be able to determine the exact value of a parameter, but obtaining the order of magnitude of a parameter or the ratio between different parameters will nevertheless give you valuable insight in the system.