diff --git a/networkComplexityBigg/qc/iSMenergy.R b/networkComplexityBigg/qc/iSMenergy.R
new file mode 100644
index 0000000000000000000000000000000000000000..3e22ff0f5ab7f684baef9e148b7adc0ce661449d
--- /dev/null
+++ b/networkComplexityBigg/qc/iSMenergy.R
@@ -0,0 +1,73 @@
+#!/usr/bin/Rscript
+
+library(sybil)
+
+
+versions <- c("ver1", "ver1.1", "ver1.2", "ver2")
+names(versions) <- versions
+
+data <- lapply(versions, function(v){
+	print(load(paste0("../models/universalBiGG.", v, ".Rdata")))
+	return(list(map=modelReactMap, uni=get(grep("^uni", ls(), value=T))))
+})
+
+data <- lapply(data, function(d){
+	within(d, {
+				ism <- rmReact(uni, setdiff(react_id(uni), map$iSM199))
+				ism <- changeObjFunc(ism, "ENERGY_BOF")
+				})
+})
+
+print(sapply(lapply(data, "[[", "ism"), optimizeProb))
+
+stopifnot(setequal(met_id(data$ver1$ism), met_id(data$ver1.1$ism)))
+
+s1 <- S(data$ver1$ism)
+s2 <- S(data$ver1.1$ism)
+
+s2 <- s2[match(met_id(data$ver1$ism), met_id(data$ver1.1$ism)),]
+
+
+matching <- sapply(1:ncol(s1), function(i){
+	which(sapply(1:ncol(s2), function(j){
+		c1 <- s1[,i]
+		c2 <- s2[,j]
+		if(sum(c1!=0) != sum(c2!=0)){
+			return(F)
+		}
+		return(all(c1 == c2))
+	}))
+})
+
+
+matchingNames <- lapply(matching, function(x) react_id(data$ver1.1$ism)[x])
+names(matchingNames) <- react_id(data$ver1$ism)
+
+for(i in seq(along=matching)){
+	if(any(lowbnd(data$ver1$ism)[i] != lowbnd(data$ver1.1$ism)[matching[[i]]])){
+		print(matchingNames[i])
+		print(lowbnd(data$ver1$ism)[i])
+		print(lowbnd(data$ver1.1$ism)[matching[[i]]])
+	}
+	if(any(uppbnd(data$ver1$ism)[i] != uppbnd(data$ver1.1$ism)[matching[[i]]])){
+		print(matchingNames[i])
+		print(uppbnd(data$ver1$ism)[i])
+		print(uppbnd(data$ver1.1$ism)[matching[[i]]])
+	}
+	
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+