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:

- A time,
*t*. Our model will describe how the system changes with*t*. - 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*. - 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. - 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 *r _{1}* and of reaction 2

*r*. Both these rates are functions of [

_{2}*X*], [

*Y*] and [

*XY*]. We can now rewrite our system as:

*d*[*X*](*t*)/*dt* = *r _{2}*([

*X*],[

*Y*],[

*XY*]) –

*r*([

_{1}*X*],[

*Y*],[

*XY*])

*d*[

*Y*](

*t*)/

*dt*=

*r*([

_{2}*X*],[

*Y*],[

*XY*]) –

*r*([

_{1}*X*],[

*Y*],[

*XY*])

*d*[

*XY*](

*t*)/

*dt*=

*r*([

_{1}*X*],[

*Y*],[

*XY*]) –

*r*([

_{2}*X*],[

*Y*],[

*XY*])

Now we only need to figure out what the rates *r _{1}* and

*r*are. When does reaction 1 occur? The molecules are randomly diffusing around inside the cell. For the reaction

_{2}*X*+

*Y*→

*XY*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,

*p*, that there will be a reaction (

_{1}*p*in turn depends on several different factors). Next we observe that the reaction rate should be linear to both the concentrations [

_{1}*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

*c*·[

_{1}*X*]·[

*Y*], were

*c*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

_{1}*r*([

_{1}*X*],[

*Y*],[

*XY*]) =

*k*·[

_{1}*X*]·[

*Y*], where

*k*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.

_{1}Next we will look at the rate of the second reaction (this is simpler). Every *XY* molecule have a certain probability, *p _{2}*, 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

*r*([

_{2}*X*],[

*Y*],[

*XY*]) =

*k*·[

_{2}*XY*].

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

*d*[*X*](*t*)/*dt* = *k _{2}*·[

*XY*] –

*k*·[

_{1}*X*]·[

*Y*]

*d*[

*Y*](

*t*)/

*dt*=

*k*·[

_{2}*XY*] –

*k*·[

_{1}*X*]·[

*Y*]

*d*[

*XY*](

*t*)/

*dt*=

*k*·[

_{1}*X*]·[

*Y*] –

*k*·[

_{2}*XY*]

Now given values of the parameters *k _{1 }*and

*k*, (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).

_{2}Lets say that [*X*]_{o} = *7*, [Y]_{0} = *3*, and that [*XY*]_{0} = *5*. We also set *k _{1}* =

*1.0*and

*k*=

_{2 }*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

*0* = *k _{2}*·[

*XY*] –

*k*·[

_{1}*X*]·[

*Y*]

*0*=

*k*·[

_{2}*XY*] –

*k*·[

_{1}*X*]·[

*Y*]

*0*=

*k*·[

_{1}*X*]·[

*Y*] –

*k*·[

_{2}*XY*]

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

*k _{1}*·[

*X*]·[

*Y*] =

*k*·[

_{2}*XY*]

which if we insert *k _{1}* =

*1.0*and

*k*=

_{2 }*2.0*is

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

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

*[X] _{tot} = [X]_{o }+ [XY]_{o} = 7 + 5 = 12*

*[Y]*

_{tot}= [Y]_{o }+ [XY]_{o}= 3 + 5 = 8We can then calculate or equilibrium state using these three equations:

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

*[X] _{ }+ [XY] = 12
*

*[Y]*

_{ }+ [XY] = 8Solving 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.