
Using Systems Biology to Analyze ReactionsFrom MathBioContentsIntroductionSystems biology is a quickly emerging approach in biological research. Technological advances in highthroughput technologies are enabling the biological community to define whole catalogs of the molecular parts of cellular components. Thus, it is possible to reconstruct the networks of biochemical reactions at the genomescale: this is the essence of systems biology. To familiarize myself with modern mathematical biology, I studied Bernhard Palsson's texts on systems biology: Simulation of Dynamic Network States and Properties of Reconstructed Networks. In particular, I focused on the stoichiometric matrix, which forms the mathematical crux of systems biology, as because functions of reconstructed networks are defined by the connections of their parts. I then applied his methods to analyzing the general behavior of fundamental chemical reactions like the reversible reaction and bilinear reversible reaction. The Stoichiometric MatrixThe stoichiometric matrix represents the reaction topology of a network. Mathematically speaking, it represents a linear transformation of the flux vector to the time derivatives of the concentration vector .
Typically, the stoichiometric matrix is sparse, and its nonzero elements have either a value of 1 or 1. Additionally, the matrix is considered "knowable", as it is composed of integers with no associated error. For more information on the stoichiometric matrix and its singular value decomposition, go here. Reversible Reaction
We can write the stoichiometric matrix as , where and is formulated using the Law of Mass Action, which is more or less based on the principle that a reaction's speed is proportional to the concentrations of its reactants, because molecular collisions fuel the reaction itself. Thus, we have two differential equations:
Also useful in the analysis of this basic reaction is the spanning vector of the left null space, whose dimension is 1 according to the Fundamental Theorem of Linear Algebra: . As presented in Palsson's texts, the left null space gives conservation quantities  that is, combinations of reactants whose combined concentration does not change over time. The spanning vector of the left null space can be written as Interpreted biologically, gives a conservation quantity, which should not be surprising, as this basic reversible reaction is closed. The existence of a conservation quantity is illustrated in a phase portrait between the two reactants, seen at right.
MATLAB Code
function [xprime] = system_reversible(t,x,k1,kneg1) xprime(1,1) = k1*x(1) + kneg1*x(2); xprime(2,1)= k1*x(1)  kneg1*x(2);
% File: reversible_solver.m % Author: Sam Hsu % Date: 5.25.2012 % Set the parameters: the rate constants for the forward and reverse % reaction. k1 = 1; kneg1 = 2; % Defines the equilibrium constant. K = k1/kneg1; % Set the time scale to graph the concentrations and the initial % conditions for x1 and x2. timeScale = [0,2]; initialConditions = [1,0]; % Numerically solves the differential equations system. [t,x] = ode45(@system_reversible, timeScale, initialConditions, [], k1, kneg1); % Defines pooled variables. p1 is a disequilibrium quantity that % thermodynamically measures how far the reaction is from equilibrium; p2 % is a conservation quantity because it spans the left null space of the % stoichiometric matrix. p1 = x(:,1)  x(:,2)/K; p2 = x(:,1) + x(:,2); % Plots the solutions. figure % Time profile of x1 and x2. hold on plot(t,x(:,1),'k'); plot(t,x(:,2),'r'); hold off xlabel('Time') ylabel('Concentration') title('Time Profile of a Reversible Reaction') legend('x1', 'x2') figure % Phase portrait of x1 versus x2. plot(x(:,1),x(:,2)); xlabel('x1') ylabel('x2') title('Phase Portrait of a Reversible Reaction') figure % Time profile of pooled variables p1 and p2. hold on plot(t,p1,'g'); plot(t,p2,'c'); hold off xlabel('Time') ylabel('Concentration') title('Time Profile of Pooled Variables') legend('p1 = x1  x2/K', 'p2 = x1 + x2') Reversible Bilinear Reaction
We can write the stoichiometric matrix as , where and , again formulated using the Law of Mass Action. Thus, we have two differential equations:
Solving for the left null space gives two linearly independent spanning vectors in and , corresponding to the conservation quantities and , as expected. Another useful variable in the analysis of this basic reaction is the disequilibrium quantity, which measures the distance of the overall reaction from chemical equilibrium, also presented in Palsson's texts. This quantity can be computed by examining the net flux through the reaction, that is: which can be rewritten as:
MATLAB Code
function [xprime] = system_bilinear_reversible(t,x,k1,kneg1) xprime(1,1) = k1*x(1)*x(2) + kneg1*x(3); xprime(2,1) = k1*x(1)*x(2) + kneg1*x(3); xprime(3,1) = k1*x(1)*x(2)  kneg1*x(3);
% File: bilinear_reversible_solver.m % Author: Sam Hsu % Date: 5.27.2012 % Set the parameters: the rate constants for the forward and reverse % reaction. k1 = 1; kneg1 = 1; % Defines the equilibrium constant. K = k1/kneg1; % Set the time scale to graph the concentrations and the initial % conditions for x1, x2, and x3. timeScale = [0,2]; initialConditions = [3,2,0]; % Numerically solves the differential equations system. [t,x] = ode45(@system_bilinear_reversible, timeScale, initialConditions, [], k1, kneg1); % Defines pooled variables. p1 is a disequilibrium quantity that % thermodynamically measures how far the reaction is from equilibrium; p2 % and p3 are conservation quantities because they span the left null space % of the stoichiometric matrix. p1 = x(:,1).*x(:,2)  x(:,3)/K; p2 = x(:,1) + x(:,3); p3 = x(:,2) + x(:,3); % Plots the solutions. figure % Time profile of x1, x2, and x3. hold on plot(t,x(:,1),'k'); plot(t,x(:,2),'r'); plot(t,x(:,3),'b'); hold off xlabel('Time') ylabel('Concentration') title('Time Profile of a Reversible Bilinear Reaction') legend('x1', 'x2', 'x3') figure % Phase portraits of x1 versus x3. plot(x(:,1),x(:,3)); xlabel('x1') ylabel('x3') title('Phase Portrait of One Reactant and One Product in a Reversible Bilinear Reaction') figure % Time profile of pooled variables p1 and p2. hold on plot(t,p1,'g'); plot(t,p2,'c'); plot(t,p3,'b'); hold off xlabel('Time') ylabel('Concentration') title('Time Profile of Pooled Variables in a Reversible Bilinear Reaction') legend('p1 = x1*x2  x3/K', 'p2 = x1 + x3', 'p3 = x2 + x3') 