--TUTORIAL:
-- Using computer algebra software to analyze reaction networks
--Anne Shiu
--AIM workshop on "Mathematical problems arising from biochemical reaction networks"
--Wednesday, March 27, 2013
--(This file will be available on my website.)
--OUTLINE:
-- Basics of Macaulay2
-- Example 1: steady states
-- Example 2: computation by networks
-- Example 3: siphons (with application to Global Attractor Conjecture)
--Mathematical background on computational algebraic geometry:
-- 1. Ideals, Varieties, and Algorithms, third edition (Cox, Little, O'Shea), 2007.
-- 2. In the context of reaction networks:
-- Part A of "Bistability and oscillations in chemical reaction networks",
-- (Mirela Domijan, Markus Kirkilionis),
-- Journal of Mathematical Biology, October 2009, Volume 59, Issue 4, pp 467-501
--Computer algebra software:
-- 1. Macaulay2 (today's focus)
-- 2. Singular
-- 3. Magma, CoCoA, Mathematica, Maple, ...
--References for Macaulay2:
-- 1. Installation and documentation:
-- http://www.math.uiuc.edu/Macaulay2/
-- 2. "Ideals, Varieties, and Macaulay 2" (Bernd Sturmfels)
-- in book "Computations in algebraic geometry with Macaulay 2"
-- 3. "Using Macaulay2" (D. Eisenbud, D.R. Grayson and M.E. Stillman)
--USAGE:
-- in terminal OR
-- via emacs (fn-F12 to start Macaulay2, fn-F11 to send line)
----------------------------------------------------------------
printWidth = 60
--Simple arithmetic
-- (note: "QQ" tells us that the output is a rational number)
3/5 + 7/11
--End a line with a semi-colon to print no output:
2*3;
--Access the previous output using "oo":
oo+1
----------------------------------------------------------------
--EXAMPLE ONE:
-- glycolysis network
--GOAL:
-- analyze steady states
--3 NETWORK SPECIES:
-- X_1 = fructose-1,6-biphosphate
-- X_2 = fructose-6-phosphate
-- X_3 = intermediate species that reacts reversibly with X_2
--NETWORK has 7 reactions:
-- 2 X_1 + X_2 --> 3 X_1
-- X_2 <--> 0 <--> X_1 <--> X_3
--Reference: Network related to network of (Sel'kov 1968) in
-- "Deformed Toric Ideal Constraints on Stoichiometric Networks"
-- (Masamichi Sato, Kenji Fukumizu), arXiv:1202.5092
--Treat reaction rate constants as unknowns:
Coeffs = QQ[k_21, k_34, k_43, k_46, k_64, k_56, k_65]
--Allow rational functions in the rate constants:
CoeffsField = frac Coeffs
--example of such a rational function:
k_21^2/(k_65-k_43)
--Work polynomial ring over the field defined above,
-- with one variable for each species:
R = CoeffsField[x_1..x_3]
--example of a polynomial in this polynomial ring:
k_21^2/(k_65-k_43) * x_2^2 - x_3
--The ODEs:
g_1 = k_21 * x_1^2*x_2 + k_46 - k_64 * x_1 - k_34 * x_1 + k_43 * x_3
g_2 = -k_21 * x_1^2*x_2 + k_56 - k_65 * x_2
g_3 = k_34 * x_1 - k_43 * x_3
--Define the ideal generated by the ODE's
-- and the conservation laws, that is,
-- the set of a polynomial-combinations of
-- the ODEs:
ODEsIdeal = ideal(g_1,g_2,g_3)
--A Grobner basis of the ideal is a "nice" set of generators
-- of the ideal, that is,
-- a (hopefully) easier system of equations to solve
-- to obtain the steady states:
ODEsGB = gb ODEsIdeal
generators ODEsGB
expression generators gb ODEsIdeal
-- See:
-- basis element #1 is linear in x_2 and x_3
-- basis element #2 is linear in x_1 and x_3
-- basis element #3 is a degree-3 polynomial in only x_3,
-- and by inspection, has 3 sign changes in the list of
-- coefficients, so by Descartes' Rule of Signs, there
-- are 1 or 3 positive roots
-- To do: further analysis of the coefficients to
-- determine the number of steady states.
--Note: this Grobner basis is with respect to the default
-- monomial ordering, graded reverse lexicographic order)
--Note: this network is not weakly reversible, so
-- deficiency theory does not apply.
----------------------------------------------------------------
--Exercise -- SOME NEXT STEPS:
-- 1. Another (computer algebra) approach to analyzing steady states:
-- "stoichiometric network analysis" / "flux-balance analysis"
-- (see Example 3 in Domijan and Kirkilionis 2009)
-- 2. Sufficient conditions for bistability/oscillations:
-- (see examples in Sections 4-5 of Domijan and Kirkilionis 2009)
----------------------------------------------------------------
--EXAMPLE TWO:
-- 2X --> Y
--GOAL: What does this CRN compute? (Doty & Soloveichik)
-- (assuming we start with initial conditions x = xInit, y = 0)
--Treat reaction rate constant and conservation value as unknowns:
Coeffs' = QQ[k, xInit];
CoeffsField' = frac Coeffs';
--Define ring, and impose elimination order (to eliminate x)
R' = CoeffsField'[x,y, MonomialOrder => Eliminate 1]
--The ODE:
ODE_1 = -2*k*x^2;
ODE_2 = k*x^2;
--Define the ideal generated by the ODE's
-- and the conservation laws:
ODEsIdeal' = ideal(ODE_1,ODE_2, 2*y+x-xInit);
--A Grobner basis of the ideal is a "nice" set of generators
-- of the ideal, that is,
-- a (hopefully) easier system of equations to solve
-- to obtain the steady states:
ODEsGB' = gb ODEsIdeal'
G = generators ODEsGB'
expression generators gb ODEsIdeal'
--Elimination ideal is:
ideal selectInSubring(1, G)
-- See: this is (y - xInit/2)^2
-- Conclusion: the function is y = f(x) = xInit/2.
-- Question: How to use computer algebra to determine that
-- the following network encodes f(x_1,x_2) = x_1(0) - 2*x_2(0) ?
-- X1 --> Y
-- X2 --> 2K
-- Y+K --> 0
----------------------------------------------------------------
--Exercise -- SOME NEXT STEPS using ELIMINATION:
-- 1. steady-state invariants (and Absolute Concentration Robustness)
-- 2. QSSA
----------------------------------------------------------------
--EXAMPLE THREE:
-- Receptor-ligand dimer model network
--GOAL: Analyze boundary steady states (with application to
-- Global Attractor Conjecture)
--TECHNIQUE: use binomial primary decomposition;
-- see Shiu and Sturmfels, "Siphons in chemical reaction networks"
-- Bulletin of Mathematical Biology (2010).
--Biochemical species are as follows:
-- A = receptor
-- B = ``dimer'' state of A
-- C = ligand that can bind either to A (to form D) or to B (to form E)
--"Square" Network:
-- 2A+C <--> A+D <--> E <--> B+C <--> 2A+C
--Reference:
-- PhD thesis of M. Chavez (2003), Section 7.2
-- see also: D. Anderson (2008), Ex. 4.1, SIAM J. Appl. Math.
-- (def) a "siphon" in a chemical reaction network is a subset of
-- the species whose absence is forward-invariant with respect to
-- the dynamics.
-- Application: only faces defined by siphons (zero-coordinate set =
-- a siphon) may contain (boundary) steady states
--According to Theorem 3.1 of Shiu & Sturmfels (2010),
-- the minimal siphons can be computed as follows:
R1 = QQ[A,B,C,D,E]
-- This ideal contains all the monomials:
M = ideal(A^2*C, A*D, E, B*C)
--(Note: for networks not strongly connected, M is replaced by
-- a certain binomial ideal.)
-- Here the command decompose outputs the minimal primes, which
-- correspond to the minimal siphons:
decompose(M)
--(Note: in general, the minimal siphons are the inclusion-minimal
-- subsets of species contained in minimal primes.)
--Next goal:
-- which siphons actually define faces of a compatibility class?
--There are 3 types of compatibility classes. Let's look at the type
-- of chamber Omega(1) (Figure 1 in Shiu & Sturmfels), represented by
-- the point c0= (0,0,1,1,0).
--According to Algorithm 3.6 (Shiu & Sturmfels), we can determine the
-- faces as follows:
c0 = {0,0,1,1,0};
R = QQ[A,B,C,D,E, Weights => c0];
--Ideal containing binomials representing the reactions and the product of species:
IG = ideal(A^2*C-A*D, A*D-E, E-B*C, A*B*C*D*E);
--Ideal representing conservation relations:
ICons = ideal(C*D*E-1, A*B^2*D*E^2-1);
--Here the output is the unique minimal siphon for type-1 compatibility classes:
Bc0 = dual radical monomialIdeal leadTerm ICons;
decompose saturate(IG,Bc0)
-- Conclusion: only one face, the one with E, B, A = 0, may admit a (boundary)
-- steady state in its interior. (It is a vertex.)
--Application to the Global Attractor Conjecture:
-- Determine whether a network admits any boundary steady states,
-- and if so, which faces of a stoichiometric class they lie in.
-- Example: in weakly reversible, deficiency zero network above,
-- GAC holds because only vertices contain steady states.
----------------------------------------------------------------
-- if time permits: motivate Dan's tutorial by LRAT example
----------------------------------------------------------------
exit