September 22, 2017, Friday

# Using Systems Biology to Analyze Reactions

## Introduction

Systems biology is a quickly emerging approach in biological research. Technological advances in high-throughput 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 genome-scale: 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 Matrix

The stoichiometric matrix ${\textbf {S}}$ represents the reaction topology of a network. Mathematically speaking, it represents a linear transformation of the flux vector ${\textbf {v}}$ to the time derivatives of the concentration vector ${\frac {d{\textbf {x}}}{dt}}$.

${\textbf {Sv}}={\frac {d{\textbf {x}}}{dt}}.$

Thus, rows of the stoichiometric matrix represent reactions, while columns represent metabolites in the network.

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

$x_{1}{\xleftarrow {v_{{-1}}}}{\xrightarrow {v_{1}}}x_{2}$

Time profile of the reversible reaction. Conditions used: time scale = 0 - 2; k1 = 1, k−1 = 2; x1(0) = 1, x2(0) = 0

We can write the stoichiometric matrix as ${\textbf {S}}={\begin{pmatrix}-1&1\\1&-1\end{pmatrix}}$, where ${\textbf {v}}={\begin{pmatrix}k_{1}x_{1}\\k_{{-1}}x_{2}\end{pmatrix}}$ and ${\textbf {x}}={\begin{pmatrix}x_{1}&x_{2}\end{pmatrix}}.$ ${\textbf {v}}$ 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:

${\begin{cases}{\frac {dx_{1}}{dt}}=-k_{1}x_{1}+k_{{-1}}x_{2}\\{\frac {dx_{2}}{dt}}=k_{1}x_{1}-k_{{-1}}x_{2}\end{cases}}$

Conservation Quantity

The linear relationship between x1 and x2 demonstrates the existence of a conservation relationship between the two reactants.

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: $rank({\textbf {S}})=1$. 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 ${\begin{pmatrix}1&1\end{pmatrix}}.$

Interpreted biologically, $x_{1}+x_{2}$ 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

system_reversible.m contains the reversible reaction differential equations to be simulated.

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);


reversible_solver.m will solve the reversible system and generate three plots: a time profile of the reactants, a phase portrait of x1 versus x2, and a time profile of the pooled variables. The rate constants, time scale, and initial conditions can easily be changed in this file.

% 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

$x_{1}+x_{2}{\xleftarrow {v_{{-1}}}}{\xrightarrow {v_{1}}}x_{3}$

Time profile of the bilinear reversible reaction. Conditions used: time scale = 0 - 2; k1 = 1, k−1 = 1; x1(0) = 3, x2(0) = 2, x3(0) = 0

We can write the stoichiometric matrix as ${\textbf {S}}={\begin{pmatrix}-1&1\\-1&1\\1&-1\end{pmatrix}}$, where ${\textbf {v}}={\begin{pmatrix}k_{1}x_{1}x_{2}\\k_{{-1}}x_{3}\end{pmatrix}}$ and ${\textbf {x}}={\begin{pmatrix}x_{1}&x_{2}&x_{3}\end{pmatrix}}.$ ${\textbf {v}}$, again formulated using the Law of Mass Action. Thus, we have two differential equations:

${\begin{cases}{\frac {dx_{1}}{dt}}=-k_{1}x_{1}x_{2}+k_{{-1}}x_{3}\\{\frac {dx_{2}}{dt}}=-k_{1}x_{1}x_{2}+k_{{-1}}x_{3}\\{\frac {dx_{3}}{dt}}=k_{1}x_{1}x_{2}-k_{{-1}}x_{3}\end{cases}}$

Conservation and Disequilibrium Quantities

Solving for the left null space gives two linearly independent spanning vectors in ${\begin{pmatrix}1&0&1\end{pmatrix}}$ and ${\begin{pmatrix}0&1&1\end{pmatrix}}$, corresponding to the conservation quantities $x_{1}+x_{3}$ and $x_{2}+x_{3}$, 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:

$v_{{1,net}}=k_{1}x_{1}x_{2}-k_{{-1}}x_{3}$

which can be rewritten as:

Behavior of the conservation and disequlibrium quantities of the bilinear reversible reaction over time.
$v_{{1,net}}=k_{1}(x_{1}x_{2}-{\frac {x_{3}}{K}}),$
where $K={\frac {k_{1}}{k_{{-1}}}},$ the equilibrium constant. In this expression, $k_{1}$ represents the kinetics of the system, whereas $x_{1}x_{2}-{\frac {x_{3}}{K}}$ represents the thermodynamics of the system and can be taken to measure the distance from equilibrium.

The behavior of the conservation and disequlibrium quantities over time can be seen in the figure to the left. Observe that the conservation quantities remain unchanged over time, whereas the disequilibrium quantity decays to zero as the reaction progresses. These quantities, together, are sometimes known as pooled variables.

MATLAB Code

system_bilinear_reversible.m contains the bilinear reversible reaction differential equations to be simulated.

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);


bilinear_reversible_solver.m will solve the bilinear reversible system and generate three plots: a time profile of the reactants, a phase portrait of x1 versus x3, and a time profile of the pooled variables. The rate constants, time scale, and initial conditions can easily be changed in this file.

% 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')