diff --git a/DESCRIPTION b/DESCRIPTION
index bc7a4d2cf6e6ef3da728e4f944db641d2809569d..3d90c4d2362e59cbb820b0e790dbc304ce0c9f03 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sybil
 Type: Package
 Title: Efficient Constrained Based Modelling in R
-Version: 2.0.3
-Date: 2017-04-20
+Version: 2.0.4
+Date: 2017-07-21
 Authors@R: c(
 	     person(c("C.", "Jonathan"), "Fritzemeier", role = c("cre", "ctb"), email = "clausjonathan.fritzemeier@uni-duesseldorf.de"),
 	     person("Gabriel", "Gelius-Dietrich", role = c("aut")),
@@ -21,7 +21,8 @@ URL:
 Description: This Systems Biology Library for R implements algorithms for constraint based analyses of metabolic networks (e.g. flux-balance analysis (FBA), minimization of metabolic adjustment (MOMA), regulatory on/off minimization (ROOM), robustness analysis and flux variability analysis). Most of the current LP/MILP solvers are supported via additional packages.
 LazyLoad: yes
 License: GPL-3 | file LICENSE
-Collate: generics.R validmodelorg.R validoptsol.R validreactId.R
+Collate: generics.R validmodelorg.R validoptsol.R validreactId.R validreact.R
+        reactClass.R
         validreactId_Exch.R validsysBiolAlg.R addAlgorithm.R
         addExchReact.R addReact.R addSolver.R blockedReact.R
         bracket_pairs.R ceilValues.R changeBounds.R changeGPR.R
@@ -53,7 +54,7 @@ Collate: generics.R validmodelorg.R validoptsol.R validreactId.R
         sysBiolAlg_lmomaClass.R sysBiolAlg_momaClass.R
         sysBiolAlg_mtfClass.R sysBiolAlg_mtfEasyConstraintClass.R
         sysBiolAlg_roomClass.R sybilLogClass.R upgradeModelorg.R
-Packaged: 2017-04-20 12:34:14 UTC; jonathan
+Packaged: 2017-07-21 12:34:14 UTC; jonathan
 Author: C. Jonathan Fritzemeier [cre, ctb],
   Gabriel Gelius-Dietrich [aut],
   Rajen Piernikarczyk [ctb],
diff --git a/NAMESPACE b/NAMESPACE
index e782ef4835465a07be574fe5c2c7610ae3cb85da..faaef08ab95c9700db55d457e6fe661d1a6090d1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -30,6 +30,7 @@ optsol_robAna,
 optsol,
 pointerToProb,
 ppProc,
+react,
 reactId,
 sybilError,
 sybilLog
diff --git a/R/generics.R b/R/generics.R
index af64b93398d579728264339f40fbcd40d648d6bd..9ba875409869a84ed6f650ef0c00d2e17d50302e 100644
--- a/R/generics.R
+++ b/R/generics.R
@@ -351,6 +351,10 @@ setGeneric(name = "getObjVal",
            def  = function(lp) { standardGeneric("getObjVal") }
 )
 
+setGeneric(name = "getReaction",
+           def  = function(X, ...) { standardGeneric("getReaction") }
+)
+
 setGeneric(name = "getRedCosts",
            def  = function(lp) { standardGeneric("getRedCosts") }
 )
@@ -389,6 +393,13 @@ setGeneric(name = "gprRules<-",
            def  = function(object, value) { standardGeneric("gprRules<-") }
 )
 
+setGeneric(name = "gprRule",
+           def  = function(object) { standardGeneric("gprRule") }
+)
+setGeneric(name = "gprRule<-",
+           def  = function(object, value) { standardGeneric("gprRule<-") }
+)
+
 setGeneric(name = "hasEffect",
            def  = function(object) { standardGeneric("hasEffect") }
 )
@@ -869,6 +880,13 @@ setGeneric(name = "S<-",
            def  = function(object, value) { standardGeneric("S<-") }
 )
 
+setGeneric(name = "s",
+           def  = function(object) { standardGeneric("s") }
+)
+setGeneric(name = "s<-",
+           def  = function(object, value) { standardGeneric("s<-") }
+)
+
 setGeneric(name = "scaleProb",
            def  = function(lp, ...) { standardGeneric("scaleProb") }
 )
diff --git a/R/modelorgClass.R b/R/modelorgClass.R
index 5cd4063e5cd24ab01b9f406e35558bea526576d4..a8da37cc8627a08dafdc603d0bacd476722dbb14 100644
--- a/R/modelorgClass.R
+++ b/R/modelorgClass.R
@@ -965,7 +965,44 @@ setMethod("printMetabolite", signature(object = "modelorg"),
 
     }
 )
