diff --git a/cobra/core/__init__.py b/cobra/core/__init__.py
index f23d109a3..f9147b5f6 100644
--- a/cobra/core/__init__.py
+++ b/cobra/core/__init__.py
@@ -8,5 +8,6 @@
 from cobra.core.model import Model
 from cobra.core.object import Object
 from cobra.core.reaction import Reaction
+from cobra.core.group import Group
 from cobra.core.solution import Solution, LegacySolution, get_solution
 from cobra.core.species import Species
diff --git a/cobra/core/gene.py b/cobra/core/gene.py
index b869b06c4..921ff0d6e 100644
--- a/cobra/core/gene.py
+++ b/cobra/core/gene.py
@@ -248,6 +248,11 @@ def remove_from_model(self, model=None,
         the_gene_re = re.compile('(^|(?<=( |\()))%s(?=( |\)|$))' %
                                  re.escape(self.id))
 
+        # remove reference to the gene in all groups
+        associated_groups = self._model.get_associated_groups(self)
+        for group in associated_groups:
+            group.remove(self)
+
         self._model.genes.remove(self)
         self._model = None
 
diff --git a/cobra/core/group.py b/cobra/core/group.py
new file mode 100644
index 000000000..079eada4a
--- /dev/null
+++ b/cobra/core/group.py
@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+
+"""Define the group class."""
+
+from __future__ import absolute_import
+
+from warnings import warn
+
+from six import string_types
+
+from cobra.core.object import Object
+
+
+KIND_TYPES = ("collection", "classification", "partonomy")
+
+
+class Group(Object):
+    """
+    Manage groups via this implementation of the SBML group specification.
+
+    `Group` is a class for holding information regarding a pathways, subsystems,
+    or other custom groupings of objects within a cobra.Model object.
+
+    Parameters
+    ----------
+    id : str
+        The identifier to associate with this group
+    name : str, optional
+        A human readable name for the group
+    members : iterable, optional
+        A list object containing references to cobra.Model-associated objects
+        that belong to the group.
+    kind : {"collection", "classification", "partonomy"}, optional
+        The kind of group, as specified for the Groups feature in the SBML level
+        3 package specification. Can be any of "classification", "partonomy", or
+        "collection". The default is "collection". Please consult the SBML level
+        3 package specification to ensure you are using the proper value for
+        kind. In short, members of a "classification" group should have an
+        "is-a" relationship to the group (e.g. member is-a polar compound, or
+        member is-a transporter). Members of a "partonomy" group should have a
+        "part-of" relationship (e.g. member is part-of glycolysis). Members of a
+        "collection" group do not have an implied relationship between the
+        members, so use this value for kind when in doubt (e.g. member is a
+        gap-filled reaction, or member is involved in a disease phenotype).
+
+    """
+
+    def __init__(self, id, name='', members=None, kind=None):
+        Object.__init__(self, id, name)
+
+        self._members = set() if members is None else set(members)
+        self._kind = None
+        self.kind = "collection" if kind is None else kind
+        # self.model is None or refers to the cobra.Model that
+        # contains self
+        self._model = None
+
+    # read-only
+    @property
+    def members(self):
+        return self._members
+
+    @property
+    def kind(self):
+        return self._kind
+
+    @kind.setter
+    def kind(self, kind):
+        if kind in KIND_TYPES:
+            self._kind = kind
+        else:
+            raise ValueError(
+                "Kind can only by one of: {}.".format(", ".join(KIND_TYPES)))
+
+    def add_members(self, new_members):
+        """
+        Add objects to the group.
+
+        Parameters
+        ----------
+        new_members : list
+            A list of cobrapy objects to add to the group.
+
+        """
+
+        if isinstance(new_members, string_types) or \
+                hasattr(new_members, "id"):
+            warn("need to pass in a list")
+            new_members = [new_members]
+
+        self._members.update(new_members)
+
+    def remove(self, to_remove):
+        """
+        Remove objects from the group.
+
+        Parameters
+        ----------
+        to_remove : list
+            A list of cobrapy objects to remove from the group
+        """
+
+        if isinstance(to_remove, string_types) or \
+                hasattr(to_remove, "id"):
+            warn("need to pass in a list")
+            to_remove = [to_remove]
+
+        self._members.difference_update(to_remove)
diff --git a/cobra/core/model.py b/cobra/core/model.py
index 27e6a9f46..64d9b652c 100644
--- a/cobra/core/model.py
+++ b/cobra/core/model.py
@@ -17,6 +17,9 @@
 from cobra.core.dictlist import DictList
 from cobra.core.object import Object
 from cobra.core.reaction import separate_forward_and_reverse_bounds, Reaction
+from cobra.core.metabolite import Metabolite
+from cobra.core.gene import Gene
+from cobra.core.group import Group
 from cobra.core.solution import get_solution
 from cobra.util.context import HistoryManager, resettable, get_context
 from cobra.util.solver import (
@@ -51,6 +54,9 @@ class Model(Object):
     genes : DictList
         A DictList where the key is the gene identifier and the value a
         Gene
+    groups : DictList
+        A DictList where the key is the group identifier and the value a
+        Group
     solution : Solution
         The last obtained solution from optimizing the model.
     """
@@ -90,6 +96,7 @@ def __init__(self, id_or_model=None, name=None):
             self.genes = DictList()
             self.reactions = DictList()  # A list of cobra.Reactions
             self.metabolites = DictList()  # A list of cobra.Metabolites
+            self.groups = DictList()  # A list of cobra.Groups
             # genes based on their ids {Gene.id: Gene}
             self._compartments = dict()
             self._contexts = []
@@ -278,7 +285,7 @@ def copy(self):
         """
         new = self.__class__()
         do_not_copy_by_ref = {"metabolites", "reactions", "genes", "notes",
-                              "annotation"}
+                              "annotation", "groups"}
         for attr in self.__dict__:
             if attr not in do_not_copy_by_ref:
                 new.__dict__[attr] = self.__dict__[attr]
@@ -324,6 +331,28 @@ def copy(self):
                 new_gene = new.genes.get_by_id(gene.id)
                 new_reaction._genes.add(new_gene)
                 new_gene._reaction.add(new_reaction)
+
+        new.groups = DictList()
+        do_not_copy_by_ref = {"_model", "_members"}
+        for group in self.groups:
+            new_group = group.__class__()
+            for attr, value in iteritems(group.__dict__):
+                if attr not in do_not_copy_by_ref:
+                    new_group.__dict__[attr] = copy(value)
+            new_group._model = new
+            new.groups.append(new_group)
+            # update awareness, as in the reaction copies
+            new_objects = []
+            for member in group.members:
+                if isinstance(member, Metabolite):
+                    new_object = new.metabolites.get_by_id(member.id)
+                elif isinstance(member, Reaction):
+                    new_object = new.reactions.get_by_id(member.id)
+                elif isinstance(member, Gene):
+                    new_objext = new.genes.get_by_id(member.id)
+                new_objects.append(new_object)
+            new_group.add_members(new_objects)
+
         try:
             new._solver = deepcopy(self.solver)
             # Cplex has an issue with deep copies
@@ -406,6 +435,11 @@ def remove_metabolites(self, metabolite_list, destructive=False):
         for x in metabolite_list:
             x._model = None
 
+            # remove reference to the metabolite in all groups
+            associated_groups = self.get_associated_groups(x)
+            for group in associated_groups:
+                group.remove(x)
+
             if not destructive:
                 for the_reaction in list(x._reaction):
                     the_coefficient = the_reaction._metabolites[x]
@@ -644,6 +678,100 @@ def remove_reactions(self, reactions, remove_orphans=False):
                             if context:
                                 context(partial(self.genes.add, gene))
 
+                # remove reference to the reaction in all groups
+                associated_groups = self.get_associated_groups(reaction)
+                for group in associated_groups:
+                    group.remove(reaction)
+
+    def add_groups(self, group_list):
+        """Add groups to the model.
+
+        Groups with identifiers identical to a group already in the model are
+        ignored.
+
+        If any group contains members that are not in the model, these members
+        are added to the model as well. Only metabolites, reactions, and genes
+        can have groups.
+
+        Parameters
+        ----------
+        group_list : list
+            A list of `cobra.Group` objects to add to the model.
+        """
+
+        def existing_filter(group):
+            if group.id in self.groups:
+                LOGGER.warning(
+                    "Ignoring group '%s' since it already exists.", group.id)
+                return False
+            return True
+
+        if isinstance(group_list, string_types) or \
+                hasattr(group_list, "id"):
+            warn("need to pass in a list")
+            group_list = [group_list]
+
+        pruned = DictList(filter(existing_filter, group_list))
+
+        for group in pruned:
+            group._model = self
+            for member in group.members:
+                # If the member is not associated with the model, add it
+                if isinstance(member, Metabolite):
+                    if member not in self.metabolites:
+                        self.add_metabolites([member])
+                if isinstance(member, Reaction):
+                    if member not in self.reactions:
+                        self.add_reactions([member])
+                # TODO(midnighter): `add_genes` method does not exist.
+                # if isinstance(member, Gene):
+                #     if member not in self.genes:
+                #         self.add_genes([member])
+
+            self.groups += [group]
+
+    def remove_groups(self, group_list):
+        """Remove groups from the model.
+
+        Members of each group are not removed
+        from the model (i.e. metabolites, reactions, and genes in the group
+        stay in the model after any groups containing them are removed).
+
+        Parameters
+        ----------
+        group_list : list
+            A list of `cobra.Group` objects to remove from the model.
+        """
+
+        if isinstance(group_list, string_types) or \
+                hasattr(group_list, "id"):
+            warn("need to pass in a list")
+            group_list = [group_list]
+
+        for group in group_list:
+            # make sure the group is in the model
+            if group.id not in self.groups:
+                LOGGER.warning("%r not in %r. Ignored.", group, self)
+            else:
+                self.groups.remove(group)
+                group._model = None
+
+    def get_associated_groups(self, element):
+        """Returns a list of groups that an element (reaction, metabolite, gene)
+        is associated with.
+
+        Parameters
+        ----------
+        element: `cobra.Reaction`, `cobra.Metabolite`, or `cobra.Gene`
+
+        Returns
+        -------
+        list of `cobra.Group`
+            All groups that the provided object is a member of
+        """
+        # check whether the element is associated with the model
+        return [g for g in self.groups if element in g.members]
+
     def add_cons_vars(self, what, **kwargs):
         """Add constraints and variables to the model's mathematical problem.
 
@@ -891,6 +1019,7 @@ def repair(self, rebuild_index=True, rebuild_relationships=True):
             self.reactions._generate_index()
             self.metabolites._generate_index()
             self.genes._generate_index()
+            self.groups._generate_index()
         if rebuild_relationships:
             for met in self.metabolites:
                 met._reaction.clear()
@@ -901,8 +1030,9 @@ def repair(self, rebuild_index=True, rebuild_relationships=True):
                     met._reaction.add(rxn)
                 for gene in rxn._genes:
                     gene._reaction.add(rxn)
+
         # point _model to self
-        for l in (self.reactions, self.genes, self.metabolites):
+        for l in (self.reactions, self.genes, self.metabolites, self.groups):
             for e in l:
                 e._model = self
 
@@ -1077,6 +1207,9 @@ def _repr_html_(self):
             </tr><tr>
                 <td><strong>Number of reactions</strong></td>
                 <td>{num_reactions}</td>
+            </tr><tr>
+                <td><strong>Number of groups</strong></td>
+                <td>{num_groups}</td>
             </tr><tr>
                 <td><strong>Objective expression</strong></td>
                 <td>{objective}</td>
@@ -1089,6 +1222,7 @@ def _repr_html_(self):
             address='0x0%x' % id(self),
             num_metabolites=len(self.metabolites),
             num_reactions=len(self.reactions),
+            num_groups=len(self.groups),
             objective=format_long_string(str(self.objective.expression), 100),
             compartments=", ".join(
                 v if v else k for k, v in iteritems(self.compartments)
diff --git a/cobra/manipulation/delete.py b/cobra/manipulation/delete.py
index ea940e21f..4a5deec06 100644
--- a/cobra/manipulation/delete.py
+++ b/cobra/manipulation/delete.py
@@ -231,4 +231,8 @@ def remove_genes(cobra_model, gene_list, remove_reactions=True):
                 reaction.gene_reaction_rule = new_rule
     for gene in gene_set:
         cobra_model.genes.remove(gene)
+        # remove reference to the gene in all groups
+        associated_groups = cobra_model.get_associated_groups(gene)
+        for group in associated_groups:
+            group.remove(gene)
     cobra_model.remove_reactions(target_reactions)
diff --git a/cobra/test/data/iJO1366.pickle b/cobra/test/data/iJO1366.pickle
index 2e89067fa..f90e88ee8 100644
Binary files a/cobra/test/data/iJO1366.pickle and b/cobra/test/data/iJO1366.pickle differ
diff --git a/cobra/test/data/mini.mat b/cobra/test/data/mini.mat
index 6726342d3..7653038e3 100644
Binary files a/cobra/test/data/mini.mat and b/cobra/test/data/mini.mat differ
diff --git a/cobra/test/data/mini.pickle b/cobra/test/data/mini.pickle
index af9ef9ccc..c90dd6c5e 100644
Binary files a/cobra/test/data/mini.pickle and b/cobra/test/data/mini.pickle differ
diff --git a/cobra/test/data/mini_cobra.xml b/cobra/test/data/mini_cobra.xml
index 48b02294b..4e7486177 100644
--- a/cobra/test/data/mini_cobra.xml
+++ b/cobra/test/data/mini_cobra.xml
@@ -205,8 +205,8 @@
           <speciesReference species="M_h2o_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_adp_c" stoichiometry="1"/>
           <speciesReference species="M_pi_c" stoichiometry="1"/>
+          <speciesReference species="M_adp_c" stoichiometry="1"/>
           <speciesReference species="M_h_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
@@ -216,20 +216,20 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="8.39" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="1" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_D_LACt2" name="D_LACt2" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3603 or b2975</p>
+            <p>GENE ASSOCIATION: b3603 or b2975</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_h_e" stoichiometry="1"/>
           <speciesReference species="M_lac__D_e" stoichiometry="1"/>
+          <speciesReference species="M_h_e" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_lac__D_c" stoichiometry="1"/>
@@ -242,23 +242,23 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_ENO" name="enolase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2779</p>
+            <p>GENE ASSOCIATION: b2779</p>
           </html>
         </notes>
         <listOfReactants>
           <speciesReference species="M_2pg_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_h2o_c" stoichiometry="1"/>
           <speciesReference species="M_pep_c" stoichiometry="1"/>
+          <speciesReference species="M_h2o_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -267,8 +267,8 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
@@ -286,8 +286,8 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-10" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
@@ -305,8 +305,8 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
@@ -324,23 +324,23 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="0" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_FBA" name="fructose-bisphosphate aldolase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1773 or b2097 or b2925</p>
+            <p>GENE ASSOCIATION: b1773 or b2097 or b2925</p>
           </html>
         </notes>
         <listOfReactants>
           <speciesReference species="M_fdp_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_g3p_c" stoichiometry="1"/>
           <speciesReference species="M_dhap_c" stoichiometry="1"/>
+          <speciesReference species="M_g3p_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -349,26 +349,26 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_GAPD" name="glyceraldehyde-3-phosphate dehydrogenase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1779</p>
+            <p>GENE ASSOCIATION: b1779</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_nad_c" stoichiometry="1"/>
           <speciesReference species="M_pi_c" stoichiometry="1"/>
+          <speciesReference species="M_nad_c" stoichiometry="1"/>
           <speciesReference species="M_g3p_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_13dpg_c" stoichiometry="1"/>
-          <speciesReference species="M_h_c" stoichiometry="1"/>
           <speciesReference species="M_nadh_c" stoichiometry="1"/>
+          <speciesReference species="M_h_c" stoichiometry="1"/>
+          <speciesReference species="M_13dpg_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -377,20 +377,20 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_GLCpts" name="D-glucose transport via PEP:Pyr PTS" reversible="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: ( b2417 and b1621 and b2415 and b2416 ) or ( b2417 and b1101 and b2415 and b2416 ) or ( b1817 and b1818 and b1819 and b2415 and b2416 )</p>
+            <p>GENE ASSOCIATION: ( b2417 and b1621 and b2415 and b2416 ) or ( b2417 and b1101 and b2415 and b2416 ) or ( b1817 and b1818 and b1819 and b2415 and b2416 )</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_pep_c" stoichiometry="1"/>
           <speciesReference species="M_glc__D_e" stoichiometry="1"/>
+          <speciesReference species="M_pep_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_g6p_c" stoichiometry="1"/>
@@ -403,15 +403,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="0" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_H2Ot" name="R H2O transport via - diffusion" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b0875 or s0001</p>
+            <p>GENE ASSOCIATION: b0875 or s0001</p>
           </html>
         </notes>
         <listOfReactants>
@@ -427,20 +427,20 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_LDH_D" name="D-lactate dehydrogenase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2133 or b1380</p>
+            <p>GENE ASSOCIATION: b2133 or b1380</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_nad_c" stoichiometry="1"/>
           <speciesReference species="M_lac__D_c" stoichiometry="1"/>
+          <speciesReference species="M_nad_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_nadh_c" stoichiometry="1"/>
@@ -454,25 +454,25 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PFK" name="phosphofructokinase" reversible="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3916 or b1723</p>
+            <p>GENE ASSOCIATION: b3916 or b1723</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_atp_c" stoichiometry="1"/>
           <speciesReference species="M_f6p_c" stoichiometry="1"/>
+          <speciesReference species="M_atp_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_adp_c" stoichiometry="1"/>
-          <speciesReference species="M_fdp_c" stoichiometry="1"/>
           <speciesReference species="M_h_c" stoichiometry="1"/>
+          <speciesReference species="M_fdp_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -481,15 +481,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="0" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="1" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PGI" name="glucose-6-phosphate isomerase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b4025</p>
+            <p>GENE ASSOCIATION: b4025</p>
           </html>
         </notes>
         <listOfReactants>
@@ -505,15 +505,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PGK" name="phosphoglycerate kinase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2926</p>
+            <p>GENE ASSOCIATION: b2926</p>
           </html>
         </notes>
         <listOfReactants>
@@ -521,8 +521,8 @@
           <speciesReference species="M_3pg_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_13dpg_c" stoichiometry="1"/>
           <speciesReference species="M_adp_c" stoichiometry="1"/>
+          <speciesReference species="M_13dpg_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -531,15 +531,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PGM" name="phosphoglycerate mutase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b4395 or b3612 or b0755</p>
+            <p>GENE ASSOCIATION: b4395 or b3612 or b0755</p>
           </html>
         </notes>
         <listOfReactants>
@@ -555,15 +555,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PIt2r" name="R phosphate reversible transport via - symport" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2987 or b3493</p>
+            <p>GENE ASSOCIATION: b2987 or b3493</p>
           </html>
         </notes>
         <listOfReactants>
@@ -581,25 +581,25 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_PYK" name="pyruvate kinase" reversible="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1854 or b1676</p>
+            <p>GENE ASSOCIATION: b1854 or b1676</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_adp_c" stoichiometry="1"/>
           <speciesReference species="M_pep_c" stoichiometry="1"/>
+          <speciesReference species="M_adp_c" stoichiometry="1"/>
           <speciesReference species="M_h_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_atp_c" stoichiometry="1"/>
           <speciesReference species="M_pyr_c" stoichiometry="1"/>
+          <speciesReference species="M_atp_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -608,15 +608,15 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="0" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
       <reaction id="R_TPI" name="triose-phosphate isomerase" reversible="true">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3919</p>
+            <p>GENE ASSOCIATION: b3919</p>
           </html>
         </notes>
         <listOfReactants>
@@ -632,8 +632,8 @@
           <listOfParameters>
             <parameter id="LOWER_BOUND" value="-1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="FLUX_VALUE" value="0" units="mmol_per_gDW_per_hr"/>
-            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
             <parameter id="OBJECTIVE_COEFFICIENT" value="0" units="dimensionless"/>
+            <parameter id="UPPER_BOUND" value="1000" units="mmol_per_gDW_per_hr"/>
           </listOfParameters>
         </kineticLaw>
       </reaction>
diff --git a/cobra/test/data/mini_fbc1.xml b/cobra/test/data/mini_fbc1.xml
index 58fb4bc9c..c286c4478 100644
--- a/cobra/test/data/mini_fbc1.xml
+++ b/cobra/test/data/mini_fbc1.xml
@@ -312,20 +312,20 @@
           <speciesReference species="M_h2o_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pi_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_D_LACt2" name="D_LACt2" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3603 or b2975</p>
+            <p>GENE ASSOCIATION: b3603 or b2975</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_h_e" stoichiometry="1" constant="true"/>
           <speciesReference species="M_lac__D_e" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_h_e" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_lac__D_c" stoichiometry="1" constant="true"/>
@@ -335,15 +335,15 @@
       <reaction id="R_ENO" name="enolase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2779</p>
+            <p>GENE ASSOCIATION: b2779</p>
           </html>
         </notes>
         <listOfReactants>
           <speciesReference species="M_2pg_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_h2o_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pep_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_h2o_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_EX_glc__D_e" name="D-Glucose exchange" reversible="true" fast="false">
@@ -373,43 +373,43 @@
       <reaction id="R_FBA" name="fructose-bisphosphate aldolase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1773 or b2097 or b2925</p>
+            <p>GENE ASSOCIATION: b1773 or b2097 or b2925</p>
           </html>
         </notes>
         <listOfReactants>
           <speciesReference species="M_fdp_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_g3p_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_dhap_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_g3p_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_GAPD" name="glyceraldehyde-3-phosphate dehydrogenase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1779</p>
+            <p>GENE ASSOCIATION: b1779</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_nad_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pi_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_nad_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_g3p_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_13dpg_c" stoichiometry="1" constant="true"/>
-          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_nadh_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_13dpg_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_GLCpts" name="D-glucose transport via PEP:Pyr PTS" reversible="false" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: ( b2417 and b1621 and b2415 and b2416 ) or ( b2417 and b1101 and b2415 and b2416 ) or ( b1817 and b1818 and b1819 and b2415 and b2416 )</p>
+            <p>GENE ASSOCIATION: ( b2417 and b1621 and b2415 and b2416 ) or ( b2417 and b1101 and b2415 and b2416 ) or ( b1817 and b1818 and b1819 and b2415 and b2416 )</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_pep_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_glc__D_e" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_pep_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_g6p_c" stoichiometry="1" constant="true"/>
@@ -419,7 +419,7 @@
       <reaction id="R_H2Ot" name="R H2O transport via - diffusion" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b0875 or s0001</p>
+            <p>GENE ASSOCIATION: b0875 or s0001</p>
           </html>
         </notes>
         <listOfReactants>
@@ -432,12 +432,12 @@
       <reaction id="R_LDH_D" name="D-lactate dehydrogenase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2133 or b1380</p>
+            <p>GENE ASSOCIATION: b2133 or b1380</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_nad_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_lac__D_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_nad_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_nadh_c" stoichiometry="1" constant="true"/>
@@ -448,23 +448,23 @@
       <reaction id="R_PFK" name="phosphofructokinase" reversible="false" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3916 or b1723</p>
+            <p>GENE ASSOCIATION: b3916 or b1723</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_f6p_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
-          <speciesReference species="M_fdp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_fdp_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_PGI" name="glucose-6-phosphate isomerase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b4025</p>
+            <p>GENE ASSOCIATION: b4025</p>
           </html>
         </notes>
         <listOfReactants>
@@ -477,7 +477,7 @@
       <reaction id="R_PGK" name="phosphoglycerate kinase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2926</p>
+            <p>GENE ASSOCIATION: b2926</p>
           </html>
         </notes>
         <listOfReactants>
@@ -485,14 +485,14 @@
           <speciesReference species="M_3pg_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_13dpg_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_13dpg_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_PGM" name="phosphoglycerate mutase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b4395 or b3612 or b0755</p>
+            <p>GENE ASSOCIATION: b4395 or b3612 or b0755</p>
           </html>
         </notes>
         <listOfReactants>
@@ -505,7 +505,7 @@
       <reaction id="R_PIt2r" name="R phosphate reversible transport via - symport" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b2987 or b3493</p>
+            <p>GENE ASSOCIATION: b2987 or b3493</p>
           </html>
         </notes>
         <listOfReactants>
@@ -520,23 +520,23 @@
       <reaction id="R_PYK" name="pyruvate kinase" reversible="false" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b1854 or b1676</p>
+            <p>GENE ASSOCIATION: b1854 or b1676</p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pep_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pyr_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <reaction id="R_TPI" name="triose-phosphate isomerase" reversible="true" fast="false">
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: b3919</p>
+            <p>GENE ASSOCIATION: b3919</p>
           </html>
         </notes>
         <listOfReactants>
diff --git a/cobra/test/data/mini_fbc2.xml.gz b/cobra/test/data/mini_fbc2.xml.gz
index 16d383c00..0fcd44c55 100644
Binary files a/cobra/test/data/mini_fbc2.xml.gz and b/cobra/test/data/mini_fbc2.xml.gz differ
diff --git a/cobra/test/data/raven.pickle b/cobra/test/data/raven.pickle
index adfaf9b65..c5e250b55 100644
Binary files a/cobra/test/data/raven.pickle and b/cobra/test/data/raven.pickle differ
diff --git a/cobra/test/data/salmonella.pickle b/cobra/test/data/salmonella.pickle
index ba0568c82..ae44cc013 100644
Binary files a/cobra/test/data/salmonella.pickle and b/cobra/test/data/salmonella.pickle differ
diff --git a/cobra/test/data/textbook_fva.json b/cobra/test/data/textbook_fva.json
index f31092a84..3ab5fe4c2 100644
--- a/cobra/test/data/textbook_fva.json
+++ b/cobra/test/data/textbook_fva.json
@@ -1 +1 @@
-{"maximum": {"EX_fum_e": 0.0, "ACALDt": 0.0, "EX_glc__D_e": -10.0, "EX_mal__L_e": -0.0, "ADK1": -0.0, "ICL": -0.0, "TALA": 1.49698, "EX_ac_e": -0.0, "PGI": 4.86086, "ACKr": 0.0, "NADTRHD": -0.0, "SUCCt2_2": -0.0, "O2t": 21.79949, "EX_co2_e": 22.80983, "PTAr": -0.0, "EX_h2o_e": 29.17583, "GLUDy": -4.54186, "ACONTa": 6.00725, "GLCpts": 10.0, "GAPD": 16.02353, "TKT1": 1.49698, "TKT2": 1.1815, "NADH16": 38.53461, "EX_etoh_e": -0.0, "ME1": -0.0, "FBP": -0.0, "GLUt2r": 0.0, "SUCDi": 1000.0, "EX_h_e": 17.53087, "ACt2r": 0.0, "GLUSy": -0.0, "TPI": 7.47738, "PYRt2": 0.0, "PGM": -14.71614, "Biomass_Ecoli_core": 0.87392, "PFL": -0.0, "RPE": 2.67848, "RPI": -2.2815, "EX_succ_e": -0.0, "ACONTb": 6.00725, "EX_lac__D_e": -0.0, "PPC": 2.50431, "ALCD2x": 0.0, "AKGDH": 5.06438, "EX_acald_e": -0.0, "EX_nh4_e": -4.76532, "GLUN": -0.0, "EX_gln__L_e": 0.0, "EX_glu__L_e": -0.0, "GND": 4.95998, "PGL": 4.95998, "PPCK": -0.0, "ENO": 14.71614, "EX_fru_e": -0.0, "AKGt2r": 0.0, "SUCCt3": -0.0, "PDH": 9.28253, "EX_pyr_e": -0.0, "EX_o2_e": -21.79949, "PPS": -0.0, "H2Ot": -29.17583, "GLNabc": 0.0, "MDH": 5.06438, "EX_akg_e": -0.0, "ME2": -0.0, "FORt2": -0.0, "EX_for_e": -0.0, "SUCOAS": -5.06438, "PIt2r": 3.2149, "CS": 6.00725, "MALS": -0.0, "FBA": 7.47738, "FRUpts2": 0.0, "PYK": 1.75818, "ATPM": 8.39, "LDH_D": 0.0, "CYTBD": 43.59899, "NH4t": 4.76532, "CO2t": -22.80983, "THD2": -0.0, "ATPS4r": 45.51401, "D_LACt2": 0.0, "FRD7": 994.93562, "GLNS": 0.22346, "G6PDH2r": 4.95998, "MALt2_2": 0.0, "FORti": -0.0, "PFK": 7.47738, "ETOHt2r": 0.0, "ICDHyr": 6.00725, "PGK": -16.02353, "ACALD": 0.0, "FUMt2_2": 0.0, "FUM": 5.06438, "EX_pi_e": -3.2149}, "minimum": {"EX_fum_e": 0.0, "ACALDt": 0.0, "EX_glc__D_e": -10.0, "EX_mal__L_e": 0.0, "ADK1": 0.0, "ICL": 0.0, "TALA": 1.49698, "EX_ac_e": 0.0, "PGI": 4.86086, "ACKr": 0.0, "NADTRHD": 0.0, "SUCCt2_2": 0.0, "O2t": 21.79949, "EX_co2_e": 22.80983, "PTAr": 0.0, "EX_h2o_e": 29.17583, "GLUDy": -4.54186, "ACONTa": 6.00725, "GLCpts": 10.0, "GAPD": 16.02353, "TKT1": 1.49698, "TKT2": 1.1815, "NADH16": 38.53461, "EX_etoh_e": 0.0, "ME1": 0.0, "FBP": 0.0, "GLUt2r": 0.0, "SUCDi": 5.06438, "EX_h_e": 17.53087, "ACt2r": 0.0, "GLUSy": 0.0, "TPI": 7.47738, "PYRt2": 0.0, "PGM": -14.71614, "Biomass_Ecoli_core": 0.87392, "PFL": 0.0, "RPE": 2.67848, "RPI": -2.2815, "EX_succ_e": 0.0, "ACONTb": 6.00725, "EX_lac__D_e": 0.0, "PPC": 2.50431, "ALCD2x": 0.0, "AKGDH": 5.06438, "EX_acald_e": 0.0, "EX_nh4_e": -4.76532, "GLUN": 0.0, "EX_gln__L_e": 0.0, "EX_glu__L_e": 0.0, "GND": 4.95998, "PGL": 4.95998, "PPCK": 0.0, "ENO": 14.71614, "EX_fru_e": 0.0, "AKGt2r": 0.0, "SUCCt3": 0.0, "PDH": 9.28253, "EX_pyr_e": 0.0, "EX_o2_e": -21.79949, "PPS": 0.0, "H2Ot": -29.17583, "GLNabc": 0.0, "MDH": 5.06438, "EX_akg_e": 0.0, "ME2": 0.0, "FORt2": 0.0, "EX_for_e": 0.0, "SUCOAS": -5.06438, "PIt2r": 3.2149, "CS": 6.00725, "MALS": 0.0, "FBA": 7.47738, "FRUpts2": 0.0, "PYK": 1.75818, "ATPM": 8.39, "LDH_D": 0.0, "CYTBD": 43.59899, "NH4t": 4.76532, "CO2t": -22.80983, "THD2": 0.0, "ATPS4r": 45.51401, "D_LACt2": 0.0, "FRD7": 0.0, "GLNS": 0.22346, "G6PDH2r": 4.95998, "MALt2_2": 0.0, "FORti": -0.0, "PFK": 7.47738, "ETOHt2r": 0.0, "ICDHyr": 6.00725, "PGK": -16.02353, "ACALD": 0.0, "FUMt2_2": 0.0, "FUM": 5.06438, "EX_pi_e": -3.2149}}
\ No newline at end of file
+{"maximum": {"PTAr": -0.0, "ACONTb": 6.00725, "ACALDt": 0.0, "CYTBD": 43.59899, "FORt2": -0.0, "SUCDi": 1000.0, "GLCpts": 10.0, "EX_etoh_e": -0.0, "ACONTa": 6.00725, "EX_nh4_e": -4.76532, "FBA": 7.47738, "CS": 6.00725, "FRD7": 994.93562, "MDH": 5.06438, "GLNS": 0.22346, "ADK1": -0.0, "PPCK": -0.0, "EX_glu__L_e": -0.0, "FRUpts2": 0.0, "EX_co2_e": 22.80983, "EX_acald_e": -0.0, "EX_pyr_e": -0.0, "FUMt2_2": 0.0, "THD2": -0.0, "O2t": 21.79949, "RPE": 2.67848, "GAPD": 16.02353, "GND": 4.95998, "ICDHyr": 6.00725, "NH4t": 4.76532, "PGK": -16.02353, "EX_fru_e": -0.0, "PGM": -14.71614, "EX_glc__D_e": -10.0, "CO2t": -22.80983, "MALt2_2": 0.0, "EX_mal__L_e": -0.0, "ATPS4r": 45.51401, "FBP": -0.0, "EX_succ_e": -0.0, "D_LACt2": 0.0, "PYRt2": 0.0, "NADTRHD": -0.0, "TALA": 1.49698, "ATPM": 8.39, "GLUN": -0.0, "AKGDH": 5.06438, "PPC": 2.50431, "PIt2r": 3.2149, "G6PDH2r": 4.95998, "ETOHt2r": 0.0, "EX_ac_e": -0.0, "GLUDy": -4.54186, "ACALD": 0.0, "EX_pi_e": -3.2149, "TKT1": 1.49698, "ACKr": 0.0, "PDH": 9.28253, "EX_for_e": -0.0, "ME2": -0.0, "PGL": 4.95998, "ENO": 14.71614, "TKT2": 1.1815, "EX_h_e": 17.53087, "GLUt2r": 0.0, "SUCOAS": -5.06438, "FUM": 5.06438, "PGI": 4.86086, "ALCD2x": 0.0, "EX_gln__L_e": 0.0, "FORti": -0.0, "ICL": -0.0, "ME1": -0.0, "GLNabc": 0.0, "PFK": 7.47738, "EX_akg_e": -0.0, "NADH16": 38.53461, "PYK": 1.75818, "AKGt2r": 0.0, "MALS": -0.0, "EX_o2_e": -21.79949, "SUCCt3": -0.0, "RPI": -2.2815, "EX_lac__D_e": -0.0, "PFL": -0.0, "ACt2r": 0.0, "H2Ot": -29.17583, "LDH_D": 0.0, "SUCCt2_2": -0.0, "TPI": 7.47738, "GLUSy": -0.0, "Biomass_Ecoli_core": 0.87392, "PPS": -0.0, "EX_h2o_e": 29.17583, "EX_fum_e": 0.0}, "minimum": {"PTAr": 0.0, "ACONTb": 6.00725, "ACALDt": 0.0, "CYTBD": 43.59899, "FORt2": 0.0, "SUCDi": 5.06438, "GLCpts": 10.0, "EX_etoh_e": 0.0, "ACONTa": 6.00725, "EX_nh4_e": -4.76532, "FBA": 7.47738, "CS": 6.00725, "FRD7": 0.0, "MDH": 5.06438, "GLNS": 0.22346, "ADK1": 0.0, "PPCK": 0.0, "EX_glu__L_e": 0.0, "FRUpts2": 0.0, "EX_co2_e": 22.80983, "EX_acald_e": 0.0, "EX_pyr_e": 0.0, "FUMt2_2": 0.0, "THD2": 0.0, "O2t": 21.79949, "RPE": 2.67848, "GAPD": 16.02353, "GND": 4.95998, "ICDHyr": 6.00725, "NH4t": 4.76532, "PGK": -16.02353, "EX_fru_e": 0.0, "PGM": -14.71614, "EX_glc__D_e": -10.0, "CO2t": -22.80983, "MALt2_2": 0.0, "EX_mal__L_e": 0.0, "ATPS4r": 45.51401, "FBP": 0.0, "EX_succ_e": 0.0, "D_LACt2": 0.0, "PYRt2": 0.0, "NADTRHD": 0.0, "TALA": 1.49698, "ATPM": 8.39, "GLUN": 0.0, "AKGDH": 5.06438, "PPC": 2.50431, "PIt2r": 3.2149, "G6PDH2r": 4.95998, "ETOHt2r": 0.0, "EX_ac_e": 0.0, "GLUDy": -4.54186, "ACALD": 0.0, "EX_pi_e": -3.2149, "TKT1": 1.49698, "ACKr": 0.0, "PDH": 9.28253, "EX_for_e": 0.0, "ME2": 0.0, "PGL": 4.95998, "ENO": 14.71614, "TKT2": 1.1815, "EX_h_e": 17.53087, "GLUt2r": 0.0, "SUCOAS": -5.06438, "FUM": 5.06438, "PGI": 4.86086, "ALCD2x": 0.0, "EX_gln__L_e": 0.0, "FORti": 0.0, "ICL": 0.0, "ME1": 0.0, "GLNabc": 0.0, "PFK": 7.47738, "EX_akg_e": 0.0, "NADH16": 38.53461, "PYK": 1.75818, "AKGt2r": 0.0, "MALS": 0.0, "EX_o2_e": -21.79949, "SUCCt3": 0.0, "RPI": -2.2815, "EX_lac__D_e": 0.0, "PFL": 0.0, "ACt2r": 0.0, "H2Ot": -29.17583, "LDH_D": 0.0, "SUCCt2_2": 0.0, "TPI": 7.47738, "GLUSy": 0.0, "Biomass_Ecoli_core": 0.87392, "PPS": 0.0, "EX_h2o_e": 29.17583, "EX_fum_e": 0.0}}
\ No newline at end of file
diff --git a/cobra/test/data/textbook_pfba_fva.json b/cobra/test/data/textbook_pfba_fva.json
index c8acb0281..d97e0177e 100644
--- a/cobra/test/data/textbook_pfba_fva.json
+++ b/cobra/test/data/textbook_pfba_fva.json
@@ -1 +1 @@
-{"maximum": {"ACALD": 0.0, "ACALDt": 0.0, "ACKr": 0.0, "ACONTa": 6.00725, "ACONTb": 6.00725, "ACt2r": 0.0, "ADK1": -0.0, "AKGDH": 5.06438, "AKGt2r": 0.0, "ALCD2x": -0.0, "ATPM": 8.39, "ATPS4r": 45.51401, "Biomass_Ecoli_core": 0.87392, "CO2t": -22.80983, "CS": 6.00725, "CYTBD": 43.59899, "D_LACt2": 0.0, "ENO": 14.71614, "ETOHt2r": -0.0, "EX_ac_e": -0.0, "EX_acald_e": -0.0, "EX_akg_e": -0.0, "EX_co2_e": 22.80983, "EX_etoh_e": -0.0, "EX_for_e": -0.0, "EX_fru_e": -0.0, "EX_fum_e": 0.0, "EX_glc__D_e": -10.0, "EX_gln__L_e": 0.0, "EX_glu__L_e": -0.0, "EX_h2o_e": 29.17583, "EX_h_e": 17.53087, "EX_lac__D_e": -0.0, "EX_mal__L_e": -0.0, "EX_nh4_e": -4.76532, "EX_o2_e": -21.79949, "EX_pi_e": -3.2149, "EX_pyr_e": -0.0, "EX_succ_e": -0.0, "FBA": 7.47738, "FBP": -0.0, "FORt2": -0.0, "FORti": -0.0, "FRD7": 25.9211, "FRUpts2": 0.0, "FUM": 5.06438, "FUMt2_2": 0.0, "G6PDH2r": 4.95998, "GAPD": 16.02353, "GLCpts": 10.0, "GLNS": 0.22346, "GLNabc": 0.0, "GLUDy": -4.54186, "GLUN": -0.0, "GLUSy": -0.0, "GLUt2r": 0.0, "GND": 4.95998, "H2Ot": -29.17583, "ICDHyr": 6.00725, "ICL": -0.0, "LDH_D": 0.0, "MALS": -0.0, "MALt2_2": 0.0, "MDH": 5.06438, "ME1": -0.0, "ME2": -0.0, "NADH16": 38.53461, "NADTRHD": -0.0, "NH4t": 4.76532, "O2t": 21.79949, "PDH": 9.28253, "PFK": 7.47738, "PFL": -0.0, "PGI": 4.86086, "PGK": -16.02353, "PGL": 4.95998, "PGM": -14.71614, "PIt2r": 3.2149, "PPC": 2.50431, "PPCK": -0.0, "PPS": -0.0, "PTAr": -0.0, "PYK": 1.75818, "PYRt2": 0.0, "RPE": 2.67848, "RPI": -2.2815, "SUCCt2_2": -0.0, "SUCCt3": -0.0, "SUCDi": 30.98548, "SUCOAS": -5.06438, "TALA": 1.49698, "THD2": -0.0, "TKT1": 1.49698, "TKT2": 1.1815, "TPI": 7.47738}, "minimum": {"ACALD": 0.0, "ACALDt": 0.0, "ACKr": 0.0, "ACONTa": 6.00725, "ACONTb": 6.00725, "ACt2r": 0.0, "ADK1": 0.0, "AKGDH": 5.06438, "AKGt2r": 0.0, "ALCD2x": 0.0, "ATPM": 8.39, "ATPS4r": 45.51401, "Biomass_Ecoli_core": 0.87392, "CO2t": -22.80983, "CS": 6.00725, "CYTBD": 43.59899, "D_LACt2": 0.0, "ENO": 14.71614, "ETOHt2r": 0.0, "EX_ac_e": 0.0, "EX_acald_e": 0.0, "EX_akg_e": 0.0, "EX_co2_e": 22.80983, "EX_etoh_e": 0.0, "EX_for_e": 0.0, "EX_fru_e": 0.0, "EX_fum_e": 0.0, "EX_glc__D_e": -10.0, "EX_gln__L_e": 0.0, "EX_glu__L_e": 0.0, "EX_h2o_e": 29.17583, "EX_h_e": 17.53087, "EX_lac__D_e": 0.0, "EX_mal__L_e": 0.0, "EX_nh4_e": -4.76532, "EX_o2_e": -21.79949, "EX_pi_e": -3.2149, "EX_pyr_e": 0.0, "EX_succ_e": 0.0, "FBA": 7.47738, "FBP": 0.0, "FORt2": 0.0, "FORti": 0.0, "FRD7": 0.0, "FRUpts2": 0.0, "FUM": 5.06438, "FUMt2_2": 0.0, "G6PDH2r": 4.95998, "GAPD": 16.02353, "GLCpts": 10.0, "GLNS": 0.22346, "GLNabc": 0.0, "GLUDy": -4.54186, "GLUN": 0.0, "GLUSy": 0.0, "GLUt2r": 0.0, "GND": 4.95998, "H2Ot": -29.17583, "ICDHyr": 6.00725, "ICL": 0.0, "LDH_D": 0.0, "MALS": 0.0, "MALt2_2": 0.0, "MDH": 5.06438, "ME1": 0.0, "ME2": 0.0, "NADH16": 38.53461, "NADTRHD": 0.0, "NH4t": 4.76532, "O2t": 21.79949, "PDH": 9.28253, "PFK": 7.47738, "PFL": 0.0, "PGI": 4.86086, "PGK": -16.02353, "PGL": 4.95998, "PGM": -14.71614, "PIt2r": 3.2149, "PPC": 2.50431, "PPCK": 0.0, "PPS": 0.0, "PTAr": 0.0, "PYK": 1.75818, "PYRt2": 0.0, "RPE": 2.67848, "RPI": -2.2815, "SUCCt2_2": 0.0, "SUCCt3": 0.0, "SUCDi": 5.06438, "SUCOAS": -5.06438, "TALA": 1.49698, "THD2": 0.0, "TKT1": 1.49698, "TKT2": 1.1815, "TPI": 7.47738}}
\ No newline at end of file
+{"maximum": {"PTAr": 0.0, "ACONTb": 6.00725, "ACALDt": 0.0, "CYTBD": 43.59899, "FORt2": 0.0, "SUCDi": 30.98548, "GLCpts": 10.0, "EX_etoh_e": 0.0, "ACONTa": 6.00725, "EX_nh4_e": -4.76532, "FBA": 7.47738, "CS": 6.00725, "FRD7": 25.9211, "MDH": 5.06438, "GLNS": 0.22346, "ADK1": 0.0, "PPCK": 0.0, "EX_glu__L_e": 0.0, "FRUpts2": 0.0, "EX_co2_e": 22.80983, "EX_acald_e": 0.0, "EX_pyr_e": 0.0, "FUMt2_2": 0.0, "THD2": 0.0, "O2t": 21.79949, "RPE": 2.67848, "GAPD": 16.02353, "GND": 4.95998, "ICDHyr": 6.00725, "NH4t": 4.76532, "PGK": -16.02353, "EX_fru_e": 0.0, "PGM": -14.71614, "EX_glc__D_e": -10.0, "CO2t": -22.80983, "MALt2_2": 0.0, "EX_mal__L_e": 0.0, "ATPS4r": 45.51401, "FBP": 0.0, "EX_succ_e": 0.0, "D_LACt2": 0.0, "PYRt2": 0.0, "NADTRHD": 0.0, "TALA": 1.49698, "ATPM": 8.39, "GLUN": 0.0, "AKGDH": 5.06438, "PPC": 2.50431, "PIt2r": 3.2149, "G6PDH2r": 4.95998, "ETOHt2r": 0.0, "EX_ac_e": 0.0, "GLUDy": -4.54186, "ACALD": 0.0, "EX_pi_e": -3.2149, "TKT1": 1.49698, "ACKr": 0.0, "PDH": 9.28253, "EX_for_e": 0.0, "ME2": 0.0, "PGL": 4.95998, "ENO": 14.71614, "TKT2": 1.1815, "EX_h_e": 17.53087, "GLUt2r": 0.0, "SUCOAS": -5.06438, "FUM": 5.06438, "PGI": 4.86086, "ALCD2x": 0.0, "EX_gln__L_e": 0.0, "FORti": 0.0, "ICL": 0.0, "ME1": 0.0, "GLNabc": 0.0, "PFK": 7.47738, "EX_akg_e": 0.0, "NADH16": 38.53461, "PYK": 1.75818, "AKGt2r": 0.0, "MALS": 0.0, "EX_o2_e": -21.79949, "SUCCt3": 0.0, "RPI": -2.2815, "EX_lac__D_e": 0.0, "PFL": 0.0, "ACt2r": 0.0, "H2Ot": -29.17583, "LDH_D": 0.0, "SUCCt2_2": 0.0, "TPI": 7.47738, "GLUSy": 0.0, "Biomass_Ecoli_core": 0.87392, "PPS": 0.0, "EX_h2o_e": 29.17583, "EX_fum_e": 0.0}, "minimum": {"PTAr": 0.0, "ACONTb": 6.00725, "ACALDt": -0.0, "CYTBD": 43.59899, "FORt2": 0.0, "SUCDi": 5.06438, "GLCpts": 10.0, "EX_etoh_e": 0.0, "ACONTa": 6.00725, "EX_nh4_e": -4.76532, "FBA": 7.47738, "CS": 6.00725, "FRD7": 0.0, "MDH": 5.06438, "GLNS": 0.22346, "ADK1": 0.0, "PPCK": 0.0, "EX_glu__L_e": 0.0, "FRUpts2": -0.0, "EX_co2_e": 22.80983, "EX_acald_e": 0.0, "EX_pyr_e": 0.0, "FUMt2_2": -0.0, "THD2": 0.0, "O2t": 21.79949, "RPE": 2.67848, "GAPD": 16.02353, "GND": 4.95998, "ICDHyr": 6.00725, "NH4t": 4.76532, "PGK": -16.02353, "EX_fru_e": 0.0, "PGM": -14.71614, "EX_glc__D_e": -10.0, "CO2t": -22.80983, "MALt2_2": -0.0, "EX_mal__L_e": 0.0, "ATPS4r": 45.51401, "FBP": 0.0, "EX_succ_e": 0.0, "D_LACt2": -0.0, "PYRt2": -0.0, "NADTRHD": 0.0, "TALA": 1.49698, "ATPM": 8.39, "GLUN": 0.0, "AKGDH": 5.06438, "PPC": 2.50431, "PIt2r": 3.2149, "G6PDH2r": 4.95998, "ETOHt2r": -0.0, "EX_ac_e": 0.0, "GLUDy": -4.54186, "ACALD": -0.0, "EX_pi_e": -3.2149, "TKT1": 1.49698, "ACKr": -0.0, "PDH": 9.28253, "EX_for_e": 0.0, "ME2": 0.0, "PGL": 4.95998, "ENO": 14.71614, "TKT2": 1.1815, "EX_h_e": 17.53087, "GLUt2r": -0.0, "SUCOAS": -5.06438, "FUM": 5.06438, "PGI": 4.86086, "ALCD2x": -0.0, "EX_gln__L_e": 0.0, "FORti": 0.0, "ICL": 0.0, "ME1": 0.0, "GLNabc": 0.0, "PFK": 7.47738, "EX_akg_e": 0.0, "NADH16": 38.53461, "PYK": 1.75818, "AKGt2r": -0.0, "MALS": 0.0, "EX_o2_e": -21.79949, "SUCCt3": 0.0, "RPI": -2.2815, "EX_lac__D_e": 0.0, "PFL": 0.0, "ACt2r": -0.0, "H2Ot": -29.17583, "LDH_D": -0.0, "SUCCt2_2": 0.0, "TPI": 7.47738, "GLUSy": 0.0, "Biomass_Ecoli_core": 0.87392, "PPS": 0.0, "EX_h2o_e": 29.17583, "EX_fum_e": 0.0}}
\ No newline at end of file
diff --git a/cobra/test/data/textbook_solution.pickle b/cobra/test/data/textbook_solution.pickle
index a9633a37a..8bebff339 100644
Binary files a/cobra/test/data/textbook_solution.pickle and b/cobra/test/data/textbook_solution.pickle differ
diff --git a/cobra/test/test_model.py b/cobra/test/test_model.py
index ecb98fd68..9460374a5 100644
--- a/cobra/test/test_model.py
+++ b/cobra/test/test_model.py
@@ -12,7 +12,7 @@
 from optlang.symbolics import Zero
 
 import cobra.util.solver as su
-from cobra.core import Metabolite, Model, Reaction
+from cobra.core import Metabolite, Model, Reaction, Group
 from cobra.util import create_stoichiometric_matrix
 from cobra.exceptions import OptimizationError
 
@@ -238,13 +238,14 @@ def test_copy(self, model):
         assert PGI is not copied
         assert PGI._model is model
         assert copied._model is not model
-        # the copy should refer to different metabolites and genes
+        # the copy should refer to different metabolites, genes, and groups
         for met in copied.metabolites:
             assert met is not model.metabolites.get_by_id(met.id)
             assert met.model is not model
         for gene in copied.genes:
             assert gene is not model.genes.get_by_id(gene.id)
             assert gene.model is not model
+        assert len(model.get_associated_groups(copied.id)) == 0
 
     def test_iadd(self, model):
         PGI = model.reactions.PGI
@@ -331,6 +332,33 @@ def test_repr_html_(self, model):
         assert '<table>' in model.genes[0]._repr_html_()
 
 
+class TestCobraGroups:
+    def test_group_add_elements(self, model):
+        num_members = 5
+        reactions_for_group = model.reactions[0:num_members]
+        group = Group("arbitrary_group1")
+        group.add_members(reactions_for_group)
+        group.kind = "collection"
+        # number of member sin group should equal the number of reactions
+        # assigned to the group
+        assert len(group.members) == num_members
+
+        # Choose an overlapping, larger subset of reactions for the group
+        num_total_members = 12
+        reactions_for_larger_group = model.reactions[0:num_total_members]
+        group.add_members(reactions_for_larger_group)
+        assert len(group.members) == num_total_members
+
+    def test_group_kind(self):
+        group = Group("arbitrary_group1")
+        with pytest.raises(ValueError) as excinfo:
+            group.kind = "non-SBML compliant group kind"
+        assert "Kind can only by one of:" in str(excinfo.value)
+
+        group.kind = "collection"
+        assert group.kind == "collection"
+
+
 class TestCobraModel:
     """test core cobra functions"""
 
@@ -602,6 +630,70 @@ def test_remove_gene(self, model):
         for reaction in gene_reactions:
             assert target_gene not in reaction.genes
 
+    def test_group_model_reaction_association(self, model):
+        num_members = 5
+        reactions_for_group = model.reactions[0:num_members]
+        group = Group("arbitrary_group1")
+        group.add_members(reactions_for_group)
+        group.kind = "collection"
+        model.add_groups([group])
+        # group should point to and be associated with the model
+        assert group._model is model
+        assert group in model.groups
+
+        # model.get_associated_groups should find the group for each reaction
+        # we added to the group
+        for reaction in reactions_for_group:
+            assert group in model.get_associated_groups(reaction)
+
+        # remove the group from the model and check that reactions are no
+        # longer associated with the group
+        model.remove_groups([group])
+        assert group not in model.groups
+        assert group._model is not model
+        for reaction in reactions_for_group:
+            assert group not in model.get_associated_groups(reaction)
+
+    def test_group_members_add_to_model(self, model):
+        # remove a few reactions from the model and add them to a new group
+        num_members = 5
+        reactions_for_group = model.reactions[0:num_members]
+        model.remove_reactions(reactions_for_group, remove_orphans=False)
+        group = Group("arbitrary_group1")
+        group.add_members(reactions_for_group)
+        group.kind = "collection"
+        # the old reactions should not be in the model
+        for reaction in reactions_for_group:
+            assert reaction not in model.reactions
+
+        # add the group to the model and check that the reactions were added
+        model.add_groups([group])
+        assert group in model.groups
+        for reaction in reactions_for_group:
+            assert reaction in model.reactions
+
+    def test_group_loss_of_elements(self, model):
+        # when a metabolite, reaction, or gene is removed from a model, it
+        # should no longer be a member of any groups
+        num_members_each = 5
+        elements_for_group = model.reactions[0:num_members_each]
+        elements_for_group.extend(model.metabolites[0:num_members_each])
+        elements_for_group.extend(model.genes[0:num_members_each])
+        group = Group("arbitrary_group1")
+        group.add_members(elements_for_group)
+        group.kind = "collection"
+        model.add_groups([group])
+
+        remove_met = model.metabolites[0]
+        model.remove_metabolites([remove_met])
+        remove_rxn = model.reactions[0]
+        model.remove_reactions([remove_rxn])
+        remove_gene = model.genes[0]
+        remove_gene.remove_from_model()
+        assert remove_met not in group.members
+        assert remove_rxn not in group.members
+        assert remove_gene not in group.members
+
     def test_exchange_reactions(self, model):
         assert set(model.exchanges) == set([rxn for rxn in model.reactions
                                             if rxn.id.startswith("EX")])