I spend most of my time modelling biochemical reaction networks, these are typically signal transduction systems. As a cell senses some change in its environment it usually want to respond in some way, this could e.g. be the activation of a certain protein. The pathway which moves the information from the outside of the cell (a change in environmental conditions) into the cell and causes a change somewhere (e.g. activating a protein) is usually transduced by a set of biochemical reaction. E.g. the changing environment on the outside might activate a membrane protein (reaction 1), which then may modify a cytoplasmic protein (reaction 2), which then can bind to a third protein (reaction 3), which then act as a transcription faction which activates the production of the response protein (reaction 4).

Here I will just briefly introduce how we model a simple biochemical reaction using mathematics. Given are two biochemical species, X and Y. These two can react to produce a new biochemical species, XY. In addition, XY can dissociate back into its two components, X and Y. The reactions is denoted using the following biochemical equation.

X + Y XY

I have previously described how systems of ordinary differential equations can be used in modelling and here we will use this approach to model our reactions. We want to model how the concentration (in a cell) of our three biochemical species evolves with time. The model will consist of:

1. A time, t. Our model will describe how the system changes with t.
2. A set of three variables, [X], [Y], and [XY]. These describes the concentrations of out three biochemical species. As the concentration changes with time these are all function of t.
3. Some initial conditions; [X]0, [Y]0, and [XY]0. These describes the starting concentration of X, Y, and XY. The model will then simulate how these starting conditions evolve with time.
4. A rule stating how our system (that is, [X], [Y], and [XY]) develops with time.

We will start with (4) since this is the hardest part (but still very simple). Since [X] depends on time it can be described as a function of time, [X](t), were for each value of t we get a corresponding value of [X]. The change with time is simply the derivative of the function [X](t) with respect to t. This is typically denoted as d[X](t)/dt. The rate of change in [X] will depend on the current state of the system, that is the three concentrations [X], [Y] and [XY]. That is, the rate of change in [X] is a function of the variables [X], [Y], and [XY]. We will denote this function f([X]([X],[Y],[XY]). Thus we have the equation d[X](t)/dt = f[X]([X],[Y],[XY]). This is called a differential equation, since it is an equation containing a derivative. Similarly holds for the biochemical species [Y] and [XY] and we can write down:

d[X](t)/dt = f[X]([X],[Y],[XY])
d[Y](t)/dt = f[Y]([X],[Y],[XY])
d[XY](t)/dt = f[XY]([X],[Y],[XY])

this is a system of ordinary differential equations (since it contains several differential equations).

Next we will try to write down expressions for our function f[X], f[Y], and f[XY]. What causes change in the [X]? Two things; [X] is reduced as X and Y combine to form XY (reaction 1), [X] is also increased as XY disassociates back into X and Y (reaction 2). The same goes for [Y] while for [XY] it is the reverse as reaction 2 removes XY while reaction 1 adds XY. Lets call the rate of reaction 1 r1 and of reaction 2 r2. Both these rates are functions of [X], [Y] and [XY]. We can now rewrite our system as:

d[X](t)/dt = r2([X],[Y],[XY]) – r1([X],[Y],[XY])
d[Y](t)/dt = r2([X],[Y],[XY]) – r1([X],[Y],[XY])
d[XY](t)/dt = r1([X],[Y],[XY]) – r2([X],[Y],[XY])

Now we only need to figure out what the rates r1 and r2 are. When does reaction 1 occur? The molecules are randomly diffusing around inside the cell. For the reaction X + YXY to occur a molecule of X and a molecule of Y must bounce into each other. Each time they do bounce into each other there is a certain probability, p1, that there will be a reaction (p1 in turn depends on several different factors). Next we observe that the reaction rate should be linear to both the concentrations [X] and [Y]. That is, if we double the concentration [X] we can expect about double the amount of events were a X and a Y bounce into each other. We will model the rate of such bouncing into each other events as c1·[X]·[Y], were c1 is a constant containing such things as the rate of diffusion etc. (if the molecules are moving around faster in the cell more bouncing events are likely to occur). We will sum all of this up and say r1([X],[Y],[XY]) = k1·[X]·[Y], where k1 is a constant containing both information on how common bouncing events are, as well as how likely a reaction is once such an event occur. Exactly what this constant is composed of is not important, what is interesting is that the reaction rate is some constant times the concentration of the substrates of the reaction.

Next we will look at the rate of the second reaction (this is simpler). Every XY molecule have a certain probability, p2, of dissociating at every time point. If we double the number of XY molecules in the cell we can also expect to see a doubling in the number of dissociation events, hence the reaction rate is linear to the concentration of XY. The concentration of X and Y does not affect the rate of the reaction. We will set r2([X],[Y],[XY]) = k2·[XY].

Finally we can combine all of this information to create the following model

d[X](t)/dtk2·[XY] – k1·[X]·[Y]
d[Y](t)/dt = k2·[XY] – k1·[X]·[Y]
d[XY](t)/dt = k1·[X]·[Y] – k2·[XY]

Now given values of the parameters kand k2, (we call these parameters, as these components of the model, unlike the variables, [X], [Y], and [XY], do not change with time) as well as some initial conditions we can simulate the time development of the system (see previous blog post about simulations of systems of ordinary differential equations).

Lets say that [X]o = 7, [Y]0 = 3, and that [XY]0 = 5. We also set k1 = 1.0 and k= 2.0. We then simulate the system from t = 0 to t = 1. The result can be seen in the following figure: Here we can see that the concentrations of our three biochemical species after a short time reaches some steady states from which the concentrations do not change. It is actually possible to calculate this concentration analytically. One the system have reached the steady state there rate of change in the three concentrations are 0. That is, all the derivatives, d[X](t)/dt, d[Y](t)/dt, and d[XY](t)/dt are all equal to 0. We now have the equation system

0k2·[XY] – k1·[X]·[Y]
0 = k2·[XY] – k1·[X]·[Y]
0 = k1·[X]·[Y] – k2·[XY]

However, all these three equations are identical, so we have

k1·[X]·[Y] = k2·[XY]

which if we insert k1 = 1.0 and k= 2.0 is

[X]·[Y] = 2·[XY]

But we also know that molecules are never created nor destroyed, merely transformed between the two species and Y, and their complex, XY. We can thus calculate the total amount of  and Y in the system from the initial conditions:

[X]tot = [X]+ [XY]o = 7 + 5 = 12
[Y]tot = [Y]+ [XY]o = 3 + 5 = 8

We can then calculate or equilibrium state using these three equations:

[X]·[Y] = 2·[XY]
[X] + [XY] = 12
[Y] + [XY] = 8

Solving this equation system gives us the equilibrium solution

([X]e, [Y]e[XY]e) = (6,2,6)

This can be confirmed by checking the steady state solution of the simulation.