+#------------------------------------------------------------------------------#
 
+setMethod("getReaction", signature(X = "modelorg"),
+	function(X, j = NULL, drop=T, tol = SYBIL_SETTINGS("TOLERANCE")) {
+		# translate reaction id's to indices
+		cj <- checkReactId(X, react = j)
+		if (!is(cj, "reactId")) {
+			stop("check argument j")
+		}
+		else {
+			cn <- react_pos(cj)
+		}
+		rl <- lapply(cn, function(r){
+			s <- S(X)[,r]
+			new("react",
+				id=react_id(X)[r],
+				name=react_name(X)[r],
+				rev=react_rev(X)[r],
+				met_id=met_id(X)[abs(s) > tol],
+				met_name=met_name(X)[abs(s) > tol],
+				met_comp=met_comp(X)[abs(s) > tol],
+				s=s[abs(s) > tol],
+				lowbnd=lowbnd(X)[r],
+				uppbnd=uppbnd(X)[r],
+				obj_coef=obj_coef(X)[r],
+				gprRule=gprRules(X)[r],
+				genes=genes(X)[[r]],
+				gpr = gpr(X)[r],
+				subSys = colnames(subSys(X))[subSys(X)[r,]]
+			)
+		})
+		
+		if(length(rl) == 1 && drop){
+			return(rl[[1]])
+		}
+		return(rl)
+	}
+)
 
 #------------------------------------------------------------------------------#
 
