diff --git a/R/readSBMLmod.R b/R/readSBMLmod.R
index 56e9a4125cef71a77b5a0b31e686bdd60d1e1ea8..c407a8897b8fcff8dce0eaeb58751ca407a0f5ea 100644
--- a/R/readSBMLmod.R
+++ b/R/readSBMLmod.R
@@ -258,9 +258,9 @@ parseNotesReact <- function(notes) {
           subSyst <- sub("SUBSYSTEM: *", "", fields_str[j])
           subSyst <- sub("^S_", "", subSyst, perl = TRUE)
           subSyst <- gsub("[_]+", " ", subSyst)
-          if (nchar(subSyst) == 0) {
-              subSyst <- "Exchange"
-          }
+          #if (nchar(subSyst) == 0) {
+          #    subSyst <- "Exchange"
+          #}
           #print(subSyst)
       }
   }
@@ -646,6 +646,9 @@ for (i in 1 : numreact) {
     
 }
 
+# get subsystem properties from the sbml groups plugin
+subSysGroups <- getSBMLGroupsList(sbmlmod)
+
 # ---------------------------------------------------------------------------- #
 # search for unused metabolites and unused reactions
 
@@ -964,7 +967,9 @@ else {
         sybil::gprRules(mod)   <- rules
         sybil::rxnGeneMat(mod) <- rxnGeneMat
         #sybil::subSys(mod)     <- subSys
-        sybil::subSys(mod)     <- sybil:::.prepareSubSysMatrix(subSys, numreact)
+        if(is.null(subSysGroups)){
+        	sybil::subSys(mod)     <- sybil:::.prepareSubSysMatrix(subSys, numreact)
+        }
 
         #sbml@gprRules <- rules
         #sbml@genes <- genes
@@ -1040,6 +1045,24 @@ if(newSybil)
   if( !is.null(reactnotes) && length(reactnotes)==numreact  )sybil::react_attr(mod)[['notes']]<-reactnotes
 }
 
+#------------------------------------------------------------------------------#
+#                         subSys from groups                                   #
+#------------------------------------------------------------------------------#
+
+if(!is.null(subSysGroups)){
+	subSysMat <- Matrix::Matrix(FALSE, nrow = numreact, ncol = length(subSysGroups), sparse = TRUE)
+	colnames(subSysMat) <- names(subSysGroups)
+	
+	subSysMat <- Matrix::Matrix( do.call("cbind",
+		lapply(subSysGroups, function(x){
+			ids <- formatSBMLid(sub( "^R[_]+", "", x))
+			sybil::react_id(mod) %in% ids
+		})),
+		sparse = TRUE
+	)
+	sybil::subSys(mod) <- subSysMat
+}
+
 
 #------------------------------------------------------------------------------#
 #                         Model Attr @Ardalan                                  #
