MATLAB approaches to modelling biological systems with differential equations

This site accompanies the tutorial in the Advanced Masters/MSc IT Bioinformatics class on the 24th of February, 2005.

Example 1

The following MATLAB code defines the differential equation describing a simple irreversible monomolecular reaction. Relevant examples in biological experiments are the decay of a radioactive tracer or a fluorescent probe. Also, in first approximation, an equation like this can describe the transition of an activated cellular receptor molecule to its inactive resting state. Save the function in a separate m.file, named decay.m, inside your MATLAB path.

function dydt = decay(t, y)
% A->B
k = 1;
dydt = [-k*y(1)
k*y(1)];

Now type the following instructions at the MATLAB command line to solve the system of differential equations (i.e. simulate the system behavior) and plot the results.

[t y] = ode45(@decay, [0 10], [5 1]);
plot (t, y);
legend ('[A]', '[B]');

Try changing the initial concentrations and the time intervall for which the simulation is performed. What happens to the results when you change the rate constant k in the decay() function?

Example 2

Describing a reversible reaction of a single molecule is only slightly more complex. Now there are two rate constants, one for the forward and one for the backward reaction. Again, safe the function as an m.file, this time named isomerisation.m.

function dydt = isomerisation(t, y)
% A <-> B
k1 = 1;
k2 = 0.5;
dydt = [-k1*y(1)+k2*y(2) % d[A]/dt
k1*y(1)-k2*y(2) % d[B]/dt
];

Simulate and plot the results as before, but referencing the isomerisation() function.

[t y] = ode45(@isomerisation, [0 10], [5 1]);
plot (t, y);
legend ('[A]', '[B]');

How do the results change if you set k1 = 2 and k2 = 1? What stays the same?

Example 3

The following example is a differential equation description of a larger biological system containing 11 different molecular species participating in eleven reactions, some of them involving more than two participants, based on work by Cho et al. (2003). Can you figure out the topology of the reaction network from the equations? Try to sketch a graphical representation of the pathway as a bigraph before you simulate the system.

function dydt = erk_pathway_wolkenhauer(t, y)
% from Kwang-Hyun Cho et al., Mathematical Modeling...
k1 = 0.53;
k2 = 0.0072;
k3 = 0.625;
k4 = 0.00245;
k5 = 0.0315;
k6 = 0.8;
k7 = 0.0075;
k8 = 0.071;
k9 = 0.92;
k10 = 0.00122;
k11 = 0.87;

dydt = [
-k1*y(1)*y(2) + k2*y(3) + k5*y(4)
-k1*y(1)*y(2) + k2*y(3) + k11*y(11)
k1*y(1)*y(2) - k2*y(3) - k3*y(3)*y(9) + k4*y(4)
k3*y(3)*y(9) - k4*y(4) - k5*y(4)
k5*y(4) - k6*y(5)*y(7) + k7*y(8)
k5*y(4) - k9*y(6)*y(10) + k10*y(11)
-k6*y(5)*y(7) + k7*y(8) + k8*y(8)
k6*y(5)*y(7) - k7*y(8) - k8*y(8)
-k3*y(3)*y(9) + k4*y(4) + k8*y(8)
-k9*y(6)*y(10) + k10*y(11) + k11*y(11)
k9*y(6)*y(10) - k10*y(11) - k11*y(11)
];

The simulation uses initial concentrations of the molecular players based on real experimental results from Walter Kolch's lab at the Beatson Institute (Glasgow).

[t y] = ode45(@erk_pathway_wolkenhauer, [0 10], [2.5 2.5 0 0 0 0 2.5 0 2.5 3 0]); %(initial values!)
plot (t, y);
legend ('[Raf1*]', '[RKIP]', '[Raf1/RKIP]', '[RAF/RKIP/ERK]', '[ERK]', '[RKIP-P]', '[MEK-PP]', '[MEK-PP/ERK]', '[ERK-PP]', '[RP]', '[RKIP-P/RP]' );

The ERK pathway simulated in this example is central to many human disorders, including cancer. Having a mathematical description helps biologists understand systems behavior and predict the effect of experimental interventions without having to perform all possible time-consuming experiments. Try to do a knock-out mutation of the tumor suppressor RKIP and examine its effect on system function. In the wet lab, an experiment like this may take several months to perform - how long does it take in MATLAB?

Example 4

As a final example, perform a simple sensitivity analysis of the system. The following function allows one to choose a compound in the ERK pathway, the initial concentration of which will be varied systematically, while tracing the concentration of selected other players.

function [tt,yy] = sensitivity(f, range, initvec, which_stuff_vary, ep, step, which_stuff_show, timeres);

timevec = range(1):timeres:range(2);

vec = [initvec];
[tt y] = ode45(f, timevec, vec);
yy = y(:,which_stuff_show);

for i=initvec(which_stuff_vary)+step:step:ep;
vec(which_stuff_vary) = i;
[t y] = ode45(f, timevec, vec);
tt = [t];
yy = [yy y(:,which_stuff_show)];
end

Call the function in the same general way as before, defining a time range to analyse ([0 1]), the initial concentrations, the compound to vary (number 5, the ERK molecule) and the final concentration (6) and stepsize (1) of the range to be examined. The concentration of compound 8 (the MEK-PP/ERK complex) will be traced in steps of 0.05 secs for the selected time range and the results are plotted as a surface graph.

[t y] = sensitivity(@erk_pathway_wolkenhauer, [0 1], [2.5 2.5 0 0 0 0 2.5 0 2.5 3 0], 5, 6, 1, 8, 0.05);
surf (y);

Try simulating larger time intervals or tracing other compounds. Where do qualitative changes in systems behavior occur that could indicate new therapeutic options?


Rainer Breitling, Bioinformatics Research Centre, Glasgow, February 2005.