diff --git a/R/reactClass.R b/R/reactClass.R
new file mode 100644
index 0000000000000000000000000000000000000000..13c7337870810314d3d7e06f570128f4d0819f03
--- /dev/null
+++ b/R/reactClass.R
@@ -0,0 +1,421 @@
+#  reactClass.R
+#  FBA and friends with R.
+#
+#  Copyright (C) 2010-2017 Claus Jonathan Fritzemeier, Dpt. for Computational Cell Biology,
+#  Institute for Informatics, Heinrich-Heine-University, Duesseldorf, Germany.
+#  All right reserved.
+#  Email: clausjonathan.fritzemeier@uni-duesseldorf.de
+#
+#  This file is part of sybil.
+#
+#  Sybil is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  Sybil is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with sybil.  If not, see <http://www.gnu.org/licenses/>.
+
+
+# reactClass
+
+
+#------------------------------------------------------------------------------#
+#					   definition of the class react						#
+#------------------------------------------------------------------------------#
+
+setClass("react",
+	representation(
+		 react_rev	  = "logical",	   # vector reversibilities
+		 react_id	  = "character",   # reaction id
+		 react_name	  = "character",   # reaction name
+		 react_single = "logical",	   # reaction using metabolites appearing only once in S
+		 react_de	  = "logical",	   # reaction using dead end metabolites
+		 met_id		  = "character",   # metabolites used in this reaction
+		 met_comp	  = "integer",   # compartments of metabolites
+		 met_name	  = "character",   # metabolite names
+		 s			  = "numeric",	   # matrix S
+		 lowbnd		  = "numeric",	   # reaction lower bound
+		 uppbnd		  = "numeric",	   # reaction upper bound
+		 obj_coef	  = "numeric",	   # objective coefficient
+		 gprRule	  = "character",
+		 genes		  = "character",
+		 gpr		  = "character",
+		 subSys		  = "character"
+
+	),
+	validity = .validreact
+)
+
+
+#------------------------------------------------------------------------------#
+#							 default constructor							   #
+#------------------------------------------------------------------------------#
+
+setMethod(f = "initialize",
+			signature = "react",
+			definition = function(.Object, 
+								id,
+								name="",
+								rev=TRUE,
+								single=NA,
+								de=NA,
+								met_id,
+								met_name=NULL,
+								met_comp=NULL,
+								s,
+								lowbnd=-1000,
+								uppbnd=1000,
+								obj_coef=0,
+								gprRule=NULL,
+								genes=NULL,
+								gpr = NULL,
+								subSys = NULL
+								) {
+			stopifnot(!missing(id))
+			stopifnot(!missing(met_id))
+			stopifnot(!missing(s))
+			
+			.Object@react_id <- id
+			.Object@react_name <- name
+			.Object@react_rev <- rev
+			.Object@react_single <- single
+			.Object@react_de <- de
+			.Object@met_id <- met_id
+			.Object@met_comp <- met_comp
+			.Object@met_name <- met_name
+			.Object@s <- s
+			.Object@lowbnd <- lowbnd
+			.Object@uppbnd <- uppbnd
+			.Object@obj_coef <- obj_coef
+			.Object@gprRule <- gprRule
+			.Object@genes <- genes
+			.Object@gpr <- gpr
+			.Object@subSys <- subSys
+			return(.Object)
+		  }
+)
+
+
+#------------------------------------------------------------------------------#
+#							 setters and getters							   #
+#------------------------------------------------------------------------------#
+
+# metabolite id's
+setMethod("met_id", signature(object = "react"),
+		  function(object) {
+			  return(object@met_id)
+		  }
+)
+
+setReplaceMethod("met_id", signature(object = "react"),
+		  function(object, value) {
+			  object@met_id <- value
+			  return(object)
+		  }
+)
+
+
+# metabolite names
+setMethod("met_name", signature(object = "react"),
+		  function(object) {
+			  return(object@met_name)
+		  }
+)
+
+setReplaceMethod("met_name", signature(object = "react"),
+		  function(object, value) {
+			  object@met_name <- value
+			  return(object)
+		  }
+)
+
+
+# metabolites compartments
+setMethod("met_comp", signature(object = "react"),
+		  function(object) {
+			  return(object@met_comp)
+		  }
+)
+
+setReplaceMethod("met_comp", signature(object = "react"),
+		  function(object, value) {
+			  object@met_comp <- value
+			  return(object)
+		  }
+)
+
+
+# reversibilities
+setMethod("react_rev", signature(object = "react"),
+		  function(object) {
+			  return(object@react_rev)
+		  }
+)
+
+setReplaceMethod("react_rev", signature(object = "react"),
+		  function(object, value) {
+			  object@react_rev <- value
+			  return(object)
+		  }
+)
+
+
+# reaction id's
+setMethod("react_id", signature(object = "react"),
+		  function(object) {
+			  return(object@react_id)
+		  }
+)
+
+setReplaceMethod("react_id", signature(object = "react"),
+		  function(object, value) {
+			  object@react_id <- value
+			  return(object)
+		  }
+)
+
+
+# reaction names
+setMethod("react_name", signature(object = "react"),
+		  function(object) {
+			  return(object@react_name)
+		  }
+)
+
+setReplaceMethod("react_name", signature(object = "react"),
+		  function(object, value) {
+			  object@react_name <- value
+			  return(object)
+		  }
+)
+
+
+# singletons
+setMethod("react_single", signature(object = "react"),
+		  function(object) {
+			  return(object@react_single)
+		  }
+)
+
+setReplaceMethod("react_single", signature(object = "react"),
+		  function(object, value) {
+			  object@react_single <- value
+			  return(object)
+		  }
+)
+
+
+# dead ends
+setMethod("react_de", signature(object = "react"),
+		  function(object) {
+			  return(object@react_de)
+		  }
+)
+
+setReplaceMethod("react_de", signature(object = "react"),
+		  function(object, value) {
+			  object@react_de <- value
+			  return(object)
+		  }
+)
+
+
+# stoichiometric matrix
+setMethod("s", signature(object = "react"),
+		  function(object) {
+			  return(object@s)
+		  }
+)
+
+setReplaceMethod("s", signature(object = "react"),
+		  function(object, value) {
+			  object@s <- value
+			  return(object)
+		  }
+)
+
+# lower bounds
+setMethod("lowbnd", signature(object = "react"),
+		  function(object) {
+			  return(object@lowbnd)
+		  }
+)
+
+setReplaceMethod("lowbnd", signature(object = "react"),
+		  function(object, value) {
+			  object@lowbnd <- value
+			  return(object)
+		  }
+)
+
+# upper bounds
+setMethod("uppbnd", signature(object = "react"),
+		  function(object) {
+			  return(object@uppbnd)
+		  }
+)
+
+setReplaceMethod("uppbnd", signature(object = "react"),
+		  function(object, value) {
+			  object@uppbnd <- value
+			  return(object)
+		  }
+)
+
+
+# objective coefficient
+setMethod("obj_coef", signature(object = "react"),
+		  function(object) {
+			  return(object@obj_coef)
+		  }
+)
+
+setReplaceMethod("obj_coef", signature(object = "react"),
+		  function(object, value) {
+			  object@obj_coef <- value
+			  return(object)
+		  }
+)
+
+
+# gprRules
+setMethod("gprRule", signature(object = "react"),
+		  function(object) {
+			  return(object@gprRule)
+		  }
+)
+
+setReplaceMethod("gprRule", signature(object = "react"),
+		  function(object, value) {
+			  object@gprRule <- value
+			  return(object)
+		  }
+)
+
+
+# genes
+setMethod("genes", signature(object = "react"),
+		  function(object) {
+			  return(object@genes)
+		  }
+)
+
+setReplaceMethod("genes", signature(object = "react"),
+		  function(object, value) {
+			  object@genes <- value
+			  return(object)
+		  }
+)
+
+
+# gpr associations
+setMethod("gpr", signature(object = "react"),
+		  function(object) {
+			  return(object@gpr)
+		  }
+)
+
+setReplaceMethod("gpr", signature(object = "react"),
+		  function(object, value) {
+			  object@gpr <- value
+			  return(object)
+		  }
+)
+
+
+# reaction sub systems
+setMethod("subSys", signature(object = "react"),
+		  function(object) {
+			  return(object@subSys)
+		  }
+)
+
+setReplaceMethod("subSys", signature(object = "react"),
+		  function(object, value) {
+			  object@subSys <- value
+			  return(object)
+		  }
+)
+
+
+#------------------------------------------------------------------------------#
+#								other methods								   #
+#------------------------------------------------------------------------------#
+
+setMethod("show", signature(object = "react"),
+	function(object) {
+		cat("react id:              ", react_id(object), "\n")
+		cat("react name:            ", react_name(object), "\n")
+		cat("stoichiometry:\n")
+		print(data.frame(met_id=met_id(object), s_coef=s(object)))
+		cat("lower bound:           ", lowbnd(object), "\n")
+		cat("upper bound:           ", uppbnd(object), "\n")
+		cat("objective function:    ", obj_coef(object), "\n")
+	}
+)
+
+
+#------------------------------------------------------------------------------#
+
+# print reactions
+setMethod("printReaction", signature(object = "react"),
+	function(object, printOut = TRUE, ...) {
+	
+		mat <- s(object)
+		reaction <- character(1)
+
+		for (j in seq(along = cind)) {
+		
+			met <- met_id(object)
+			nzv <- mat
+			
+			ed <- nzv < 0
+			pd <- nzv > 0
+
+			if (sum(ed) > 0) {
+				educt	<- paste(paste("(", abs(nzv[ed]), ")", sep = ""),
+								 met[ed], collapse = " + ")
+			}
+			else {
+				educt = ""
+			}
+
+			if (sum(pd) > 0) {
+				product <- paste(paste("(", nzv[pd], ")", sep = ""),
+								 met[pd], collapse = " + ")
+			}
+			else {
+				product = ""
+			}
+			
+			arrow	<- ifelse(react_rev(object)[cind[j]], " <==> ", " --> ")
+			
+			reaction[j] <- paste(react_id(check)[j],
+								 paste(educt, product, sep = arrow), sep = "\t")
+		}
+
+		if (isTRUE(printOut)) {
+		   cat("abbreviation\tequation", reaction, sep = "\n", ...)
+		}
+		
+		return(invisible(reaction))
+
+	}
+)
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/validreact.R b/R/validreact.R
new file mode 100644
index 0000000000000000000000000000000000000000..0036dd2923c5e198c89c3f2053862ee09d1a6752
--- /dev/null
+++ b/R/validreact.R
@@ -0,0 +1,112 @@
+#  validreact.R
+#  FBA and friends with R.
+#
+#  Copyright (C) 2010-2017 Claus Jonathan Fritzemeier, Dpt. for Computational Cell Biology,
+#  Institute for Informatics, Heinrich-Heine-University, Duesseldorf, Germany.
+#  All right reserved.
+#  Email: clausjonathan.fritzemeier@uni-duesseldorf.de
+#
+#  This file is part of sybil.
+#
+#  Sybil is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  Sybil is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with sybil.  If not, see <http://www.gnu.org/licenses/>.
+
+
+# validreact
+
+
+################################################
+# Function: .validreact
+#
+# Validity checking of an object of class react.
+#
+# Returns TRUE if the model is valid, otherwise
+# a character String containing a description of
+# the error.
+
+
+.validreact <- function(object) {
+	# data has to have the same length or to be NULL
+	if(!is(object, "react")){
+		"object is not of react class"
+	}
+	
+	if(length(object@id) != 1){
+		return("id has to be length 1")
+	}
+	if(length(object@name) != 1){
+		return("id has to be length 1")
+	}
+	
+	met_count <- length(object@s)
+	if(met_comp < 1){
+		return("reactions have to have at least one metabolite")
+	}
+	if(length(object@met_id) == met_count){
+		return("s, met_id, met_name, and met_comp have to have the same length")
+	}
+	if(length(object@met_name) == met_count && !is.null(object@met_name)){
+		return("s, met_id, met_name, and met_comp have to have the same length")
+	}
+	if(length(object@met_comp) == met_count && !is.null(object@met_comp)){
+		return("s, met_id, met_name, and met_comp have to have the same length")
+	}
+	
+	if(length(object@id)!=1){
+		return("length of id has to be 1")
+	}
+	if(length(object@rev)!=1){
+		return("length of rev has to be 1")
+	}
+	if(length(object@name)!=1){
+		return("length of name has to be 1")
+	}
+	if(length(object@lowbnd)!=1){
+		return("length of lowbnd has to be 1")
+	}
+	if(length(object@uppbnd)!=1){
+		return("length of uppbnd has to be 1")
+	}
+	if(length(object@obj_coef)!=1){
+		return("length of obj_coef has to be 1")
+	}
+	
+	if(!is.null(gprRule) && length(object@gprRule)!=1){
+		return("if not NULL, the length of gprRule has to be 1")
+	}
+	if(!is.null(genes) && length(object@genes)!=1){
+		return("if not NULL, the length of genes has to be 1")
+	}
+	if(!is.null(gpr) && length(object@gpr)!=1){
+		return("if not NULL, the length of gpr has to be 1")
+	}
+	if(!is.null(subSys) && length(object@subSys)!=1){
+		return("if not NULL, the length of subSys has to be 1")
+	}
+	
+	return(TRUE)
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+