diff --git a/R/sybilSBML.R b/R/sybilSBML.R
index b3eb40bee432ac7871489f4779c8fd2d086de6b1..8c1b02297e5f5c1753473670a87e1e4eff5c7141 100644
--- a/R/sybilSBML.R
+++ b/R/sybilSBML.R
@@ -222,6 +222,17 @@ getSBMLmodAnnotation <- function(sbmlm) {
 
 #------------------------------------------------------------------------------#
 
+getSBMLGroupsList <- function(sbmlm) {
+  
+  modGroups <- .Call("getSBMLGroupsList", PACKAGE = "sybilSBML",
+                sbmlPointer(sbmlm)
+  )
+  
+  return(modGroups)
+}
+
+
+#------------------------------------------------------------------------------#
 getSBMLnumCompart <- function(sbmlm) {
 
     num <- .Call("getSBMLnumCompart", PACKAGE = "sybilSBML",
@@ -519,15 +530,21 @@ writeSBML<- function(morg=NULL,level=2,version=4,fbcLevel=0,filename="export.xml
     
   } ####END newSybil attr
   
-  # Subsystem
-  if(is.null(newsubS) && !(modhasubSys) ) for ( i in 1:react_num(morg)) {newsubS[i]<- paste(names(which(subSys(morg)[i,])), collapse=", ")}
-  
+  # format ids for sbml
   newmet_id <- paste0("M_", (deformatSBMLid(met_id(morg))))
-  #newmet_id <- sub("\\[(\\w)\\]$", "_\\1",  newmet_id) # append compartment id, if in postfix with square brkts
-  
   newreact_id <- paste0("R_", deformatSBMLid(react_id(morg)))
   newmet_comp<-mod_compart(morg)[met_comp(morg)]
   
+  # Subsystem
+  if(is.null(newsubS) && !(modhasubSys) ) for ( i in 1:react_num(morg)) {newsubS[i]<- paste(names(which(subSys(morg)[i,])), collapse=", ")}
+  subSysGroups <- NULL
+  if(fbcLevel >= 2 && modhasubSys){
+  	subSysGroups <- apply(subSys(morg), 2, function(x){
+  		newreact_id[x]
+  	})
+  	names(subSysGroups) <- colnames(subSys(morg))
+  }
+  
   
   success <-.Call("exportSBML", PACKAGE = "sybilSBML",
                   as.integer(version),
@@ -550,6 +567,7 @@ writeSBML<- function(morg=NULL,level=2,version=4,fbcLevel=0,filename="export.xml
                   as.numeric(uppbnd(morg)),
                   as.integer(obj_coef(morg)),
                   as.character(newsubS),
+                  subSysGroups,
                   as.character(deformatGene(gpr(morg))),
                   as.numeric(shrinkMatrix(morg,j=1:react_num(morg))),
                   mod_notes,
@@ -572,4 +590,4 @@ writeSBML<- function(morg=NULL,level=2,version=4,fbcLevel=0,filename="export.xml
   }  
   else  message(paste("Could not write file ",filename,"\n",sep=""), appendLF = FALSE);
   return (success)
-}
\ No newline at end of file
+}
diff --git a/src/init.c b/src/init.c
index 7e4b174aafb9654417111f7c416117caf8ca2610..ef92c062b6264fd79e9641d5d914498c0d0d21eb 100644
--- a/src/init.c
+++ b/src/init.c
@@ -42,6 +42,7 @@ static const R_CallMethodDef callMethods[] = {
     {"getSBMLversion",             (DL_FUNC) &getSBMLversion,             1},
     {"validateDocument",           (DL_FUNC) &validateDocument,           1},
     {"getSBMLerrors",              (DL_FUNC) &getSBMLerrors,              1},
+    {"getSBMLGroupsList",          (DL_FUNC) &getSBMLGroupsList,          1},
     {"getSBMLmodel",               (DL_FUNC) &getSBMLmodel,               2},
     {"getSBMLmodId",               (DL_FUNC) &getSBMLmodId,               1},
     {"getSBMLmodName",             (DL_FUNC) &getSBMLmodName,             1},
@@ -54,7 +55,7 @@ static const R_CallMethodDef callMethods[] = {
     {"getSBMLCompartList",         (DL_FUNC) &getSBMLCompartList,         1},
     {"getSBMLSpeciesList",         (DL_FUNC) &getSBMLSpeciesList,         1},
     {"getSBMLReactionsList",       (DL_FUNC) &getSBMLReactionsList,       1},
-    {"exportSBML",                 (DL_FUNC) &exportSBML,                 33},
+    {"exportSBML",                 (DL_FUNC) &exportSBML,                 34},
     {"getSBMLFbcversion",          (DL_FUNC) &getSBMLFbcversion,          1},
     {NULL, NULL, 0}
 };
@@ -65,4 +66,4 @@ static const R_CallMethodDef callMethods[] = {
 void R_init_sybilSBML(DllInfo *info) {
     R_registerRoutines(info, NULL, callMethods, NULL, NULL);
     R_useDynamicSymbols(info, FALSE);
-}
\ No newline at end of file
+}
diff --git a/src/sybilSBML.c b/src/sybilSBML.c
index 0eee3b09a9bb8259d8c3b962876ff95f36bf043c..61b3d174e65d8735baa7b8eb0ca24be09b87b3a2 100644
--- a/src/sybilSBML.c
+++ b/src/sybilSBML.c
@@ -58,6 +58,17 @@ along with sybilSBML.  If not, see <http://www.gnu.org/licenses/>.
 #include <sbml/packages/fbc/sbml/FbcAnd.h>
 #include <sbml/packages/fbc/sbml/FbcOr.h>
 
+/*groups plugin*/
+#include <sbml/packages/groups/common/GroupsExtensionTypes.h>
+#include <sbml/packages/groups/extension/GroupsSBMLDocumentPlugin.h>
+#include <sbml/packages/groups/extension/GroupsExtension.h>
+#include <sbml/packages/groups/extension/GroupsModelPlugin.h>
+
+#include <sbml/packages/groups/sbml/Group.h>
+#include <sbml/packages/groups/sbml/ListOfGroups.h>
+#include <sbml/packages/groups/sbml/ListOfMembers.h>
+#include <sbml/packages/groups/sbml/Member.h>
+
 
 static SEXP tagSBMLmodel;
 static SEXP tagSBMLdocument;
@@ -1161,6 +1172,41 @@ SEXP getSBMLSpeciesList(SEXP sbmlmod) {
 }
 
 
+SEXP getSBMLGroupsList(SEXP sbmlmod) {
+	GroupsModelPlugin_t  * modelPlug = NULL;
+	modelPlug = (GroupsModelPlugin_t  *) SBase_getPlugin((SBase_t *)(R_ExternalPtrAddr(sbmlmod)), "groups");
+	
+	if(modelPlug != NULL){
+		int n = GroupsModelPlugin_getNumGroups(modelPlug);
+		SEXP rgroups = PROTECT(Rf_allocVector(VECSXP, n));
+		SEXP groupnames = PROTECT(Rf_allocVector(STRSXP, n));
+		for(int i=0; i<n; i++){
+			Group_t* group = GroupsModelPlugin_getGroup(modelPlug, i);
+			if(Group_isSetName(group) == 1){ // skip group if no name is set
+				SET_STRING_ELT(groupnames, i, Rf_mkChar(Group_getName(group)));
+				
+				int m = Group_getNumMembers(group);
+				SEXP rmembers = PROTECT(Rf_allocVector(STRSXP, m));
+				for(int j=0; j < m; j++){
+					Member_t * member = Group_getMember(group, j);
+					char * memref = Member_getIdRef(member);
+					SET_STRING_ELT(rmembers, j, Rf_mkChar(memref));
+				}
+				SET_VECTOR_ELT(rgroups, i, rmembers);
+				UNPROTECT(1);
+			}else{
+				SET_VECTOR_ELT(rgroups, i, R_NilValue);
+			}
+		}
+		Rf_namesgets(rgroups, groupnames);
+		UNPROTECT(2);
+		return rgroups;
+	}else{
+		return R_NilValue;
+	}
+	return R_NilValue;
+}
+
 /* -------------------------------------------------------------------------- */
 /* get list of reactions */
 SEXP getSBMLReactionsList(SEXP sbmlmod) {
@@ -1605,7 +1651,7 @@ void ParseModtoAnno  (SBase_t* comp , char* Mannocopy)
 
 
 
-SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybil_max, SEXP mod_desc, SEXP mod_name, SEXP mod_compart, SEXP met_id, SEXP met_name, SEXP met_comp, SEXP met_form,SEXP met_charge, SEXP react_id, SEXP react_name, SEXP react_rev, SEXP lowbnd, SEXP uppbnd, SEXP obj_coef, SEXP subSys, SEXP gpr, SEXP SMatrix, SEXP mod_notes, SEXP mod_anno, SEXP com_notes , SEXP com_anno, SEXP met_notes, SEXP met_anno, SEXP met_bnd , SEXP react_notes, SEXP react_anno, SEXP ex_react,SEXP allgenes)
+SEXP exportSBML (SEXP version, SEXP level, SEXP FbcLevel, SEXP filename, SEXP sybil_max, SEXP mod_desc, SEXP mod_name, SEXP mod_compart, SEXP met_id, SEXP met_name, SEXP met_comp, SEXP met_form, SEXP met_charge, SEXP react_id, SEXP react_name, SEXP react_rev, SEXP lowbnd, SEXP uppbnd, SEXP obj_coef, SEXP subSys, SEXP subSysGroups, SEXP gpr, SEXP SMatrix, SEXP mod_notes, SEXP mod_anno, SEXP com_notes , SEXP com_anno, SEXP met_notes, SEXP met_anno, SEXP met_bnd , SEXP react_notes, SEXP react_anno, SEXP ex_react, SEXP allgenes)
 {
   //Varaibles from R
   const char* fname = CHAR(STRING_ELT(filename, 0));
@@ -1615,6 +1661,7 @@ SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybi
   int SBMLlevel = INTEGER(level)[0];
   int SBMLversion = INTEGER(version)[0];
   int SBMLfbcversion = INTEGER(FbcLevel)[0];
+  int SBMLgroupsversion = 1;
   double sybilmax = REAL(sybil_max)[0];
   double sybilmin = sybilmax*(-1);
   
@@ -1680,12 +1727,17 @@ SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybi
       SBMLExtension_t *sbmlext = SBMLExtensionRegistry_getExtension("fbc");
       
       /* create the sbml namespaces object with fbc */
-      fbc = XMLNamespaces_create();
+      XMLNamespaces_t * fbc = XMLNamespaces_create();
       XMLNamespaces_add(fbc, SBMLExtension_getURI(sbmlext, 3, 1,  SBMLfbcversion), "fbc");
       
       sbmlns = SBMLNamespaces_create(3, 1);
       SBMLNamespaces_addNamespaces(sbmlns, fbc);
       
+      /* add groups extention */
+      SBMLExtension_t *sbmlgext = SBMLExtensionRegistry_getExtension("groups");
+      XMLNamespaces_t * groups = XMLNamespaces_create();
+      XMLNamespaces_add(groups, SBMLExtension_getURI(sbmlgext, SBMLlevel, SBMLversion, SBMLgroupsversion), "groups");
+      SBMLNamespaces_addNamespaces(sbmlns, groups);
       
       /* create the document */
       sbmlDoc = SBMLDocument_createWithSBMLNamespaces(sbmlns);
@@ -2231,7 +2283,7 @@ SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybi
       const double *lower_bnd = REAL(lowbnd);
       const double *upper_bnd = REAL(uppbnd);
       
-      char buf[20];
+      char buf[21]; // changed from 20 to 21 to avoid buffer overflow
       // FBC1 FLUXBOUNDS
       sprintf(buf, "LOWER_BOUND%d", i);
       if (INTEGER(obj_coef)[i] != 1)
@@ -2275,6 +2327,29 @@ SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybi
     }  
   }
   
+  /* add subsystem as groups if fbc is >= 2 */
+  if(SBMLfbcversion >= 2){
+  	if(!Rf_isNull(subSysGroups)){
+		GroupsModelPlugin_t* groupsPlug = NULL;
+  		groupsPlug = (GroupsModelPlugin_t*) SBase_getPlugin((SBase_t *)(model), "groups");
+  		
+  		for(int i=0; i < Rf_length(subSysGroups); i++){
+  			Group_t* newGroup = GroupsModelPlugin_createGroup(groupsPlug);
+			
+			Group_setKindAsString(newGroup, "partonomy");
+			Group_setName(newGroup, CHAR(STRING_ELT(Rf_getAttrib(subSysGroups, R_NamesSymbol), i)));
+			SBase_setSBOTerm((SBase_t *) newGroup, 0000633);
+  			
+  			for(int j=0; j < Rf_length(VECTOR_ELT(subSysGroups, i)); j++){
+				Member_t* newMember = Member_create(SBMLlevel, SBMLversion, SBMLgroupsversion);
+				Member_setIdRef(newMember, CHAR(STRING_ELT(VECTOR_ELT(subSysGroups, i), j)));
+				Group_addMember(newGroup, newMember);
+  			}
+  			//GroupsModelPlugin_addGroup(groupsPlug, newGroup);
+  		}
+  	}
+  }
+  
   // write SBML file
   int result = writeSBML(sbmlDoc, fname);
   SEXP out = R_NilValue;
@@ -2287,4 +2362,4 @@ SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybi
 
 
 
-/* -------------------------------------------------------------------------- */
\ No newline at end of file
+/* -------------------------------------------------------------------------- */
diff --git a/src/sybilSBML.h b/src/sybilSBML.h
index 9fac9ea3f09297661487956e0deca29afa6a2976..6393177a8b161b24de6c767be6dd20a71008bedd 100644
--- a/src/sybilSBML.h
+++ b/src/sybilSBML.h
@@ -105,8 +105,11 @@ SEXP getSBMLCompartList(SEXP sbmlmod);
 /* get list of species (metabolites) */
 SEXP getSBMLSpeciesList(SEXP sbmlmod);
 
+/* get list of groups */
+SEXP getSBMLGroupsList(SEXP sbmlmod);
+
 /* get list of reactions */
 SEXP getSBMLReactionsList(SEXP sbmlmod);
 
 /* export Modelorg to SBML*/
-SEXP exportSBML (SEXP version, SEXP level,SEXP FbcLevel, SEXP filename,SEXP sybil_max, SEXP mod_desc, SEXP mod_name, SEXP mod_compart, SEXP met_id, SEXP met_name, SEXP met_comp, SEXP met_form,SEXP met_charge, SEXP react_id, SEXP react_name, SEXP react_rev, SEXP lowbnd, SEXP uppbnd, SEXP obj_coef, SEXP subSys, SEXP gpr, SEXP SMatrix, SEXP mod_notes, SEXP mod_anno, SEXP com_notes , SEXP com_anno, SEXP met_notes, SEXP met_anno, SEXP met_bnd , SEXP react_notes, SEXP react_anno, SEXP ex_react, SEXP allgenes);
\ No newline at end of file
+SEXP exportSBML (SEXP version, SEXP level, SEXP FbcLevel, SEXP filename, SEXP sybil_max, SEXP mod_desc, SEXP mod_name, SEXP mod_compart, SEXP met_id, SEXP met_name, SEXP met_comp, SEXP met_form, SEXP met_charge, SEXP react_id, SEXP react_name, SEXP react_rev, SEXP lowbnd, SEXP uppbnd, SEXP obj_coef, SEXP subSys, SEXP subSysGroups, SEXP gpr, SEXP SMatrix, SEXP mod_notes, SEXP mod_anno, SEXP com_notes , SEXP com_anno, SEXP met_notes, SEXP met_anno, SEXP met_bnd , SEXP react_notes, SEXP react_anno, SEXP ex_react, SEXP allgenes);