Read the stoichiometric matrix from a tab-delimited data file with M rows and N columns
In [1]:
S <- as.matrix(read.table("shared/ucsd_ecoli_stoich_matrix.txt",
header=TRUE,
row.names=1,
sep="\t",
comment.char="",
quote="",
stringsAsFactors=FALSE))
M <- nrow(S)
N <- ncol(S)
cat(sprintf("The stoichiometrix matrix has %d rows (metabolites) and %d columns (reactions)\n",
M, N))
Read the min and max flux values, and the objective function coefficients, from a tab-delimited data file with N rows and 3 columns
In [2]:
v_minmax <- as.matrix(read.table("shared/ucsd_ecoli_reaction_maxmin.txt",
sep="\t",
header=TRUE,
comment.char="",
row.names=1,
quote="",
stringsAsFactors=FALSE))
stopifnot(all(make.names(rownames(v_minmax)) == colnames(S)))
head(v_minmax)
Create a constraint matrix (with C rows and N columns) and a corresponding constant vector of length C, giving the right-hand-side of the inequality
In [3]:
Svmin <- S %*% matrix(v_minmax[,"min"], ncol=1)
constraint_matrix <- rbind(S,
S,
diag(rep(1,N)),
diag(rep(1,N)))
constraint_rhs <- c(-Svmin,
-Svmin,
rep(0, N),
v_minmax[,"max"] - v_minmax[,"min"])
constraint_dir <- c(rep(">=", M),
rep("<=", M),
rep(">=", N),
rep("<=", N))
stopifnot(length(constraint_rhs) == nrow(constraint_matrix))
stopifnot(ncol(constraint_matrix) == N)
stopifnot(nrow(v_minmax) == N)
dim(constraint_matrix)
Find the flux value that maximizes the "biomass" reaction flux, consistent with the steady-state and thermodynamic constraints
In [5]:
library(lpSolve)
lpres <- lp("max",
objective.in=v_minmax[,"objective"],
const.mat=constraint_matrix,
const.dir=constraint_dir,
const.rhs=constraint_rhs)
print(lpres$solution + v_minmax[,"min"])
In [ ]: