[med-svn] [python-cobra] 02/05: Imported Upstream version 0.4.0b6

Afif Elghraoui afif-guest at moszumanska.debian.org
Tue Dec 22 06:21:24 UTC 2015


This is an automated email from the git hooks/post-receive script.

afif-guest pushed a commit to branch master
in repository python-cobra.

commit 24c228544e0d6649fc5167d7ad241c02235594fb
Author: Afif Elghraoui <afif at ghraoui.name>
Date:   Mon Dec 21 20:10:17 2015 -0800

    Imported Upstream version 0.4.0b6
---
 .gitignore                                         |   1 +
 .travis.yml                                        |  16 +-
 CONTRIBUTING.md                                    |  25 ++
 INSTALL.md                                         |   3 +-
 README.md                                          |   2 +-
 appveyor.yml                                       |  23 +-
 appveyor/run_with_compiler.cmd                     |  47 ---
 appveyor/run_with_env.cmd                          |  88 +++++
 cobra/VERSION                                      |   2 +-
 cobra/__init__.py                                  |   2 +-
 cobra/core/DictList.py                             |   7 +-
 cobra/core/Gene.py                                 |   2 +-
 cobra/core/Metabolite.py                           |   2 +-
 cobra/core/Model.py                                |   1 +
 cobra/core/Object.py                               |   2 +-
 cobra/core/Reaction.py                             |  13 +-
 cobra/design/__init__.py                           |   1 +
 cobra/design/design_algorithms.py                  | 400 +++++++++++++++++++++
 cobra/flux_analysis/phenotype_phase_plane.py       |  15 +-
 cobra/flux_analysis/reaction.py                    |   5 +-
 cobra/io/json.py                                   |   4 +-
 cobra/io/mat.py                                    |  14 +-
 cobra/io/sbml.py                                   |  68 +++-
 cobra/io/sbml3.py                                  | 110 +++---
 cobra/manipulation/__init__.py                     |   9 +-
 cobra/manipulation/annotate.py                     |  24 ++
 cobra/manipulation/modify.py                       | 110 +++++-
 cobra/manipulation/validate.py                     |  54 +++
 cobra/solvers/__init__.py                          |   6 +-
 cobra/solvers/cglpk.pyx                            |  11 +-
 cobra/solvers/cplex_solver.py                      |   3 +-
 cobra/solvers/cplex_solver_java.py                 |  16 +-
 cobra/solvers/glpk_solver.py                       |  13 +-
 cobra/solvers/glpk_solver_java.py                  |  30 +-
 cobra/solvers/gurobi_solver.py                     |  36 +-
 cobra/solvers/gurobi_solver_java.py                |  14 +-
 cobra/solvers/mosek.py                             |   6 +
 cobra/test/__init__.py                             |   3 +-
 cobra/test/data/invalid0.xml                       |  47 +++
 cobra/test/data/invalid1.xml                       |  45 +++
 cobra/test/data/invalid2.xml                       |  68 ++++
 cobra/test/data/mini.json                          |  24 +-
 cobra/test/data/mini.mat                           | Bin 14352 -> 14600 bytes
 cobra/test/data/mini.pickle                        | Bin 25348 -> 25060 bytes
 cobra/test/data/mini_cobra.xml                     | 212 ++++++-----
 cobra/test/data/mini_fbc1.xml                      | 218 +++++------
 cobra/test/data/mini_fbc2.xml                      |  44 +--
 cobra/test/data/mini_fbc2.xml.bz2                  | Bin 4939 -> 4957 bytes
 cobra/test/data/mini_fbc2.xml.gz                   | Bin 5637 -> 5670 bytes
 cobra/test/data/textbook.xml.gz                    | Bin 16329 -> 16890 bytes
 cobra/test/design.py                               |  89 +++++
 cobra/test/flux_analysis.py                        | 157 +-------
 cobra/test/io_tests.py                             |  37 +-
 cobra/test/manipulation.py                         | 250 +++++++++++++
 cobra/test/unit_tests.py                           |  36 ++
 cobra/test_all.py                                  |   7 -
 cobra/topology/reporter_metabolites.py             |   4 +-
 documentation_builder/io.ipynb                     |  12 +-
 documentation_builder/io.rst                       |  41 +--
 documentation_builder/loopless.ipynb               |   4 +-
 documentation_builder/loopless.rst                 |   4 +-
 documentation_builder/phenotype_phase_plane.ipynb  |  23 +-
 documentation_builder/phenotype_phase_plane.rst    |   7 +-
 .../phenotype_phase_plane_3_0.png                  | Bin 94579 -> 89909 bytes
 .../phenotype_phase_plane_5_0.png                  | Bin 106030 -> 100924 bytes
 .../phenotype_phase_plane_5_1.png                  | Bin 97694 -> 92540 bytes
 .../phenotype_phase_plane_7_0.png                  | Bin 68889 -> 63193 bytes
 documentation_builder/simulating.ipynb             |  80 ++++-
 documentation_builder/simulating.rst               |  43 ++-
 documentation_builder/solvers.ipynb                |   2 +-
 setup.py                                           |  20 +-
 71 files changed, 1986 insertions(+), 676 deletions(-)

diff --git a/.gitignore b/.gitignore
index dd122dc..7357a55 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,7 @@ cobra.egg-info
 build/
 dist/
 .DS_Store
+.eggs/
 setuptools-*egg
 setuptools-*.tar.gz
 cobra/solvers/cglpk.c
diff --git a/.travis.yml b/.travis.yml
index 3b51629..18c5fa3 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -6,6 +6,7 @@ cache:
 python:
   - "2.7"
   - "3.4"
+  - "3.5"
 addons:
   apt:
     packages:
@@ -22,19 +23,22 @@ env:
   - PIP_CACHE_DIR=$HOME/.cache/pip
 before_install:
   - pip install pip --upgrade
-  - export PIP_OPTS="-f https://opencobra.github.io/pypi_cobrapy_travis --no-index"
   # These get cached
   - pip install numpy scipy python-libsbml -v  # verbose for long-running compiles
-  - pip install cython coveralls jsonschema six
-  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then pip install lxml; fi
-  # These use pre-compiled wheels (defined in PIP_OPTS)
-  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then pip install glpk cylp $PIP_OPTS; fi
+  - pip install cython coveralls jsonschema six matplotlib
+  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then pip install lxml glpk pep8 palettable; fi
   # Download esolver and add it to the path
   - wget https://opencobra.github.io/pypi_cobrapy_travis/esolver.gz
   - gzip -d esolver.gz; chmod +x esolver; export PATH=$PATH:$PWD
+  - mkdir -p $HOME/.config/matplotlib
+  - "echo 'backend: Agg' >> $HOME/.config/matplotlib/matplotlibrc"
 install:
   - python setup.py develop
 # # command to run tests
-script: coverage run --source=cobra setup.py test
+script:
+  - coverage run --source=cobra setup.py test
+  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then
+    pep8 cobra/core cobra/design cobra/manipulation cobra/test cobra/version.py;
+    fi
 after_success:
   - coveralls
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000..5724522
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,25 @@
+Contribution Guidelines
+-----------------------
+
+Generally, the following practices are recommended for making contributions to
+cobrapy. These aren't all necessarily hard-and-fast rules, but should serve as
+guidelines in most cases.
+
+1. Please comment code.
+2. All new python code should be pep8 compliant.
+3. Please use git best practices, with a 50 line summary for each commit.
+   Generally, separate features should be made in separate commits so
+   they can be tested and merged independently. For example, adding a new
+   solver would be a separate commit from fixing whitespace in cobra.core.
+4. Documentation is written as IPython/jupyter notebooks in the
+   ```documentation_builder``` directory, which are then converted to
+   rst by the ```autodoc.sh``` script.
+5. Tests are in the ```cobra/test``` directory. They are automatically run
+   through continuous integration services on both python 2 and python 3
+   when pull requests are made.
+6. Please write tests for new functions. Writing documentation as well
+   would also be very helpful.
+7. Ensure code will work with both python 2 and python 3. For example,
+   instead of ```my_dict.iteritems()``` use ```six.iteritems(my_dict)```
+
+Thank you very much for contributing to cobrapy.
diff --git a/INSTALL.md b/INSTALL.md
index 8839752..125b291 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -45,7 +45,8 @@ index](https://pypi.python.org/pypi/cobra/).
 
 ## Hacking version installation
 Use pip to install [Cython](http://cython.org/). Install libglpk 
-using your package manger. This would be ```brew install glpk``` on a Mac
+using your package manger. This would be
+```brew install homebrew/science/glpk``` on a Mac
 and ```sudo apt-get install libglpk-dev``` on debian-based systems
 (including Ubuntu and Mint). GLPK can also be compiled from the
 released source.
diff --git a/README.md b/README.md
index 7810886..5b96ea6 100644
--- a/README.md
+++ b/README.md
@@ -16,7 +16,7 @@ flux variability analysis, and gene deletion analyses.
 To install, please follow the [instructions](INSTALL.md).
 
 The documentation is browseable online at
-[readthedocs](https://cobrapy.readthedocs.org/en/latest/)
+[readthedocs](https://cobrapy.readthedocs.org/en/stable/)
 and can also be
 [downloaded](https://readthedocs.org/projects/cobrapy/downloads/).
 
diff --git a/appveyor.yml b/appveyor.yml
index c0feca7..bbd790c 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -4,11 +4,12 @@ environment:
     # SDK v7.0 MSVC Express 2008's SetEnv.cmd script will fail if the
     # /E:ON and /V:ON options are not enabled in the batch script intepreter
     # See: http://stackoverflow.com/a/13751649/163740
-    WITH_COMPILER: "cmd /E:ON /V:ON /C .\\appveyor\\run_with_compiler.cmd"
+    WITH_COMPILER: "cmd /E:ON /V:ON /C .\\appveyor\\run_with_env.cmd"
+    UCRT: " "
 
   matrix:
     - PYTHON: "C:\\Python27"
-      PYTHON_VERSION: "2.7.9"
+      PYTHON_VERSION: "2.7.10"
       PYTHON_ARCH: "32"
 
     - PYTHON: "C:\\Python34"
@@ -16,27 +17,30 @@ environment:
       PYTHON_ARCH: "32"
 
     - PYTHON: "C:\\Python27-x64"
-      PYTHON_VERSION: "2.7.9"
+      PYTHON_VERSION: "2.7.10"
       PYTHON_ARCH: "64"
-      WINDOWS_SDK_VERSION: "v7.0"
 
     - PYTHON: "C:\\Python34-x64"
       PYTHON_VERSION: "3.4.3"
       PYTHON_ARCH: "64"
-      WINDOWS_SDK_VERSION: "v7.1"
+
+    - PYTHON: "C:\\Python35-x64"
+      PYTHON_VERSION: "3.5.0"
+      PYTHON_ARCH: "64"
+      UCRT: "UCRT"
 
 clone_depth: 25
 
 init:
-  - "ECHO %PYTHON% %PYTHON_VERSION% %PYTHON_ARCH%"
+  - "ECHO %PYTHON% %PYTHON_VERSION% %PYTHON_ARCH% %UCRT%"
 
 install:
   - "powershell appveyor\\install.ps1"
   - ps: Start-FileDownload 'https://bitbucket.org/gutworth/six/raw/default/six.py'
   # Download glpk.h and glpk.lib (separate files for 32/64 bit)
   - ps: Start-FileDownload 'https://opencobra.github.io/pypi_cobrapy_travis/glpk.h'
-  - "appveyor DownloadFile https://opencobra.github.io/pypi_cobrapy_travis/glpk.lib%PYTHON_ARCH% -FileName glpk.lib"
-  - "%PYTHON%/python -m pip install Cython -f https://opencobra.github.io/pypi_cobrapy_travis/ --no-index"
+  - "appveyor DownloadFile https://opencobra.github.io/pypi_cobrapy_travis/glpk.lib%PYTHON_ARCH%%UCRT% -FileName glpk.lib"
+  - "%PYTHON%/python -m pip install Cython"
 
 build: off
 
@@ -44,8 +48,7 @@ test_script:
   - "%WITH_COMPILER% %PYTHON%/python setup.py test"
 
 after_test:
-  - "%WITH_COMPILER% %PYTHON%/python setup.py bdist_wheel"
-  - "%WITH_COMPILER% %PYTHON%/python setup.py bdist_wininst"
+  - "%WITH_COMPILER% %PYTHON%/python setup.py bdist_wheel bdist_wininst"
 
 artifacts:
   - path: dist\*
diff --git a/appveyor/run_with_compiler.cmd b/appveyor/run_with_compiler.cmd
deleted file mode 100644
index 3a472bc..0000000
--- a/appveyor/run_with_compiler.cmd
+++ /dev/null
@@ -1,47 +0,0 @@
-:: To build extensions for 64 bit Python 3, we need to configure environment
-:: variables to use the MSVC 2010 C++ compilers from GRMSDKX_EN_DVD.iso of:
-:: MS Windows SDK for Windows 7 and .NET Framework 4 (SDK v7.1)
-::
-:: To build extensions for 64 bit Python 2, we need to configure environment
-:: variables to use the MSVC 2008 C++ compilers from GRMSDKX_EN_DVD.iso of:
-:: MS Windows SDK for Windows 7 and .NET Framework 3.5 (SDK v7.0)
-::
-:: 32 bit builds do not require specific environment configurations.
-::
-:: Note: this script needs to be run with the /E:ON and /V:ON flags for the
-:: cmd interpreter, at least for (SDK v7.0)
-::
-:: More details at:
-:: https://github.com/cython/cython/wiki/64BitCythonExtensionsOnWindows
-:: http://stackoverflow.com/a/13751649/163740
-::
-:: Author: Olivier Grisel
-:: License: CC0 1.0 Universal: http://creativecommons.org/publicdomain/zero/1.0/
- at ECHO OFF
-
-SET COMMAND_TO_RUN=%*
-SET WIN_SDK_ROOT=C:\Program Files\Microsoft SDKs\Windows
-
-SET MAJOR_PYTHON_VERSION="%PYTHON_VERSION:~0,1%"
-IF %MAJOR_PYTHON_VERSION% == "2" (
-    SET WINDOWS_SDK_VERSION="v7.0"
-) ELSE IF %MAJOR_PYTHON_VERSION% == "3" (
-    SET WINDOWS_SDK_VERSION="v7.1"
-) ELSE (
-    ECHO Unsupported Python version: "%MAJOR_PYTHON_VERSION%"
-    EXIT 1
-)
-
-IF "%PYTHON_ARCH%"=="64" (
-    ECHO Configuring Windows SDK %WINDOWS_SDK_VERSION% for Python %MAJOR_PYTHON_VERSION% on a 64 bit architecture
-    SET DISTUTILS_USE_SDK=1
-    SET MSSdk=1
-    "%WIN_SDK_ROOT%\%WINDOWS_SDK_VERSION%\Setup\WindowsSdkVer.exe" -q -version:%WINDOWS_SDK_VERSION%
-    "%WIN_SDK_ROOT%\%WINDOWS_SDK_VERSION%\Bin\SetEnv.cmd" /x64 /release
-    ECHO Executing: %COMMAND_TO_RUN%
-    call %COMMAND_TO_RUN% || EXIT 1
-) ELSE (
-    ECHO Using default MSVC build environment for 32 bit architecture
-    ECHO Executing: %COMMAND_TO_RUN%
-    call %COMMAND_TO_RUN% || EXIT 1
-)
diff --git a/appveyor/run_with_env.cmd b/appveyor/run_with_env.cmd
new file mode 100644
index 0000000..5da547c
--- /dev/null
+++ b/appveyor/run_with_env.cmd
@@ -0,0 +1,88 @@
+:: To build extensions for 64 bit Python 3, we need to configure environment
+:: variables to use the MSVC 2010 C++ compilers from GRMSDKX_EN_DVD.iso of:
+:: MS Windows SDK for Windows 7 and .NET Framework 4 (SDK v7.1)
+::
+:: To build extensions for 64 bit Python 2, we need to configure environment
+:: variables to use the MSVC 2008 C++ compilers from GRMSDKX_EN_DVD.iso of:
+:: MS Windows SDK for Windows 7 and .NET Framework 3.5 (SDK v7.0)
+::
+:: 32 bit builds, and 64-bit builds for 3.5 and beyond, do not require specific
+:: environment configurations.
+::
+:: Note: this script needs to be run with the /E:ON and /V:ON flags for the
+:: cmd interpreter, at least for (SDK v7.0)
+::
+:: More details at:
+:: https://github.com/cython/cython/wiki/64BitCythonExtensionsOnWindows
+:: http://stackoverflow.com/a/13751649/163740
+::
+:: Author: Olivier Grisel
+:: License: CC0 1.0 Universal: http://creativecommons.org/publicdomain/zero/1.0/
+::
+:: Notes about batch files for Python people:
+::
+:: Quotes in values are literally part of the values:
+::      SET FOO="bar"
+:: FOO is now five characters long: " b a r "
+:: If you don't want quotes, don't include them on the right-hand side.
+::
+:: The CALL lines at the end of this file look redundant, but if you move them
+:: outside of the IF clauses, they do not run properly in the SET_SDK_64==Y
+:: case, I don't know why.
+ at ECHO OFF
+
+SET COMMAND_TO_RUN=%*
+SET WIN_SDK_ROOT=C:\Program Files\Microsoft SDKs\Windows
+SET WIN_WDK=c:\Program Files (x86)\Windows Kits\10\Include\wdf
+
+:: Extract the major and minor versions, and allow for the minor version to be
+:: more than 9.  This requires the version number to have two dots in it.
+SET MAJOR_PYTHON_VERSION=%PYTHON_VERSION:~0,1%
+IF "%PYTHON_VERSION:~3,1%" == "." (
+    SET MINOR_PYTHON_VERSION=%PYTHON_VERSION:~2,1%
+) ELSE (
+    SET MINOR_PYTHON_VERSION=%PYTHON_VERSION:~2,2%
+)
+
+:: Based on the Python version, determine what SDK version to use, and whether
+:: to set the SDK for 64-bit.
+IF %MAJOR_PYTHON_VERSION% == 2 (
+    SET WINDOWS_SDK_VERSION="v7.0"
+    SET SET_SDK_64=Y
+) ELSE (
+    IF %MAJOR_PYTHON_VERSION% == 3 (
+        SET WINDOWS_SDK_VERSION="v7.1"
+        IF %MINOR_PYTHON_VERSION% LEQ 4 (
+            SET SET_SDK_64=Y
+        ) ELSE (
+            SET SET_SDK_64=N
+            IF EXIST "%WIN_WDK%" (
+                :: See: https://connect.microsoft.com/VisualStudio/feedback/details/1610302/
+                REN "%WIN_WDK%" 0wdf
+            )
+        )
+    ) ELSE (
+        ECHO Unsupported Python version: "%MAJOR_PYTHON_VERSION%"
+        EXIT 1
+    )
+)
+
+IF %PYTHON_ARCH% == 64 (
+    IF %SET_SDK_64% == Y (
+        ECHO Configuring Windows SDK %WINDOWS_SDK_VERSION% for Python %MAJOR_PYTHON_VERSION% on a 64 bit architecture
+        SET DISTUTILS_USE_SDK=1
+        SET MSSdk=1
+        "%WIN_SDK_ROOT%\%WINDOWS_SDK_VERSION%\Setup\WindowsSdkVer.exe" -q -version:%WINDOWS_SDK_VERSION%
+        "%WIN_SDK_ROOT%\%WINDOWS_SDK_VERSION%\Bin\SetEnv.cmd" /x64 /release
+        ECHO Executing: %COMMAND_TO_RUN%
+        call %COMMAND_TO_RUN% || EXIT 1
+    ) ELSE (
+        ECHO Using default MSVC build environment for 64 bit architecture
+        ECHO Executing: %COMMAND_TO_RUN%
+        call %COMMAND_TO_RUN% || EXIT 1
+    )
+) ELSE (
+    ECHO Using default MSVC build environment for 32 bit architecture
+    ECHO Executing: %COMMAND_TO_RUN%
+    call %COMMAND_TO_RUN% || EXIT 1
+)
diff --git a/cobra/VERSION b/cobra/VERSION
index d4e083c..c259ae3 100644
--- a/cobra/VERSION
+++ b/cobra/VERSION
@@ -1 +1 @@
-0.4.0b4
+0.4.0b6
diff --git a/cobra/__init__.py b/cobra/__init__.py
index 4df31aa..e57f01a 100644
--- a/cobra/__init__.py
+++ b/cobra/__init__.py
@@ -20,7 +20,7 @@ from .version import get_version
 __version__ = get_version()
 from .core import Object, Metabolite, Gene, Reaction, Model, \
     DictList, Species
-from . import io, flux_analysis
+from . import io, flux_analysis, design
 
 try:
     from .core import ArrayBasedModel
diff --git a/cobra/core/DictList.py b/cobra/core/DictList.py
index c61dc46..70619a5 100644
--- a/cobra/core/DictList.py
+++ b/cobra/core/DictList.py
@@ -1,6 +1,6 @@
 from copy import copy, deepcopy
 import re
-from six import string_types, iteritems
+from six import string_types, iteritems, PY3
 from itertools import islice
 
 try:
@@ -288,7 +288,10 @@ class DictList(list):
         if key is None:
             def key(i):
                 return i.id
-        list.sort(self, cmp=cmp, key=key, reverse=reverse)
+        if PY3:
+            list.sort(self, key=key, reverse=reverse)
+        else:
+            list.sort(self, cmp=cmp, key=key, reverse=reverse)
         self._generate_index()
 
     def __getitem__(self, i):
diff --git a/cobra/core/Gene.py b/cobra/core/Gene.py
index 9a608e8..7e97534 100644
--- a/cobra/core/Gene.py
+++ b/cobra/core/Gene.py
@@ -131,7 +131,7 @@ def parse_gpr(str_expr):
 
 class Gene(Species):
 
-    def __init__(self, id=id, name=None, functional=True):
+    def __init__(self, id=None, name="", functional=True):
         """
         id: A string.
 
diff --git a/cobra/core/Metabolite.py b/cobra/core/Metabolite.py
index 1a3b139..e1a46fb 100644
--- a/cobra/core/Metabolite.py
+++ b/cobra/core/Metabolite.py
@@ -16,7 +16,7 @@ class Metabolite(Species):
     """
 
     def __init__(self, id=None, formula=None,
-                 name=None, compartment=None):
+                 name="", compartment=None):
         """
         id: str
 
diff --git a/cobra/core/Model.py b/cobra/core/Model.py
index 9ab90db..2a62c4f 100644
--- a/cobra/core/Model.py
+++ b/cobra/core/Model.py
@@ -47,6 +47,7 @@ class Model(Object):
             # genes based on their ids {Gene.id: Gene}
             self.compartments = {}
             self.solution = Solution(None)
+            self.media_compositions = {}
 
     @property
     def description(self):
diff --git a/cobra/core/Object.py b/cobra/core/Object.py
index 03ce0bb..667c4e1 100644
--- a/cobra/core/Object.py
+++ b/cobra/core/Object.py
@@ -4,7 +4,7 @@ from six import iteritems
 class Object(object):
     """Defines common behavior of object in cobra.core"""
 
-    def __init__(self, id=None, name=None):
+    def __init__(self, id=None, name=""):
         """
         id: None or a string
 
diff --git a/cobra/core/Reaction.py b/cobra/core/Reaction.py
index 96346d8..119db1f 100644
--- a/cobra/core/Reaction.py
+++ b/cobra/core/Reaction.py
@@ -47,26 +47,27 @@ class Reaction(Object):
 
     """
 
-    def __init__(self, name=None):
+    def __init__(self, id=None, name='', subsystem='', lower_bound=0.,
+                 upper_bound=1000., objective_coefficient=0.):
         """An object for housing reactions and associated information
         for cobra modeling.
 
         """
-        Object.__init__(self, name)
+        Object.__init__(self, id, name)
         self._gene_reaction_rule = ''
-        self.subsystem = ''
+        self.subsystem = subsystem
         # The cobra.Genes that are used to catalyze the reaction
         self._genes = set()
         # A dictionary of metabolites and their stoichiometric coefficients in
         # this reaction.
         self._metabolites = {}
-        self.name = name
         # self.model is None or refers to the cobra.Model that
         # contains self
         self._model = None
 
-        self.objective_coefficient = self.lower_bound = 0.
-        self.upper_bound = 1000.
+        self.objective_coefficient = objective_coefficient
+        self.upper_bound = upper_bound
+        self.lower_bound = lower_bound
         # Used during optimization.  Indicates whether the
         # variable is modeled as continuous, integer, binary, semicontinous, or
         # semiinteger.
diff --git a/cobra/design/__init__.py b/cobra/design/__init__.py
new file mode 100644
index 0000000..fa6a8c1
--- /dev/null
+++ b/cobra/design/__init__.py
@@ -0,0 +1 @@
+from .design_algorithms import *
diff --git a/cobra/design/design_algorithms.py b/cobra/design/design_algorithms.py
new file mode 100644
index 0000000..4a6c345
--- /dev/null
+++ b/cobra/design/design_algorithms.py
@@ -0,0 +1,400 @@
+from ..core import Model, Reaction, Metabolite
+from ..manipulation.modify import convert_to_irreversible, canonical_form
+
+from six import iteritems
+from collections import defaultdict
+from itertools import chain
+from copy import deepcopy
+
+
+def _add_decision_variable(model, reaction_id):
+    """Add an integer decision variable for the given reaction."""
+    reaction = model.reactions.get_by_id(reaction_id)
+    # add integer variable
+    var = Reaction("%s_decision_var" % reaction_id)
+    var.lower_bound = 0
+    var.upper_bound = 1
+    var.variable_kind = "integer"
+    var.decision_reaction_id = reaction_id
+    model.add_reaction(var)
+    # add constraints
+    # v <= ub * y  -->  v - ub * y <= 0
+    ub_constr = Metabolite("%s_upper_bound" % var.id)
+    ub_constr._constraint_sense = "L"
+    # v >= lb * y  -->  - v + lb * y <= 0
+    lb_constr = Metabolite("%s_lower_bound" % var.id)
+    lb_constr._constraint_sense = "L"
+    reaction.add_metabolites({lb_constr: - 1,
+                              ub_constr:   1})
+    var.add_metabolites({lb_constr:   reaction.lower_bound,
+                         ub_constr: - reaction.upper_bound})
+    return var
+
+
+def set_up_optknock(model, chemical_objective, knockable_reactions,
+                    biomass_objective=None, n_knockouts=5,
+                    n_knockouts_required=True, dual_maximum=1000, copy=True):
+    """Set up the OptKnock problem described by Burgard et al., 2003:
+
+        Burgard AP, Pharkya P, Maranas CD. Optknock: a bilevel programming
+        framework for identifying gene knockout strategies for microbial strain
+        optimization.  Biotechnol Bioeng. 2003;84(6):647-57.
+        http://dx.doi.org/10.1002/bit.10803.
+
+
+    model : :class:`~cobra.core.Model` object.
+
+    chemical_objective: str. The ID of the reaction to maximize in the outer
+    problem.
+
+    knockable_reactions: [str]. A list of reaction IDs that can be knocked out.
+
+    biomass_objective: str. The ID of the reaction to maximize in the inner
+    problem. By default, this is the existing objective function in the passed
+    model.
+
+    n_knockouts: int. The number of knockouts allowable.
+
+    n_knockouts_required: bool. Require exactly the number of knockouts
+    specified by n_knockouts.
+
+    dual_maximum: float or int. The upper bound for dual variables.
+
+    copy: bool. Copy the model before making any modifications.
+
+
+    Zachary King 2015
+
+    """
+
+    if copy:
+        model = model.copy()
+
+    # add the integer decision variables
+    decision_variable_ids = [_add_decision_variable(model, r_id).id
+                             for r_id in knockable_reactions]
+
+    # inner problem
+    inner_problem = model.copy()
+    if biomass_objective:
+        found = False
+        for reaction in inner_problem.reactions:
+            obj = reaction.id == biomass_objective
+            reaction.objective_coefficient = 1 if obj else 0
+            if obj:
+                found = True
+        if not found:
+            raise Exception("Could not find biomass_objective %s in model" %
+                            biomass_objective)
+
+    # dual of inner problem
+    inner_dual = dual_problem(inner_problem,
+                              integer_vars_to_maintain=decision_variable_ids,
+                              already_irreversible=False, copy=False,
+                              dual_maximum=dual_maximum)
+
+    # add constraints and variables from inner problem to outer problem
+    inner_objectives = {}
+    for reaction in inner_dual.reactions:
+        inner_objectives[reaction.id] = reaction.objective_coefficient
+        reaction.objective_coefficient = 0
+        if reaction.id in model.reactions:
+            existing_reaction = model.reactions.get_by_id(reaction.id)
+            for met, coeff in iteritems(reaction.metabolites):
+                if met.id in model.metabolites:
+                    existing_reaction.add_metabolites(
+                        {model.metabolites.get_by_id(met.id): coeff})
+                else:
+                    existing_reaction.add_metabolites({deepcopy(met): coeff})
+        else:
+            model.add_reaction(reaction)
+
+    # constraint to set outer and inner objectives equal, and set chemical
+    # objective
+    equal_objectives_constr = Metabolite("equal_objectives_constraint")
+    equal_objectives_constr._constraint_sense = "E"
+    equal_objectives_constr._bound = 0
+    for reaction in model.reactions:
+        if reaction.objective_coefficient != 0:
+            reaction.add_metabolites({equal_objectives_constr:
+                                      reaction.objective_coefficient})
+        inner_objective = inner_objectives.get(reaction.id, 0)
+        if inner_objective:
+            reaction.add_metabolites(
+                {equal_objectives_constr: - inner_objective})
+        # set chemical objective
+        reaction.objective_coefficient = 1 \
+            if reaction.id == chemical_objective else 0
+
+    # add the n_knockouts constraint
+    n_knockouts_constr = Metabolite("n_knockouts_constraint")
+    n_knockouts_constr._constraint_sense = "E" if n_knockouts_required else "G"
+    n_knockouts_constr._bound = len(decision_variable_ids) - n_knockouts
+    for r_id in decision_variable_ids:
+        reaction = model.reactions.get_by_id(r_id)
+        reaction.add_metabolites({n_knockouts_constr: 1})
+
+    return model
+
+
+def run_optknock(optknock_problem, solver=None, tolerance_integer=1e-9,
+                 **kwargs):
+    """Run the OptKnock problem created with set_up_optknock.
+
+
+    optknock_problem: :class:`~cobra.core.Model` object. The problem generated
+    by set_up_optknock.
+
+    solver: str. The name of the preferred solver.
+
+    tolerance_integer: float. The integer tolerance for the MILP.
+
+    **kwargs: Keyword arguments are passed to Model.optimize().
+
+
+    Zachary King 2015
+
+    """
+    solution = optknock_problem.optimize(solver=solver,
+                                         tolerance_integer=tolerance_integer,
+                                         **kwargs)
+    solution.knockouts = []
+    for reaction in optknock_problem.reactions:
+        if solution.x_dict.get(reaction.id, None) == 0:
+            r_id = getattr(reaction, "decision_reaction_id", None)
+            if r_id is not None:
+                solution.knockouts.append(r_id)
+    return solution
+
+
+# This function will generalize the set_up_optknock code to other MILPs:
+# def dual_embed(outer_model, inner_model, ..., objective_sense="maximize",
+#                integer_vars_to_maintain=[], already_irreversible=False,
+#                copy=True, dual_maximum=1000):
+#     """Embed the dual of the inner model within the outer model"""
+
+
+def dual_problem(model, objective_sense="maximize",
+                 integer_vars_to_maintain=[],
+                 already_irreversible=False, copy=True, dual_maximum=1000):
+    """Return a new model representing the dual of the model.
+
+    Make the problem irreversible, then take the dual. Convert the problem:
+
+        Maximize (c^T)x subject to Ax <= b, x >= 0
+
+    which is something like this in COBRApy:
+
+        Maximize sum(objective_coefficient_j * reaction_j for all j)
+            s.t.
+            sum(coefficient_i_j * reaction_j for all j) <= metabolite_bound_i
+            reaction_j <= upper_bound_j
+            reaction_j >= 0
+
+    to the problem:
+
+        Minimize (b^T)w subject to (A^T)w >= c, w >= 0
+
+    which is something like this in COBRApy (S matrix is m x n):
+
+        Minimize sum( metabolite_bound_i * dual_i   for all i ) +
+                 sum( upper_bound_j *      dual_m+j for all j ) +
+            s.t.
+             sum( coefficient_i_j * dual_i for all i ) +
+             sum( dual_2m+j' for all j' ) >= objective_coefficient_j
+            dual_k >= 0
+
+
+    Arguments
+    ---------
+
+    model : :class:`~cobra.core.Model` object.
+
+    objective_sense: str. The objective sense of the starting problem, either
+    'maximize' or 'minimize'. A minimization problems will be converted to a
+    maximization before taking the dual. This function always returns a
+    minimization problem.
+
+    iteger_vars_to_maintain: [str]. A list of IDs for Boolean integer variables
+    to be maintained in the dual problem. See 'Maintaining integer variables'
+    below for more details
+
+    already_irreversible: Boolean. If True, then do not convert the model to
+    irreversible.
+
+    copy: bool. If True, then make a copy of the model before modifying
+    it. This is not necessary if already_irreversible is True.
+
+    dual_maximum: float or int. The upper bound for dual variables.
+
+
+    Maintaining integer variables
+    -----------------------------
+
+    The argument integer_vars_to_maintain can be used to specify certin Boolean
+    integer variables that will be maintained in the dual problem. This makes
+    it possible to join outer and inner problems in a bi-level MILP. The method
+    for maintaining integer variables is described by Tepper and Shlomi, 2010:
+
+        Tepper N, Shlomi T. Predicting metabolic engineering knockout
+        strategies for chemical production: accounting for competing pathways.
+        Bioinformatics.  2010;26(4):536-43. doi:10.1093/bioinformatics/btp704.
+
+    In COBRApy, this roughly translates to transforming (decision variables p,
+    integer constraints o):
+
+        Maximize (c^T)x subject to (A_x)x + (A_y)y <= b, x >= 0
+
+        (1) Maximize sum(objective_coefficient_j * reaction_j for all j)
+                s.t.
+        (2)     sum(coeff_i_j * reaction_j for all j) +
+                sum(decision_coeff_i_j * decision_var_j for all j)
+                <= metabolite_bound_i
+        (3)     reaction_j <= upper_bound_j
+        (4)     reaction_j >= 0
+
+    to the problem:
+
+        Minimize (b - (A_y)y)^T w subject to (A_x^T)w >= c, w >= 0
+
+    which linearizes to (with auxiliary variables z):
+
+        Minimize (b^T)w - { ((A_y)y)^T w with yw --> z }
+        subject to (A_x^T)w >= c, linearization constraints, w >= 0
+          Linearization constraints: z <= w_max * y, z <= w,
+                                     z >= w - w_max * (1 - y), z >= 0
+
+        (5) Minimize sum( metabolite_bound_i *  dual_i            for all i ) +
+                      sum( upper_bound_j *      dual_m+j          for all j ) +
+                    - sum( decision_coeff_i_j * auxiliary_var_i_j
+                          for all combinations i, j )
+                s.t.
+        (6)   - sum( coefficient_i_j * dual_i for all i ) - dual_m+j
+              <= - objective_coefficient_j
+        (7)     auxiliary_var_i_j - dual_maximum * decision_var_j          <= 0
+        (8)     auxiliary_var_i_j - dual_i                                 <= 0
+        (9)   - auxiliary_var_i_j + dual_i + dual_maximum * decision_var_j
+              <= dual_maximum
+       (10)     dual_maximum >= dual_i            >= 0
+       (11)     dual_maximum >= dual_m+j          >= 0
+       (12)     dual_maximum >= auxiliary_var_i_j >= 0
+       (13)                1 >= decision_var_j    >= 0
+
+
+    Zachary King 2015
+
+    """
+
+    # convert to canonical form and copy
+    model = canonical_form(model, objective_sense=objective_sense,
+                           already_irreversible=already_irreversible,
+                           copy=copy)
+
+    # new model for the dual
+    dual = Model("%s_dual" % model.id)
+
+    # keep track of dual_i
+    dual_var_for_met = {}
+
+    # add dual variables for constraints. (2) --> dual_i
+    for metabolite in model.metabolites:
+        # add constraints based on metabolite constraint sense
+        if metabolite._constraint_sense != "L":
+            raise Exception("Not a less than or equal constraint: %s"
+                            % metabolite.id)
+
+        var = Reaction("%s__dual" % metabolite.id)
+        # Without auxiliary variables, the objective coefficient would include
+        # integer variables when present. However, we will separate out the
+        # integer parts into objective coefficients for auxiliary variables.
+        var.objective_coefficient = metabolite._bound  # (5)
+        # [dual_vars] >= 0
+        var.lower_bound = 0
+        var.upper_bound = dual_maximum
+        dual.add_reaction(var)
+        # remember
+        dual_var_for_met[metabolite.id] = var
+
+    # keep track of decision variables (integer_vars_to_maintain) as tuples:
+    # (reaction in dual problem, reaction in original problem)
+    integer_vars_added = []
+
+    # add constraints and upper bound variables
+    for reaction in model.reactions:
+        # integer vars to maintain
+        if reaction.id in integer_vars_to_maintain:
+            # keep these integer variables in the dual, with new transformed
+            # constraints
+            if (reaction.lower_bound not in [0, 1] or
+                    reaction.upper_bound not in [0, 1] or
+                    reaction.variable_kind != "integer"):
+                raise Exception("Reaction %s from integer_vars_to_maintain is "
+                                "not a Boolean integer variable" % reaction.id)
+            integer_var = Reaction(reaction.id)
+            integer_var.upper_bound = reaction.upper_bound
+            integer_var.lower_bound = reaction.lower_bound
+            integer_var.variable_kind = reaction.variable_kind
+            integer_var.objective_coefficient = 0
+            # constraints
+            dual.add_reaction(integer_var)
+            integer_vars_added.append((integer_var, reaction))
+
+        # other vars
+        else:
+            # other variables become constraints, (1) to (6)
+            constr = Metabolite("%s__dual_constrained_by_c" %
+                                reaction.id)  # (6)
+            constr._constraint_sense = "L"
+            constr._bound = - reaction.objective_coefficient
+            for met, coeff in iteritems(reaction.metabolites):
+                dual_var = dual_var_for_met[met.id]
+                dual_var.add_metabolites({constr: - coeff})
+
+            # upper bound constraints -> variables (3) to (5) and (6)
+            var_bound = Reaction("%s__dual_for_upper_bound_constraint" %
+                                 reaction.id)  # dual_m+j
+            var_bound.objective_coefficient = reaction.upper_bound  # (5)
+            # [dual_vars] >= 0
+            var_bound.lower_bound = 0
+            var_bound.upper_bound = dual_maximum
+            # add bound dual variables to dual constraints
+            var_bound.add_metabolites({constr: -1})  # (6)
+            dual.add_reaction(var_bound)
+
+    # add auxiliary variables
+    for integer_var, original_reaction in integer_vars_added:
+        for metabolite, coeff in iteritems(original_reaction.metabolites):
+            dual_var = dual_var_for_met[metabolite.id]
+            # create an auxiliary variable
+            aux_var = Reaction("%s__auxiliary__%s" % (integer_var.id,
+                                                      dual_var.id))
+            aux_var.lower_bound = 0
+            aux_var.upper_bound = dual_maximum
+            aux_var.objective_coefficient = - coeff
+            dual.add_reaction(aux_var)
+
+            # add linearization constraints
+            # (7)     auxiliary_var_i_j - dual_maximum * decision_var_j    <= 0
+            le_decision_constr = Metabolite("%s__le_decision" % aux_var.id)
+            le_decision_constr._constraint_sense = "L"
+            le_decision_constr._bound = 0
+            aux_var.add_metabolites({le_decision_constr: 1})
+            integer_var.add_metabolites({le_decision_constr: - dual_maximum})
+
+            # (8)     auxiliary_var_i_j - dual_i                           <= 0
+            le_dual_constr = Metabolite("%s__le_dual" % aux_var.id)
+            le_dual_constr._constraint_sense = "L"
+            le_dual_constr._bound = 0
+            aux_var.add_metabolites({le_dual_constr: 1})
+            dual_var.add_metabolites({le_dual_constr: -1})
+
+            # (9)   - auxiliary_var_i_j + dual_i +
+            #         dual_maximum * decision_var_j <= dual_maximum
+            g_constr = Metabolite("%s__g_dual" % aux_var.id)
+            g_constr._constraint_sense = "L"
+            g_constr._bound = dual_maximum
+            aux_var.add_metabolites({g_constr: -1})
+            dual_var.add_metabolites({g_constr: 1})
+            integer_var.add_metabolites({g_constr: dual_maximum})
+
+    return dual
diff --git a/cobra/flux_analysis/phenotype_phase_plane.py b/cobra/flux_analysis/phenotype_phase_plane.py
index 27ca550..aeb937d 100644
--- a/cobra/flux_analysis/phenotype_phase_plane.py
+++ b/cobra/flux_analysis/phenotype_phase_plane.py
@@ -1,5 +1,5 @@
 from numpy import linspace, zeros, array, meshgrid, abs, empty, arange, \
-    int32, unravel_index
+    int32, unravel_index, dtype
 from multiprocessing import Pool
 
 from ..solvers import solver_dict, get_solver_name
@@ -13,9 +13,12 @@ except ImportError:
     axes3d = None
 mlab = None  # mayavi may crash python
 try:  # for prettier colors
-    from brewer2mpl import get_map
-except:
-    get_map = None
+    from palettable.colorbrewer import get_map
+except ImportError:
+    try:
+        from brewer2mpl import get_map
+    except ImportError:
+        get_map = None
 
 
 class phenotypePhasePlaneData:
@@ -52,12 +55,12 @@ class phenotypePhasePlaneData:
     def plot_matplotlib(self, theme="Paired", scale_grid=False):
         """Use matplotlib to plot a phenotype phase plane in 3D.
 
-        theme: color theme to use (requires brewer2mpl)
+        theme: color theme to use (requires palettable)
 
         returns: maptlotlib 3d subplot object"""
         if pyplot is None:
             raise ImportError("Error importing matplotlib 3D plotting")
-        colors = empty(self.growth_rates.shape, dtype="|S7")
+        colors = empty(self.growth_rates.shape, dtype=dtype((str, 7)))
         n_segments = self.segments.max()
         # pick colors
         if get_map is None:
diff --git a/cobra/flux_analysis/reaction.py b/cobra/flux_analysis/reaction.py
index a81b7d4..c7e92d6 100644
--- a/cobra/flux_analysis/reaction.py
+++ b/cobra/flux_analysis/reaction.py
@@ -1,6 +1,7 @@
 #cobra.flux_analysis.reaction.py
 #functions for analyzing / creating objective functions
 from ..core.Reaction import Reaction
+from six import iteritems
 
 def assess(model, reaction, flux_coefficient_cutoff=0.001):
     """Assesses the capacity of the model to produce the precursors for the reaction
@@ -81,7 +82,7 @@ def assess_precursors(model, reaction,  flux_coefficient_cutoff=0.001):
 
     #Otherwise assess the ability of the model to produce each precursor individually.
     #Now assess the ability of the model to produce each reactant for a reaction
-    for sink_reaction, (component, coefficient) in sink_reactions.iteritems():
+    for sink_reaction, (component, coefficient) in iteritems(sink_reactions):
         model.optimize(new_objective=sink_reaction) #Calculate the maximum amount of the
         #metabolite that can be produced.
         if flux_coefficient_cutoff > model.solution.f:
@@ -142,7 +143,7 @@ def assess_products(model, reaction, flux_coefficient_cutoff=0.001):
         return(True)
 
     #Now assess the ability of the model to produce each reactant for a reaction
-    for source_reaction, (component, coefficient) in source_reactions.iteritems():
+    for source_reaction, (component, coefficient) in iteritems(source_reactions):
         model.optimize(new_objective=source_reaction) #Calculate the maximum amount of the
         #metabolite that can be produced.
         if flux_coefficient_cutoff > model.solution.f:
diff --git a/cobra/io/json.py b/cobra/io/json.py
index 9bdc1ca..f04a9b5 100644
--- a/cobra/io/json.py
+++ b/cobra/io/json.py
@@ -63,6 +63,8 @@ def _fix_type(value):
     # handle legacy Formula type
     if value.__class__.__name__ == "Formula":
         return str(value)
+    if value is None:
+        return ''
     return value
 
 
@@ -110,7 +112,7 @@ def _update_optional(cobra_object, new_dict, optional_attribute_dict):
     """update new_dict with optional attributes from cobra_object"""
     for key, default_value in iteritems(optional_attribute_dict):
         value = getattr(cobra_object, key)
-        if value != default_value:
+        if value is not None and value != default_value:
             new_dict[key] = _fix_type(value)
 
 
diff --git a/cobra/io/mat.py b/cobra/io/mat.py
index a7fb7f6..a7443de 100644
--- a/cobra/io/mat.py
+++ b/cobra/io/mat.py
@@ -102,13 +102,18 @@ def create_mat_dict(model):
     mat["mets"] = _cell(mets.list_attr("id"))
     mat["metNames"] = _cell(mets.list_attr("name"))
     mat["metFormulas"] = _cell([str(m.formula) for m in mets])
+    try:
+        mat["metCharge"] = array(mets.list_attr("charge")) * 1.
+    except TypeError:
+        # can't have any None entries for charge, or this will fail
+        pass
     mat["genes"] = _cell(model.genes.list_attr("id"))
     mat["grRules"] = _cell(rxns.list_attr("gene_reaction_rule"))
     mat["rxns"] = _cell(rxns.list_attr("id"))
     mat["rxnNames"] = _cell(rxns.list_attr("name"))
     mat["subSystems"] = _cell(rxns.list_attr("subsystem"))
     mat["csense"] = "".join(model._constraint_sense)
-    mat["S"] = model.S
+    mat["S"] = model.S if model.S is not None else [[]]
     # multiply by 1 to convert to float, working around scipy bug
     # https://github.com/scipy/scipy/issues/4537
     mat["lb"] = array(rxns.list_attr("lower_bound")) * 1.
@@ -155,6 +160,13 @@ def from_mat_struct(mat_struct, model_id=None):
             new_metabolite.formula = str(m["metFormulas"][0][0][i][0][0])
         except (IndexError, ValueError):
             pass
+        try:
+            new_metabolite.charge = float(m["metCharge"][0, 0][i][0])
+            int_charge = int(new_metabolite.charge)
+            if new_metabolite.charge == int_charge:
+                new_metabolite.charge = int_charge
+        except (IndexError, ValueError):
+            pass
         model.add_metabolites([new_metabolite])
     new_reactions = []
     for i, name in enumerate(m["rxns"][0, 0]):
diff --git a/cobra/io/sbml.py b/cobra/io/sbml.py
index 9002041..eb5ecde 100644
--- a/cobra/io/sbml.py
+++ b/cobra/io/sbml.py
@@ -8,6 +8,8 @@ from warnings import warn
 import re
 from math import isnan, isinf
 
+from six import iteritems
+
 #
 if __name == 'java':
     from org.sbml.jsbml import SBMLDocument, SpeciesReference, KineticLaw, Parameter
@@ -135,7 +137,8 @@ def create_cobra_model_from_sbml_file(sbml_filename, old_sbml=False, legacy_meta
         tmp_metabolite.name = sbml_metabolite.getName()
         tmp_formula = ''
         tmp_metabolite.notes = parse_legacy_sbml_notes(sbml_metabolite.getNotesString())
-        tmp_metabolite.charge = sbml_metabolite.getCharge()
+        if sbml_metabolite.isSetCharge():
+            tmp_metabolite.charge = sbml_metabolite.getCharge()
         if "CHARGE" in tmp_metabolite.notes:
             note_charge = tmp_metabolite.notes["CHARGE"][0]
             try:
@@ -145,8 +148,7 @@ def create_cobra_model_from_sbml_file(sbml_filename, old_sbml=False, legacy_meta
             except:
                 warn("charge of %s is not a number (%s)" % (tmp_metabolite.id, str(note_charge)))
             else:
-                if tmp_metabolite.charge == 0 or tmp_metabolite.charge == note_charge:  # get_charge() when unspecified is 0
-                    tmp_metabolite.charge = note_charge
+                if tmp_metabolite.charge is None or tmp_metabolite.charge == note_charge:
                     tmp_metabolite.notes.pop("CHARGE")
                 else:  # tmp_metabolite.charge != note_charge
                     msg = "different charges specified for %s (%d and %d)"
@@ -289,8 +291,9 @@ def create_cobra_model_from_sbml_file(sbml_filename, old_sbml=False, legacy_meta
                         gene_id_to_object[tmp_locus_id].name = tmp_row_dict['ABBREVIATION']
 
         if 'SUBSYSTEM' in reaction_note_dict:
-            reaction.subsystem = reaction_note_dict['SUBSYSTEM'][0]   
+            reaction.subsystem = reaction_note_dict.pop('SUBSYSTEM')[0]
 
+        reaction.notes = reaction_note_dict
 
 
     #Now, add all of the reactions to the model.
@@ -324,7 +327,7 @@ def parse_legacy_sbml_notes(note_string, note_delimiter = ':'):
             note_field = the_note[:note_delimiter_index].lstrip(' ').rstrip(' ').replace('_',' ').upper()
             note_value = the_note[(note_delimiter_index+1):].lstrip(' ').rstrip(' ')
             if note_field in note_dict:
-                note_dict[note_field ].append(note_value)
+                note_dict[note_field].append(note_value)
             else:
                 note_dict[note_field] = [note_value]
         note_string = note_string[(note_end+len(end_tag)): ]
@@ -334,7 +337,6 @@ def parse_legacy_sbml_notes(note_string, note_delimiter = ':'):
         
     return(note_dict)
 
-
 def write_cobra_model_to_sbml_file(cobra_model, sbml_filename,
                                    sbml_level=2, sbml_version=1,
                                    print_time=False,
@@ -358,6 +360,22 @@ def write_cobra_model_to_sbml_file(cobra_model, sbml_filename,
     Level 2 Version 4
 
     """
+
+    sbml_doc = get_libsbml_document(cobra_model, 
+                                   sbml_level=sbml_level, sbml_version=sbml_version,
+                                   print_time=print_time,
+                                   use_fbc_package=use_fbc_package)
+
+    writeSBML(sbml_doc, sbml_filename)
+
+def get_libsbml_document(cobra_model,
+                                   sbml_level=2, sbml_version=1,
+                                   print_time=False,
+                                   use_fbc_package=True):
+
+    """ Return a libsbml document object for writing to a file. This function
+    is used by write_cobra_model_to_sbml_file(). """
+
     note_start_tag, note_end_tag = '<p>', '</p>'
     if sbml_level > 2 or (sbml_level == 2 and sbml_version == 4):
         note_start_tag, note_end_tag = '<html:p>', '</html:p>'
@@ -475,13 +493,29 @@ def write_cobra_model_to_sbml_file(cobra_model, sbml_filename,
                 sbml_parameter.setValue(v)
             sbml_law.addParameter(sbml_parameter)
         sbml_reaction.setKineticLaw(sbml_law)
-        sbml_reaction.setNotes('<html xmlns="http://www.w3.org/1999/xhtml">%sGENE_ASSOCIATION: %s%s%sSUBSYSTEM: %s%s</html>'%(note_start_tag,
-                                                                 the_reaction.gene_reaction_rule,
-                                                                 note_end_tag,
-                                                                 note_start_tag,
-                                                                 the_reaction.subsystem,
-                                                                 note_end_tag))
 
+        #Checks if GPR and Subsystem annotations are present in the notes section and if they are the same as those in
+        #the reaction's gene_reaction_rule/ subsystem attribute
+        #If they are not identical, they are set to be identical
+        note_dict = the_reaction.notes.copy()
+        if the_reaction.gene_reaction_rule:
+            note_dict['GENE_ASSOCIATION'] = [str(the_reaction.gene_reaction_rule)]
+        if the_reaction.subsystem:
+            note_dict['SUBSYSTEM'] = [str(the_reaction.subsystem)]
+
+        #In a cobrapy model the notes section is stored as a dictionary. The following section turns the key-value-pairs
+        #of the dictionary into a string and replaces recurring symbols so that the string has the required syntax for
+        #an SBML doc.
+        note_str = str(list(iteritems(note_dict)))
+        note_start_tag, note_end_tag, note_delimiter = '<p>', '</p>', ':'
+        note_str = note_str.replace('(\'',note_start_tag)
+        note_str = note_str.replace('\']),',note_end_tag)
+        note_str = note_str.replace('\',',note_delimiter)
+        note_str = note_str.replace('\']','')
+        note_str = note_str.replace('[\'','')
+        note_str = note_str.replace('[','<html xmlns="http://www.w3.org/1999/xhtml">')
+        note_str = note_str.replace(')]',note_end_tag+'</html>')
+        sbml_reaction.setNotes(note_str)
 
     if use_fbc_package:
         try:
@@ -508,7 +542,7 @@ def write_cobra_model_to_sbml_file(cobra_model, sbml_filename,
                     error_string += "You've got libsbml %s installed.   You need 5.8.0 or later with the fbc package"
 
             raise(Exception(error_string))
-    writeSBML(sbml_doc, sbml_filename)
+    return sbml_doc
 
 def add_sbml_species(sbml_model, cobra_metabolite, note_start_tag,
                      note_end_tag, boundary_metabolite=False):
@@ -551,13 +585,17 @@ def add_sbml_species(sbml_model, cobra_metabolite, note_start_tag,
             return cobra_metabolite
     if cobra_metabolite.charge is not None:
         sbml_species.setCharge(cobra_metabolite.charge)
-    if hasattr(cobra_metabolite.formula, 'id') or \
-       hasattr(cobra_metabolite.notes, 'items'):
+    #Deal with cases where the formula in the model is not an object but as a string
+    if cobra_metabolite.formula or cobra_metabolite.notes:
         tmp_note =  '<html xmlns="http://www.w3.org/1999/xhtml">'
         if hasattr(cobra_metabolite.formula, 'id'):
             tmp_note += '%sFORMULA: %s%s'%(note_start_tag,
                                               cobra_metabolite.formula.id,
                                               note_end_tag)
+        else:
+            tmp_note += '%sFORMULA: %s%s'%(note_start_tag,
+                                  cobra_metabolite.formula,
+                                  note_end_tag)
         if hasattr(cobra_metabolite.notes, 'items'):
             for the_id_type, the_id in cobra_metabolite.notes.items():
                 if the_id_type.lower() == 'charge':
diff --git a/cobra/io/sbml3.py b/cobra/io/sbml3.py
index f699928..a9faea7 100644
--- a/cobra/io/sbml3.py
+++ b/cobra/io/sbml3.py
@@ -1,7 +1,6 @@
 from collections import defaultdict
 from warnings import warn, catch_warnings, simplefilter
 from decimal import Decimal
-from math import isinf, isnan
 from ast import parse as ast_parse, Name, Or, And, BoolOp
 from gzip import GzipFile
 from bz2 import BZ2File
@@ -12,6 +11,9 @@ from six import iteritems, string_types
 
 from .. import Metabolite, Reaction, Gene, Model
 from ..core.Gene import parse_gpr
+from ..manipulation.modify import _renames
+from ..manipulation.validate import check_reaction_bounds, \
+    check_metabolite_compartment_formula
 
 try:
     from lxml.etree import parse, Element, SubElement, \
@@ -44,27 +46,6 @@ except:
         pass
 
 
-_renames = (
-    (".", "_DOT_"),
-    ("(", "_LPAREN_"),
-    (")", "_RPAREN_"),
-    ("-", "__"),
-    ("[", "_LSQBKT"),
-    ("]", "_RSQBKT"),
-    (",", "_COMMA_"),
-    (":", "_COLON_"),
-    (">", "_GT_"),
-    ("<", "_LT"),
-    ("/", "_FLASH"),
-    ("\\", "_BSLASH"),
-    ("+", "_PLUS_"),
-    ("=", "_EQ_"),
-    (" ", "_SPACE_"),
-    ("'", "_SQUOT_"),
-    ('"', "_DQUOT_"),
-)
-
-
 # deal with namespaces
 namespaces = {"fbc": "http://www.sbml.org/sbml/level3/version1/fbc/version2",
               "sbml": "http://www.sbml.org/sbml/level3/version1/core",
@@ -129,8 +110,11 @@ class CobraSBMLError(Exception):
 def get_attrib(tag, attribute, type=lambda x: x, require=False):
     value = tag.get(ns(attribute))
     if require and value is None:
-        raise CobraSBMLError("required attribute '%s' not found in tag '%s'" %
-                             (attribute, tag.tag))
+        msg = "required attribute '%s' not found in tag '%s'" % \
+                             (attribute, tag.tag)
+        if tag.get("id") is not None:
+            msg += "with id '%s'" % tag.get("id")
+        raise CobraSBMLError(msg)
     return type(value) if value is not None else None
 
 
@@ -154,7 +138,7 @@ def parse_stream(filename):
         else:
             return parse(filename)
     except ParseError as e:
-        raise CobraSBMLError("Malformed XML file: " + e.message)
+        raise CobraSBMLError("Malformed XML file: " + str(e))
 
 
 # string utility functions
@@ -208,7 +192,7 @@ def annotate_cobra_from_sbml(cobra_element, sbml_element):
             warn("%s does not start with http://identifiers.org/" % uri)
             continue
         try:
-            provider, identifier = uri[23:].split("/")
+            provider, identifier = uri[23:].split("/", 1)
         except ValueError:
             warn("%s does not conform to http://identifiers.org/provider/id"
                  % uri)
@@ -337,6 +321,8 @@ def parse_xml_into_model(xml, number=float):
         object_stoichiometry = {}
         for met_id in stoichiometry:
             if met_id in boundary_metabolites:
+                warn("Boundary metabolite '%s' used in reaction '%s'" %
+                     (met_id, reaction.id))
                 continue
             try:
                 metabolite = model.metabolites.get_by_id(met_id)
@@ -368,7 +354,11 @@ def parse_xml_into_model(xml, number=float):
     obj_query = OBJECTIVES_XPATH % target_objective
     for sbml_objective in obj_list.findall(obj_query):
         rxn_id = clip(get_attrib(sbml_objective, "fbc:reaction"), "R_")
-        model.reactions.get_by_id(rxn_id).objective_coefficient = \
+        try:
+            objective_reaction = model.reactions.get_by_id(rxn_id)
+        except KeyError as e:
+            raise CobraSBMLError("Objective reaction '%s' not found" % rxn_id)
+        objective_reaction.objective_coefficient = \
             get_attrib(sbml_objective, "fbc:coefficient", type=number)
 
     return model
@@ -413,8 +403,13 @@ def model_to_xml(cobra_model, units=True):
     if units:
         param_attr["units"] = "mmol_per_gDW_per_hr"
     # the most common bounds are the minimum, maxmium, and 0
-    min_value = min(cobra_model.reactions.list_attr("lower_bound"))
-    max_value = max(cobra_model.reactions.list_attr("upper_bound"))
+    if len(cobra_model.reactions) > 0:
+        min_value = min(cobra_model.reactions.list_attr("lower_bound"))
+        max_value = max(cobra_model.reactions.list_attr("upper_bound"))
+    else:
+        min_value = -1000
+        max_value = 1000
+
     SubElement(parameter_list, "parameter", value=strnum(min_value),
                id="cobra_default_lb", sboTerm="SBO:0000626", **param_attr)
     SubElement(parameter_list, "parameter", value=strnum(max_value),
@@ -562,7 +557,7 @@ id_required = {ns(i) for i in ("sbml:model", "sbml:reaction:", "sbml:species",
 invalid_id_detector = re.compile("|".join(re.escape(i[0]) for i in _renames))
 
 
-def validate_sbml_model(filename):
+def validate_sbml_model(filename, check_model=True):
     """returns the model along with a list of errors"""
     xmlfile = parse_stream(filename)
     xml = xmlfile.getroot()
@@ -576,6 +571,15 @@ def validate_sbml_model(filename):
     def err(err_msg):
         sbml_errors.append(err_msg)
 
+    # make sure there is exactly one model
+    xml_models = xml.findall(ns("sbml:model"))
+    if len(xml_models) > 1:
+        err("More than 1 SBML model detected in file")
+    elif len(xml_models) == 0:
+        err("No SBML model detected in file")
+    else:
+        xml_model = xml_models[0]
+
     # make sure all sbml id's are valid
     all_ids = set()
     for element in xmlfile.iter():
@@ -602,7 +606,7 @@ def validate_sbml_model(filename):
                     (str_id[e.start:e.end], id_repr))
             if invalid_id_detector.search(str_id):
                 bad_chars = "".join(invalid_id_detector.findall(str_id))
-                err("invalid character%s %s found in %s" %
+                err("invalid character%s '%s' found in %s" %
                     ("s" if len(bad_chars) > 1 else "", bad_chars, id_repr))
             if not str_id[0].isalpha():
                 err("%s does not start with alphabet character" % id_repr)
@@ -611,31 +615,30 @@ def validate_sbml_model(filename):
     for element in xml.findall(".//*[@sboTerm]"):
         sbo_term = element.get("sboTerm")
         if not sbo_term.startswith("SBO:"):
-            warn("sboTerm '%s' does not begin with 'SBO:'" % sbo_term)
+            err("sboTerm '%s' does not begin with 'SBO:'" % sbo_term)
 
     # ensure can be made into model
+    # all warnings generated while loading will be logged as errors
     with catch_warnings(record=True) as warning_list:
         simplefilter("always")
-        model = parse_xml_into_model(xml)
-    sbml_errors.extend(i.message.message for i in warning_list)
-
-    # ensure exactly one objective
-    if len(model.objective) == 0:
-        sbml_errors.append("no objective reaction identified")
-    elif len(model.objective) > 1:
-        sbml_errors.append("only one reaction should be the objective")
-
-    # make sure there are no infinite bounds
-    if any(isinf(i) or isnan(i)
-           for i in model.reactions.list_attr("lower_bound")):
-        sbml_errors.append("infinite or NaN value detected in lower bounds")
-    if any(isinf(i) or isnan(i)
-           for i in model.reactions.list_attr("upper_bound")):
-        sbml_errors.append("infinite or NaN value detected in upper bounds")
-    for reaction in model.reactions:
-        if reaction.lower_bound > reaction.upper_bound:
-            sbml_errors.append("reaction '%s' has lower bound > upper bound" %
-                               reaction.id)
+        try:
+            model = parse_xml_into_model(xml)
+        except CobraSBMLError as e:
+            err(str(e))
+            return (None, sbml_errors)
+    sbml_errors.extend(str(i.message) for i in warning_list)
+
+    # check genes
+    xml_genes = {
+        get_attrib(i, "fbc:id").replace(SBML_DOT, ".")
+        for i in xml_model.iterfind(GENES_XPATH)}
+    for gene in model.genes:
+        if "G_" + gene.id not in xml_genes and gene.id not in xml_genes:
+            err("No gene specfied with id 'G_%s'" % gene.id)
+
+    if check_model:
+        sbml_errors.extend(check_reaction_bounds(model))
+        sbml_errors.extend(check_metabolite_compartment_formula(model))
 
     return model, sbml_errors
 
@@ -648,9 +651,10 @@ def write_sbml_model(cobra_model, filename, use_fbc_package=True, **kwargs):
         return
     # create xml
     xml = model_to_xml(cobra_model, **kwargs)
-    write_args = {"encoding": "UTF-8"}
+    write_args = {"encoding": "UTF-8", "xml_declaration": True}
     if _with_lxml:
         write_args["pretty_print"] = True
+        write_args["pretty_print"] = True
     else:
         indent_xml(xml)
     # write xml to file
diff --git a/cobra/manipulation/__init__.py b/cobra/manipulation/__init__.py
index 827c4f1..6bf72f7 100644
--- a/cobra/manipulation/__init__.py
+++ b/cobra/manipulation/__init__.py
@@ -1,3 +1,8 @@
-from .delete import delete_model_genes, undelete_model_genes, remove_genes
+from .delete import delete_model_genes, undelete_model_genes, remove_genes, \
+    find_gene_knockout_reactions
 from .modify import initialize_growth_medium, convert_to_irreversible, \
-    revert_to_reversible, escape_ID
+    revert_to_reversible, escape_ID, canonical_form, \
+    get_compiled_gene_reaction_rules
+from .annotate import add_SBO
+from .validate import check_mass_balance, check_reaction_bounds, \
+    check_metabolite_compartment_formula
diff --git a/cobra/manipulation/annotate.py b/cobra/manipulation/annotate.py
new file mode 100644
index 0000000..715107a
--- /dev/null
+++ b/cobra/manipulation/annotate.py
@@ -0,0 +1,24 @@
+from six import iteritems
+
+
+def add_SBO(model):
+    """adds SBO terms for demands and exchanges
+
+    This works for models which follow the standard convention for
+    constructing and naming these reactions.
+
+    The reaction should only contain the single metabolite being exchanged,
+    and the id should be EX_metid or DM_metid
+    """
+    for r in model.reactions:
+        # don't annotate already annotated reactions
+        if r.annotation.get("SBO"):
+            continue
+        # only doing exchanges
+        if len(r.metabolites) != 1:
+            continue
+        met_id = list(r.metabolites)[0].id
+        if r.id.startswith("EX_") and r.id == "EX_" + met_id:
+            r.annotation["SBO"] = "SBO:0000627"
+        elif r.id.startswith("DM_") and r.id == "DM_" + met_id:
+            r.annotation["SBO"] = "SBO:0000628"
diff --git a/cobra/manipulation/modify.py b/cobra/manipulation/modify.py
index d24bf6d..75a2e9d 100644
--- a/cobra/manipulation/modify.py
+++ b/cobra/manipulation/modify.py
@@ -5,10 +5,30 @@ from ast import NodeTransformer
 
 from six import iteritems
 
-from .. import Reaction
+from .. import Reaction, Metabolite
 from .delete import get_compiled_gene_reaction_rules
 from ..core.Gene import ast2str
-from ..io.sbml3 import _renames
+
+
+_renames = (
+    (".", "_DOT_"),
+    ("(", "_LPAREN_"),
+    (")", "_RPAREN_"),
+    ("-", "__"),
+    ("[", "_LSQBKT"),
+    ("]", "_RSQBKT"),
+    (",", "_COMMA_"),
+    (":", "_COLON_"),
+    (">", "_GT_"),
+    ("<", "_LT"),
+    ("/", "_FLASH"),
+    ("\\", "_BSLASH"),
+    ("+", "_PLUS_"),
+    ("=", "_EQ_"),
+    (" ", "_SPACE_"),
+    ("'", "_SQUOT_"),
+    ('"', "_DQUOT_"),
+)
 
 
 def _escape_str_id(id_str):
@@ -160,21 +180,23 @@ def convert_to_irreversible(cobra_model):
         # will be constrained to 0) will be left in the model.
         if reaction.lower_bound < 0:
             reverse_reaction = Reaction(reaction.id + "_reverse")
-            reverse_reaction.lower_bound = 0
-            reverse_reaction.upper_bound = reaction.lower_bound * -1
+            reverse_reaction.lower_bound = max(0, -reaction.upper_bound)
+            reverse_reaction.upper_bound = -reaction.lower_bound
             reverse_reaction.objective_coefficient = \
                 reaction.objective_coefficient * -1
-            reaction.lower_bound = 0
+            reaction.lower_bound = max(0, reaction.lower_bound)
+            reaction.upper_bound = max(0, reaction.upper_bound)
             # Make the directions aware of each other
             reaction.notes["reflection"] = reverse_reaction.id
             reverse_reaction.notes["reflection"] = reaction.id
-            reaction_dict = dict([(k, v*-1)
-                                  for k, v in reaction._metabolites.items()])
+            reaction_dict = {k: v * -1
+                             for k, v in iteritems(reaction._metabolites)}
             reverse_reaction.add_metabolites(reaction_dict)
             reverse_reaction._model = reaction._model
             reverse_reaction._genes = reaction._genes
             for gene in reaction._genes:
                 gene._reaction.add(reverse_reaction)
+            reverse_reaction.subsystem = reaction.subsystem
             reverse_reaction._gene_reaction_rule = reaction._gene_reaction_rule
             reactions_to_add.append(reverse_reaction)
     cobra_model.add_reactions(reactions_to_add)
@@ -205,6 +227,8 @@ def revert_to_reversible(cobra_model, update_solution=True):
         forward_id = reverse.notes.pop("reflection")
         forward = cobra_model.reactions.get_by_id(forward_id)
         forward.lower_bound = -reverse.upper_bound
+        if forward.upper_bound == 0:
+            forward.upper_bound = -reverse.lower_bound
 
         # update the solution dict
         if update_solution:
@@ -223,3 +247,75 @@ def revert_to_reversible(cobra_model, update_solution=True):
     if update_solution:
         cobra_model.solution.x_dict = x_dict
         cobra_model.solution.x = [x_dict[r.id] for r in cobra_model.reactions]
+
+
+def canonical_form(model, objective_sense='maximize',
+                   already_irreversible=False, copy=True):
+    """Return a model (problem in canonical_form).
+
+    Converts a minimization problem to a maximization, makes all variables
+    positive by making reactions irreversible, and converts all constraints to
+    <= constraints.
+
+
+    model: class:`~cobra.core.Model`. The model/problem to convert.
+
+    objective_sense: str. The objective sense of the starting problem, either
+    'maximize' or 'minimize'. A minimization problems will be converted to a
+    maximization.
+
+    already_irreversible: bool. If the model is already irreversible, then pass
+    True.
+
+    copy: bool. Copy the model before making any modifications.
+
+    """
+    if copy:
+        model = model.copy()
+
+    if not already_irreversible:
+        convert_to_irreversible(model)
+
+    if objective_sense == "minimize":
+        # if converting min to max, reverse all the objective coefficients
+        for reaction in model.reactions:
+            reaction.objective_coefficient = - reaction.objective_coefficient
+    elif objective_sense != "maximize":
+        raise Exception("Invalid objective sense '%s'. "
+                        "Must be 'minimize' or 'maximize'." % objective_sense)
+
+    # convert G and E constraints to L constraints
+    for metabolite in model.metabolites:
+        if metabolite._constraint_sense == "G":
+            metabolite._constraint_sense = "L"
+            metabolite._bound = - metabolite._bound
+            for reaction in metabolite.reactions:
+                coeff = reaction.get_coefficient(metabolite)
+                # reverse the coefficient
+                reaction.add_metabolites({metabolite: -2 * coeff})
+        elif metabolite._constraint_sense == "E":
+            # change existing constraint to L
+            metabolite._constraint_sense = "L"
+            # add new constraint
+            new_constr = Metabolite("%s__GE_constraint" % metabolite.id)
+            new_constr._constraint_sense = "L"
+            new_constr._bound = - metabolite._bound
+            for reaction in metabolite.reactions:
+                coeff = reaction.get_coefficient(metabolite)
+                reaction.add_metabolites({new_constr: -coeff})
+
+    # convert lower bounds to LE constraints
+    for reaction in model.reactions:
+        if reaction.lower_bound < 0:
+            raise Exception("Bounds of irreversible reactions should be >= 0,"
+                            " for %s" % reaction.id)
+        elif reaction.lower_bound == 0:
+            continue
+        # new constraint for lower bound
+        lb_constr = Metabolite("%s__LB_constraint" % reaction.id)
+        lb_constr._constraint_sense = "L"
+        lb_constr._bound = - reaction.lower_bound
+        reaction.add_metabolites({lb_constr: -1})
+        reaction.lower_bound = 0
+
+    return model
diff --git a/cobra/manipulation/validate.py b/cobra/manipulation/validate.py
new file mode 100644
index 0000000..7d5aa72
--- /dev/null
+++ b/cobra/manipulation/validate.py
@@ -0,0 +1,54 @@
+from math import isinf, isnan
+
+NOT_MASS_BALANCED_TERMS = {"SBO:0000627",  # EXCHANGE
+                           "SBO:0000628",  # DEMAND
+                           "SBO:0000629",  # BIOMASS
+                           "SBO:0000631",  # PSEUDOREACTION
+                           "SBO:0000632",  # SINK
+                           }
+
+
+def check_mass_balance(model):
+    warnings = []
+    unbalanced = {}
+    for reaction in model.reactions:
+        if reaction.annotation.get("SBO") not in NOT_MASS_BALANCED_TERMS:
+            balance = reaction.check_mass_balance()
+            if balance:
+                unbalanced[reaction] = balance
+    return unbalanced
+
+
+def check_reaction_bounds(model):
+    errors = []
+    for reaction in model.reactions:
+        if reaction.lower_bound > reaction.upper_bound:
+            errors.append("Reaction '%s' has lower bound > upper bound" %
+                          reaction.id)
+        if isinf(reaction.lower_bound):
+            errors.append("Reaction '%s' has infinite lower_bound" %
+                          reaction.id)
+        elif isnan(reaction.lower_bound):
+            errors.append("Reaction '%s' has NaN for lower_bound" %
+                          reaction.id)
+        if isinf(reaction.upper_bound):
+            errors.append("Reaction '%s' has infinite upper_bound" %
+                          reaction.id)
+        elif isnan(reaction.upper_bound):
+            errors.append("Reaction '%s' has NaN for upper_bound" %
+                          reaction.id)
+    return errors
+
+
+def check_metabolite_compartment_formula(model):
+    errors = []
+    for met in model.metabolites:
+        if met.compartment is not None and \
+                met.compartment not in model.compartments:
+            errors.append("Metabolite '%s' compartment '%s' not found" %
+                          (met.id, met.compartment))
+        if met.formula is not None and len(met.formula) > 0:
+            if not met.formula.isalnum():
+                errors.append("Metabolite '%s' formula '%s' not alphanumeric" %
+                              (met.id, met.formula))
+    return errors
diff --git a/cobra/solvers/__init__.py b/cobra/solvers/__init__.py
index afd7d14..0d90690 100644
--- a/cobra/solvers/__init__.py
+++ b/cobra/solvers/__init__.py
@@ -77,9 +77,9 @@ def get_solver_name(mip=False, qp=False):
     if len(solver_dict) == 0:
         raise SolverNotFound("no solvers installed")
     # glpk only does lp, not qp. Gurobi and cplex are better at mip
-    mip_order = ["gurobi", "cplex", "coin", "cglpk", "glpk"]
-    lp_order = ["cglpk", "cplex",  "glpk", "gurobi", "coin"]
-    qp_order = ["gurobi", "cplex"]
+    mip_order = ["gurobi", "cplex", "mosek", "coin", "cglpk", "glpk"]
+    lp_order = ["cglpk", "cplex",  "gurobi", "mosek", "coin", "glpk"]
+    qp_order = ["gurobi", "cplex", "mosek"]
 
     if mip is False and qp is False:
         for solver_name in lp_order:
diff --git a/cobra/solvers/cglpk.pyx b/cobra/solvers/cglpk.pyx
index b375914..eb1a989 100644
--- a/cobra/solvers/cglpk.pyx
+++ b/cobra/solvers/cglpk.pyx
@@ -150,6 +150,7 @@ cdef class GLP:
         self.exact = False
         glp_term_hook(hook, NULL)
         self.parameters.msg_lev = GLP_MSG_OFF
+        self.integer_parameters.tol_int = 1e-9
 
     def __dealloc__(self):
         glp_delete_prob(self.glp)
@@ -299,7 +300,15 @@ cdef class GLP:
             check_error(glp_simplex(glp, &self.parameters))
         self.parameters.tm_lim = time_limit
         if self.exact:
-            check_error(glp_exact(glp, &self.parameters))
+            # sigh... it looks like the exact routine doesn't fully respect
+            # the verbosity parameter
+            if self.parameters.msg_lev == GLP_MSG_OFF:
+                glp_term_hook(silent_hook, NULL)
+            try:
+                check_error(glp_exact(glp, &self.parameters))
+            finally:
+                glp_term_hook(hook, NULL)
+
         if self.is_mip():
             self.integer_parameters.tm_lim = self.parameters.tm_lim
             self.integer_parameters.msg_lev = self.parameters.msg_lev
diff --git a/cobra/solvers/cplex_solver.py b/cobra/solvers/cplex_solver.py
index 3b1ea14..9072842 100644
--- a/cobra/solvers/cplex_solver.py
+++ b/cobra/solvers/cplex_solver.py
@@ -56,6 +56,7 @@ status_dict = {'MIP_infeasible': 'infeasible',
                'MIP_optimal_tolerance': 'optimal',
                'MIP_unbounded':  'unbounded',
                'infeasible': 'infeasible',
+               'integer infeasible': 'infeasible',
                'optimal': 'optimal',
                'optimal_tolerance': 'optimal',
                'unbounded': 'unbounded',
@@ -240,7 +241,7 @@ def set_quadratic_objective(lp, quadratic_objective):
         raise Exception('quadratic component must have method todok')
     # ensure the matrix is properly read in
     nnz = quadratic_objective.nnz
-    if lp.parameters.read.qpnonzeros < nnz:
+    if lp.parameters.read.qpnonzeros.get() < nnz:
         lp.parameters.read.qpnonzeros.set(nnz + 1)
     # Reset the quadratic coefficient if it exists
     if lp.objective.get_num_quadratic_nonzeros() > 0:
diff --git a/cobra/solvers/cplex_solver_java.py b/cobra/solvers/cplex_solver_java.py
index a5ad3f2..56cd03c 100644
--- a/cobra/solvers/cplex_solver_java.py
+++ b/cobra/solvers/cplex_solver_java.py
@@ -1,6 +1,7 @@
-#cobra.solvers.cplex_solver
+# PLEASE NOTE THAT JYTHON SUPPORT (and this jython-only-solver) is deprecated
 #Interface to ilog/cplex 12.4 python / jython interfaces
 #QPs are not yet supported under jython
+from __future__ import print_function
 from os import name as __name
 from copy import deepcopy
 from warnings import warn
@@ -11,6 +12,7 @@ from .parameters import status_dict, variable_kind_dict, \
 
 from ..core.Solution import Solution
 from time import time
+from six import iteritems
 solver_name = 'cplex'
 parameter_defaults = parameter_defaults[solver_name]
 sense_dict = eval(sense_dict[solver_name])
@@ -120,7 +122,7 @@ def create_problem(cobra_model,  **kwargs):
         lp.setWarning(None)
         lp.setOut(None)
     [set_parameter(lp, parameter_mappings[k], v)
-     for k, v in the_parameters.iteritems() if k in parameter_mappings]
+     for k, v in iteritems(the_parameters) if k in parameter_mappings]
     quadratic_component = the_parameters['quadratic_component']
     new_objective = the_parameters['new_objective']
     error_reporting = the_parameters['error_reporting']
@@ -240,7 +242,7 @@ def solve_problem(lp, **kwargs):
     #Update parameter settings if provided
     if kwargs:
         [set_parameter(lp, parameter_mappings[k], v)
-         for k, v in kwargs.iteritems() if k in parameter_mappings]
+         for k, v in iteritems(kwargs) if k in parameter_mappings]
     try:
         the_problem = kwargs['the_problem']
     except:
@@ -269,10 +271,8 @@ def solve(cobra_model, **kwargs):
     if kwargs:
         the_parameters.update(kwargs)
     #Update objectives if they are new.
-    if 'new_objective' in the_parameters and \
-           the_parameters['new_objective'] not in ['update problem', None]:
-       from ..flux_analysis.objective import update_objective
-       update_objective(cobra_model, the_parameters['new_objective'])
+    if 'new_objective' in the_parameters:
+        raise ValueError("new_objective option removed")
 
     if 'the_problem' in the_parameters:
         the_problem = the_parameters['the_problem']
@@ -315,7 +315,7 @@ def solve(cobra_model, **kwargs):
 
     the_solution = format_solution(lp, cobra_model)
     if status != 'optimal' and error_reporting:
-        print '%s failed: %s'%(solver_name, status)
+        print('{:s} failed: {:s}'.format(solver_name, status))
     cobra_model.solution = the_solution
     solution = {'the_problem': lp, 'the_solution': the_solution}
     return solution
diff --git a/cobra/solvers/glpk_solver.py b/cobra/solvers/glpk_solver.py
index eabfd6e..7115d36 100644
--- a/cobra/solvers/glpk_solver.py
+++ b/cobra/solvers/glpk_solver.py
@@ -2,8 +2,13 @@
 #This script provides wrappers for pyglpk 0.3
 from warnings import warn
 from copy import deepcopy
-from itertools import izip
+try:
+    # Import izip for python versions < 3.x
+    from itertools import izip as zip
+except ImportError:
+    pass
 
+from six import iteritems
 from glpk import LPX
 
 from ..core.Solution import Solution
@@ -79,7 +84,7 @@ def create_problem(cobra_model, **kwargs):
     lp.rows.add(len(cobra_model.metabolites))
     lp.cols.add(len(cobra_model.reactions))
 
-    for r, the_metabolite in izip(lp.rows, cobra_model.metabolites):
+    for r, the_metabolite in zip(lp.rows, cobra_model.metabolites):
         r.name = the_metabolite.id
         b = float(the_metabolite._bound)
         c = the_metabolite._constraint_sense
@@ -94,12 +99,12 @@ def create_problem(cobra_model, **kwargs):
 
     objective_coefficients = []
     linear_constraints = []
-    for c, the_reaction in izip(lp.cols, cobra_model.reactions):
+    for c, the_reaction in zip(lp.cols, cobra_model.reactions):
         c.name = the_reaction.id
         c.kind = variable_kind_dict[the_reaction.variable_kind]
         c.bounds = the_reaction.lower_bound, the_reaction.upper_bound
         objective_coefficients.append(float(the_reaction.objective_coefficient))
-        for metabolite, coefficient in the_reaction._metabolites.iteritems():
+        for metabolite, coefficient in iteritems(the_reaction._metabolites):
             metabolite_index = metabolite_to_index[metabolite]
             linear_constraints.append((metabolite_index, c.index, coefficient))
 
diff --git a/cobra/solvers/glpk_solver_java.py b/cobra/solvers/glpk_solver_java.py
index 0777b36..3fa9a19 100644
--- a/cobra/solvers/glpk_solver_java.py
+++ b/cobra/solvers/glpk_solver_java.py
@@ -1,5 +1,6 @@
-##cobra.solvers.glpk_solver
+# PLEASE NOTE THAT JYTHON SUPPORT (and this jython-only-solver) is deprecated
 #This script provides wrappers for libglpk-java 1.0.22 and pyglpk 0.3
+from __future__ import print_function
 from warnings import warn
 from copy import deepcopy
 ###solver specific parameters
@@ -9,6 +10,7 @@ from .parameters import status_dict, variable_kind_dict, \
 
 from ..core.Solution import Solution
 from time import time
+from six import iteritems
 solver_name = 'glpk'
 sense_dict = eval(sense_dict[solver_name])
 #Functions that are different for java implementation of a solver
@@ -83,18 +85,18 @@ class Problem():
             try:
                 setattr(self._simplex_parameters, parameter_name,
                         parameter_value)
-            except Exception, e1:
+            except Exception as e1:
                 try:
                     setattr(self._mip_parameters, parameter_name,
                             parameter_value)
-                except Exception, e2:
+                except Exception as e2:
                     if warning:
-                        print "Could not set simplex parameter " +\
-                              "%s: %s"%(parameter_name, repr(e1))
+                        print("Could not set simplex parameter " +\
+                              "{:s}: {:s}".format(parameter_name, repr(e1)))
                         
                         if self._mip_parameters is not None:
-                            print "Could not set mip parameter " +\
-                                  "%s: %s"%(parameter_name, repr(e2))
+                            print("Could not set mip parameter " +\
+                                  "{:s}: {:s}".format(parameter_name, repr(e2)))
     def get_objective_value(self):
         if self._mip:
             tmp_value = self._g.glp_mip_obj_val(self._lp)
@@ -190,7 +192,7 @@ class Problem():
 
         g.glp_set_obj_name(lp, "z")
         [g.glp_set_obj_coef(lp, k, v)
-          for k, v in objective_dict.iteritems()]
+          for k, v in iteritems(objective_dict)]
 
         
 
@@ -222,8 +224,8 @@ def format_solution(lp, cobra_model, **kwargs):
             y_dict = dict(zip(cobra_model.metabolites, y))
         
             objective_value = lp.objective_value
-        except Exception, e:
-            print repr(e)
+        except Exception as e:
+            print(repr(e))
             y = y_dict = x = x_dict = objective_value = None
             #print status
     else:
@@ -255,7 +257,7 @@ def create_problem(cobra_model,  **kwargs):
     lp = Problem()        # Create empty problem instance
     lp.create_problem(cobra_model)
     [set_parameter(lp, parameter_mappings[k], v)
-     for k, v in the_parameters.iteritems() if k in parameter_mappings]
+     for k, v in iteritems(the_parameters) if k in parameter_mappings]
     return(lp)
 
 def update_problem(lp, cobra_model, **kwargs):
@@ -288,7 +290,7 @@ def solve_problem(lp, **kwargs):
     #Update parameter settings if provided
     if kwargs:
         [set_parameter(lp, parameter_mappings[k], v)
-         for k, v in kwargs.iteritems() if k in parameter_mappings]
+         for k, v in iteritems(kwargs) if k in parameter_mappings]
     try:
         print_solver_time = kwargs['print_solver_time']
         start_time = time()
@@ -298,7 +300,7 @@ def solve_problem(lp, **kwargs):
     lp.solve()
     status = get_status(lp)
     if print_solver_time:
-        print 'optimize time: %f'%(time() - start_time)
+        print('optimize time: {:f}'.format(time() - start_time))
     return status
 
     
@@ -359,7 +361,7 @@ def solve(cobra_model, **kwargs):
     
     the_solution = format_solution(lp, cobra_model)
     if status != 'optimal' and error_reporting:
-        print '%s failed: %s'%(solver_name, status)
+        print('{:s} failed: {:s}'.format(solver_name, status))
     cobra_model.solution = the_solution
     solution = {'the_problem': lp, 'the_solution': the_solution}
     return solution
diff --git a/cobra/solvers/gurobi_solver.py b/cobra/solvers/gurobi_solver.py
index 57ad4ee..b2953b8 100644
--- a/cobra/solvers/gurobi_solver.py
+++ b/cobra/solvers/gurobi_solver.py
@@ -1,13 +1,37 @@
 # Interface to gurobipy
 
 from warnings import warn
-from itertools import izip
+from multiprocessing import Process
+import platform
+try:
+    # Import izip for python versions < 3.x
+    from itertools import izip as zip
+except ImportError:
+    pass
+
+
+def test_import():
+    """Sometimes trying to import gurobipy can segfault. To prevent this from
+    crashing everything, ensure it can be imported in a separate process."""
+    try:
+        import gurobipy
+    except ImportError:
+        pass
+
+if platform.system() != "Windows":
+    # https://github.com/opencobra/cobrapy/issues/207
+    p = Process(target=test_import)
+    p.start()
+    p.join()
+    if p.exitcode != 0:
+        raise RuntimeError("importing gurobi causes a crash (exitcode %d)" %
+                           p.exitcode)
 
 from gurobipy import Model, LinExpr, GRB, QuadExpr
 
 from ..core.Solution import Solution
 
-from six import string_types
+from six import string_types, iteritems
 
 try:
     from sympy import Basic, Number
@@ -82,13 +106,13 @@ def format_solution(lp, cobra_model, **kwargs):
     else:
         objective_value = lp.ObjVal
         x = [v.X for v in lp.getVars()]
-        x_dict = {r.id: value for r, value in izip(cobra_model.reactions, x)}
+        x_dict = {r.id: value for r, value in zip(cobra_model.reactions, x)}
         if lp.isMIP:
             y = y_dict = None  # MIP's don't have duals
         else:
             y = [c.Pi for c in lp.getConstrs()]
             y_dict = {m.id: value for m, value
-                      in izip(cobra_model.metabolites, y)}
+                      in zip(cobra_model.metabolites, y)}
         the_solution = Solution(objective_value, x=x, x_dict=x_dict, y=y,
                                 y_dict=y_dict, status=status)
     return(the_solution)
@@ -164,7 +188,7 @@ def create_problem(cobra_model, quadratic_component=None, **kwargs):
         the_parameters = parameter_defaults.copy()
         the_parameters.update(kwargs)
 
-    for k, v in the_parameters.iteritems():
+    for k, v in iteritems(the_parameters):
         set_parameter(lp, k, v)
 
 
@@ -228,7 +252,7 @@ def solve_problem(lp, **kwargs):
 
     """
     #Update parameter settings if provided
-    for k, v in kwargs.iteritems():
+    for k, v in iteritems(kwargs):
         set_parameter(lp, k, v)
 
     lp.update()
diff --git a/cobra/solvers/gurobi_solver_java.py b/cobra/solvers/gurobi_solver_java.py
index 4cea19f..6893f1f 100644
--- a/cobra/solvers/gurobi_solver_java.py
+++ b/cobra/solvers/gurobi_solver_java.py
@@ -1,9 +1,11 @@
-##cobra.solvers.gurobi_solver
+# PLEASE NOTE THAT JYTHON SUPPORT (and this jython-only-solver) is deprecated
 #Interface to the gurobi 5.0.1 python and java solvers
 #QPs are not yet supported on java
+from __future__ import print_function
 from warnings import warn
 from os import name as __name
 from copy import deepcopy
+from six import iteritems
 ###solver specific parameters
 from .parameters import status_dict, variable_kind_dict, \
      sense_dict, parameter_mappings, parameter_defaults, \
@@ -70,7 +72,7 @@ def set_parameter(lp, parameter_name, parameter_value):
         else:
             warn("%s is not a DoubleParam, IntParam, StringParam, IntAttr"%parameter_name)
             ## raise Exception("%s is not a DoubleParam, IntParam, StringParam, IntAttr"%parameter_name)
-    except Exception, e:
+    except Exception as e:
         warn("%s %s didn't work %s"%(parameter_name, parameter_value, e))
 
 def get_objective_value(lp):
@@ -143,7 +145,7 @@ def create_problem(cobra_model,  **kwargs):
         the_parameters.update(kwargs)
 
     [set_parameter(lp, parameter_mappings[k], v)
-         for k, v in the_parameters.iteritems() if k in parameter_mappings]
+         for k, v in iteritems(the_parameters) if k in parameter_mappings]
     quadratic_component = the_parameters['quadratic_component']
     objective_sense = objective_senses[the_parameters['objective_sense']]
 
@@ -201,7 +203,7 @@ def solve_problem(lp, **kwargs):
     #Update parameter settings if provided
     if kwargs:
         [set_parameter(lp, parameter_mappings[k], v)
-         for k, v in kwargs.iteritems() if k in parameter_mappings]
+         for k, v in iteritems(kwargs) if k in parameter_mappings]
 
     try:
         print_solver_time = kwargs['print_solver_time']
@@ -213,7 +215,7 @@ def solve_problem(lp, **kwargs):
     lp.optimize()
     status = get_status(lp)
     if print_solver_time:
-        print 'optimize time: %f'%(time() - start_time)
+        print('optimize time: {:f}'.format(time() - start_time))
     return status
 
     
@@ -274,7 +276,7 @@ def solve(cobra_model, **kwargs):
     status = solve_problem(lp, **the_parameters)
     the_solution = format_solution(lp, cobra_model)
     if status != 'optimal' and error_reporting:
-        print '%s failed: %s'%(solver_name, status)
+        print('{:s} failed: {:s}'.format(solver_name, status))
     cobra_model.solution = the_solution
     solution = {'the_problem': lp, 'the_solution': the_solution}
     return solution
diff --git a/cobra/solvers/mosek.py b/cobra/solvers/mosek.py
index 077ab93..8078fa7 100644
--- a/cobra/solvers/mosek.py
+++ b/cobra/solvers/mosek.py
@@ -7,6 +7,12 @@ from six import iteritems, string_types
 
 env = mosek.Env()
 
+# make sure the mosek environment works
+test = env.Task(0, 0)
+test.optimize()
+del test
+
+
 solver_name = "mosek"
 __mosek_version__ = ".".join(str(i) for i in mosek.getversion())
 _SUPPORTS_MILP = True
diff --git a/cobra/test/__init__.py b/cobra/test/__init__.py
index 564fefc..fc59265 100644
--- a/cobra/test/__init__.py
+++ b/cobra/test/__init__.py
@@ -10,7 +10,8 @@ except:
 from ..io import read_sbml_model
 
 
-available_tests = ['unit_tests', 'solvers', 'flux_analysis', 'io_tests']
+available_tests = ['unit_tests', 'solvers', 'flux_analysis', 'io_tests',
+                   'design', 'manipulation']
 
 
 cobra_directory = abspath(join(dirname(abspath(__file__)), ".."))
diff --git a/cobra/test/data/invalid0.xml b/cobra/test/data/invalid0.xml
new file mode 100644
index 0000000..5400b79
--- /dev/null
+++ b/cobra/test/data/invalid0.xml
@@ -0,0 +1,47 @@
+<sbml xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" level="3" sboTerm="SBO:0000624" version="1" xmlns="http://www.sbml.org/sbml/level3/version1/core" fbc:required="false">
+  <model fbc:strict="true" id="invalid3">
+    <listOfUnitDefinitions>
+      <unitDefinition id="mmol_per_gDW_per_hr">
+        <listOfUnits>
+          <unit exponent="1" kind="mole" multiplier="1" scale="-3"/>
+          <unit exponent="-1" kind="gram" multiplier="1" scale="0"/>
+          <unit exponent="-1" kind="second" multiplier="3600" scale="0"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <fbc:listOfObjectives fbc:activeObjective="obj">
+      <fbc:objective fbc:id="obj" fbc:type="maximize">
+        <fbc:listOfFluxObjectives>
+          <fbc:fluxObjective fbc:reaction="R_ATPM" fbc:coefficient="1"/>
+        </fbc:listOfFluxObjectives>
+      </fbc:objective>
+    </fbc:listOfObjectives>
+    <listOfParameters>
+      <parameter constant="true" id="cobra_default_lb" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="8.39"/>
+      <parameter constant="true" id="cobra_default_ub" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="1000"/>
+      <parameter constant="true" id="cobra_0_bound" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="0"/>
+    </listOfParameters>
+    <listOfCompartments/>
+    <listOfSpecies>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_c" name="H2O" compartment="c" fbc:chemicalFormula="H2O"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_atp_c" name="ATP" compartment="c" fbc:chemicalFormula="C10H12N5O13P3"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_c" name="H+" compartment="c" fbc:chemicalFormula="H"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_adp_c" name="ADP" compartment="c" fbc:chemicalFormula="C10H12N5O10P2"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_c" name="Phosphate" compartment="c" fbc:chemicalFormula="HO4P"/>
+      <species boundaryCondition="true" id="Please_never_use_this"/>
+    </listOfSpecies>
+    <listOfReactions>
+      <reaction fast="false" id="R_ATPM" reversible="false" name="ATP maintenance requirement" metaid="R_ATPM" fbc:upperFluxBound="fake">
+        <listOfReactants>
+          <speciesReference constant="true" species="M_atp_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_h2o_c" stoichiometry="1"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference constant="true" species="M_adp_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_h_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_pi_c" stoichiometry="1"/>
+        </listOfProducts>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>
diff --git a/cobra/test/data/invalid1.xml b/cobra/test/data/invalid1.xml
new file mode 100644
index 0000000..b83d731
--- /dev/null
+++ b/cobra/test/data/invalid1.xml
@@ -0,0 +1,45 @@
+<sbml xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" level="3" sboTerm="SBO:0000624" version="1" xmlns="http://www.sbml.org/sbml/level3/version1/core" fbc:required="false">
+  <model fbc:strict="true" id="invalid1">
+    <listOfUnitDefinitions>
+      <unitDefinition id="mmol_per_gDW_per_hr">
+        <listOfUnits>
+          <unit exponent="1" kind="mole" multiplier="1" scale="-3"/>
+          <unit exponent="-1" kind="gram" multiplier="1" scale="0"/>
+          <unit exponent="-1" kind="second" multiplier="3600" scale="0"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <fbc:listOfObjectives fbc:activeObjective="obj">
+      <fbc:objective fbc:id="obj" fbc:type="maximize">
+        <fbc:listOfFluxObjectives>
+          <fbc:fluxObjective fbc:reaction="invalid" fbc:coefficient="1"/>
+        </fbc:listOfFluxObjectives>
+      </fbc:objective>
+    </fbc:listOfObjectives>
+    <listOfParameters>
+      <parameter constant="true" id="cobra_default_lb" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="-1000"/>
+      <parameter constant="true" id="cobra_default_ub" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="1000"/>
+      <parameter constant="true" id="cobra_0_bound" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="0"/>
+    </listOfParameters>
+    <listOfCompartments/>
+    <listOfSpecies>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id=" M_g6p_c" name="D-Glucose 6-phosphate" metaid="M_g6p_c" compartment="c" fbc:chemicalFormula="C6H11O9P">
+      </species>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="$M_f6p_c" name="D-Fructose 6-phosphate" metaid="M_f6p_c" compartment="c" fbc:chemicalFormula="C6H11O9P">
+      </species>
+    </listOfSpecies>
+    <listOfReactions>
+      <reaction fast="false" id="R_PGI" reversible="true" name="glucose-6-phosphate isomerase" metaid="R_PGI" fbc:upperFluxBound="cobra_default_ub" fbc:lowerFluxBound="cobra_default_lb">
+        <listOfReactants>
+          <speciesReference constant="true" species="M_g6p_c" stoichiometry="1"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference constant="true" species="M_f6p_c" stoichiometry="1"/>
+        </listOfProducts>
+        <fbc:geneProductAssociation>
+          <fbc:geneProductRef fbc:geneProduct="G_b4025"/>
+        </fbc:geneProductAssociation>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>
diff --git a/cobra/test/data/invalid2.xml b/cobra/test/data/invalid2.xml
new file mode 100644
index 0000000..23e8d73
--- /dev/null
+++ b/cobra/test/data/invalid2.xml
@@ -0,0 +1,68 @@
+<sbml xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" level="3" sboTerm="SBO:0000624" version="1" xmlns="http://www.sbml.org/sbml/level3/version1/core" fbc:required="false">
+  <model fbc:strict="true" id="invalid2">
+    <listOfUnitDefinitions>
+      <unitDefinition id="mmol_per_gDW_per_hr">
+        <listOfUnits>
+          <unit exponent="1" kind="mole" multiplier="1" scale="-3"/>
+          <unit exponent="-1" kind="gram" multiplier="1" scale="0"/>
+          <unit exponent="-1" kind="second" multiplier="3600" scale="0"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <fbc:listOfObjectives fbc:activeObjective="obj">
+      <fbc:objective fbc:id="obj" fbc:type="maximize">
+        <fbc:listOfFluxObjectives>
+          <fbc:fluxObjective fbc:reaction="R_ATPM" fbc:coefficient="1"/>
+        </fbc:listOfFluxObjectives>
+      </fbc:objective>
+    </fbc:listOfObjectives>
+    <listOfParameters>
+      <parameter constant="true" id="cobra_default_lb" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="8.39"/>
+      <parameter constant="true" id="cobra_default_ub" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="1000"/>
+      <parameter constant="true" id="cobra_0_bound" sboTerm="SBO:0000626" units="mmol_per_gDW_per_hr" value="0"/>
+    </listOfParameters>
+    <listOfCompartments/>
+    <listOfSpecies>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_c" metaid="M_h2o_c" name="H2O" compartment="c" fbc:chemicalFormula="H2O">
+        <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
+          <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
+            <rdf:Description rdf:about="#M_h2o_c">
+              <bqbiol:is xmlns:bqbiol="http://biomodels.net/biology-qualifiers/">
+                <rdf:Bag>
+                  <rdf:li rdf:resource="http://random.org/bigg.metabolite/h2o"/>
+                </rdf:Bag>
+              </bqbiol:is>
+              <bqbiol:isEncodedBy xmlns:bqbiol="http://biomodels.net/biology-qualifiers/">
+                <rdf:Bag>
+                  <rdf:li rdf:resource="http://identifiers.org/biocyc_WATER"/>
+                </rdf:Bag>
+              </bqbiol:isEncodedBy>
+            </rdf:Description>
+          </rdf:RDF>
+        </sbml:annotation>
+      </species>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_atp_c" name="ATP" compartment="c" fbc:chemicalFormula="C10H12N5O13P3"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_c" name="H+" compartment="c" fbc:chemicalFormula="H"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_adp_c" name="ADP" compartment="c" fbc:chemicalFormula="C10H12N5O10P2"/>
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_c" name="Phosphate" compartment="c" fbc:chemicalFormula="HO4P"/>
+      <species boundaryCondition="true" id="Please_never_use_this"/>
+    </listOfSpecies>
+    <listOfReactions>
+      <reaction fast="false" id="R_ATPM" reversible="false" name="ATP maintenance requirement" metaid="R_ATPM" fbc:upperFluxBound="cobra_default_lb" fbc:lowerFluxBound="cobra_default_ub">
+        <listOfReactants>
+          <speciesReference constant="true" species="M_atp_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_h2o_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="Please_never_use_this" stoichiometry="1"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference constant="true" species="M_adp_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_h_c" stoichiometry="1"/>
+          <speciesReference constant="true" species="M_pi_c" stoichiometry="1"/>
+        </listOfProducts>
+        <fbc:geneProductAssociation>
+          <fbc:geneProductRef fbc:geneProduct="G_b4025"/>
+        </fbc:geneProductAssociation>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>
diff --git a/cobra/test/data/mini.json b/cobra/test/data/mini.json
index 3e1d20b..85ea068 100644
--- a/cobra/test/data/mini.json
+++ b/cobra/test/data/mini.json
@@ -3,7 +3,6 @@
         "c": "cytosol",
         "e": "extracellular"
     },
-    "description": "mini_textbook",
     "genes": [
         {
             "id": "b0755",
@@ -171,6 +170,7 @@
                 "seed.compound": "cpd00203",
                 "unipathway.compound": "UPC00236"
             },
+            "charge": -4,
             "compartment": "c",
             "formula": "C3H4O10P2",
             "id": "13dpg_c",
@@ -200,6 +200,7 @@
                 "seed.compound": "cpd00482",
                 "unipathway.compound": "UPC00631"
             },
+            "charge": -3,
             "compartment": "c",
             "formula": "C3H4O7P",
             "id": "2pg_c",
@@ -237,6 +238,7 @@
                     "UPC00197"
                 ]
             },
+            "charge": -3,
             "compartment": "c",
             "formula": "C3H4O7P",
             "id": "3pg_c",
@@ -281,6 +283,7 @@
                 "seed.compound": "cpd00008",
                 "unipathway.compound": "UPC00008"
             },
+            "charge": -3,
             "compartment": "c",
             "formula": "C10H12N5O10P2",
             "id": "adp_c",
@@ -322,6 +325,7 @@
                 "seed.compound": "cpd00002",
                 "unipathway.compound": "UPC00002"
             },
+            "charge": -4,
             "compartment": "c",
             "formula": "C10H12N5O13P3",
             "id": "atp_c",
@@ -358,6 +362,7 @@
                 "seed.compound": "cpd00095",
                 "unipathway.compound": "UPC00111"
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "C3H5O6P",
             "id": "dhap_c",
@@ -394,6 +399,7 @@
                     "UPC00085"
                 ]
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "C6H11O9P",
             "id": "f6p_c",
@@ -429,6 +435,7 @@
                 "seed.compound": "cpd00290",
                 "unipathway.compound": "UPC00354"
             },
+            "charge": -4,
             "compartment": "c",
             "formula": "C6H10O12P2",
             "id": "fdp_c",
@@ -459,6 +466,7 @@
                     "UPC00118"
                 ]
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "C3H5O6P",
             "id": "g3p_c",
@@ -500,6 +508,7 @@
                 "seed.compound": "cpd00079",
                 "unipathway.compound": "UPC00092"
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "C6H11O9P",
             "id": "g6p_c",
@@ -515,6 +524,7 @@
                 "kegg.compound": "C00031",
                 "pubchem.substance": "3333"
             },
+            "charge": 0,
             "compartment": "e",
             "formula": "C6H12O6",
             "id": "glc__D_e",
@@ -600,6 +610,7 @@
                     "UPC01328"
                 ]
             },
+            "charge": 0,
             "compartment": "c",
             "formula": "H2O",
             "id": "h2o_c",
@@ -685,6 +696,7 @@
                     "UPC01328"
                 ]
             },
+            "charge": 0,
             "compartment": "e",
             "formula": "H2O",
             "id": "h2o_e",
@@ -731,6 +743,7 @@
                 "seed.compound": "cpd00067",
                 "unipathway.compound": "UPC00080"
             },
+            "charge": 1,
             "compartment": "c",
             "formula": "H",
             "id": "h_c",
@@ -777,6 +790,7 @@
                 "seed.compound": "cpd00067",
                 "unipathway.compound": "UPC00080"
             },
+            "charge": 1,
             "compartment": "e",
             "formula": "H",
             "id": "h_e",
@@ -830,6 +844,7 @@
                 "seed.compound": "cpd00003",
                 "unipathway.compound": "UPC00003"
             },
+            "charge": -1,
             "compartment": "c",
             "formula": "C21H26N7O14P2",
             "id": "nad_c",
@@ -864,6 +879,7 @@
                 "seed.compound": "cpd00004",
                 "unipathway.compound": "UPC00004"
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "C21H27N7O14P2",
             "id": "nadh_c",
@@ -897,6 +913,7 @@
                 "seed.compound": "cpd00061",
                 "unipathway.compound": "UPC00074"
             },
+            "charge": -3,
             "compartment": "c",
             "formula": "C3H2O6P",
             "id": "pep_c",
@@ -964,6 +981,7 @@
                 ],
                 "unipathway.compound": "UPC00009"
             },
+            "charge": -2,
             "compartment": "c",
             "formula": "HO4P",
             "id": "pi_c",
@@ -1031,6 +1049,7 @@
                 ],
                 "unipathway.compound": "UPC00009"
             },
+            "charge": -2,
             "compartment": "e",
             "formula": "HO4P",
             "id": "pi_e",
@@ -1065,6 +1084,7 @@
                 "seed.compound": "cpd00020",
                 "unipathway.compound": "UPC00022"
             },
+            "charge": -1,
             "compartment": "c",
             "formula": "C3H3O3",
             "id": "pyr_c",
@@ -1100,7 +1120,7 @@
                 "lac__D_c": 1,
                 "lac__D_e": -1
             },
-            "name": "D_LACt2",
+            "name": "",
             "upper_bound": 1000.0
         },
         {
diff --git a/cobra/test/data/mini.mat b/cobra/test/data/mini.mat
index 9e8b013..d0f5128 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 013b578..5b5cb8a 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 1d796d0..4a136b3 100644
--- a/cobra/test/data/mini_cobra.xml
+++ b/cobra/test/data/mini_cobra.xml
@@ -15,153 +15,199 @@
       <compartment id="e" name="extracellular" size="1"/>
     </listOfCompartments>
     <listOfSpecies>
-      <species id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" compartment="c">
+      <species id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" compartment="c" charge="-4">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O10P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_2pg_c" name="D-Glycerate 2-phosphate" compartment="c">
+      <species id="M_2pg_c" name="D-Glycerate 2-phosphate" compartment="c" charge="-3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O7P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_3pg_c" name="3-Phospho-D-glycerate" compartment="c">
+      <species id="M_3pg_c" name="3-Phospho-D-glycerate" compartment="c" charge="-3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O7P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_adp_c" name="ADP" compartment="c">
+      <species id="M_adp_c" name="ADP" compartment="c" charge="-3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C10H12N5O10P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_atp_c" name="ATP" compartment="c">
+      <species id="M_atp_c" name="ATP" compartment="c" charge="-4">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C10H12N5O13P3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_dhap_c" name="Dihydroxyacetone phosphate" compartment="c">
+      <species id="M_dhap_c" name="Dihydroxyacetone phosphate" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_f6p_c" name="D-Fructose 6-phosphate" compartment="c">
+      <species id="M_f6p_c" name="D-Fructose 6-phosphate" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H11O9P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" compartment="c">
+      <species id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" compartment="c" charge="-4">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H10O12P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_g3p_c" name="Glyceraldehyde 3-phosphate" compartment="c">
+      <species id="M_g3p_c" name="Glyceraldehyde 3-phosphate" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_g6p_c" name="D-Glucose 6-phosphate" compartment="c">
+      <species id="M_g6p_c" name="D-Glucose 6-phosphate" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H11O9P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_glc__D_e" name="D-Glucose" compartment="e">
+      <species id="M_glc__D_e" name="D-Glucose" compartment="e" charge="0">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H12O6</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h2o_c" name="H2O" compartment="c">
+      <species id="M_h2o_c" name="H2O" compartment="c" charge="0">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H2O</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h2o_e" name="H2O" compartment="e">
+      <species id="M_h2o_e" name="H2O" compartment="e" charge="0">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H2O</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_c" name="H+" compartment="c">
+      <species id="M_h_c" name="H+" compartment="c" charge="1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_e" name="H+" compartment="e">
+      <species id="M_h_e" name="H+" compartment="e" charge="1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
       <species id="M_lac__D_c" name="D-Lactate" compartment="c" charge="-1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
       <species id="M_lac__D_e" name="D-Lactate" compartment="e" charge="-1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_nad_c" name="Nicotinamide adenine dinucleotide" compartment="c">
+      <species id="M_nad_c" name="Nicotinamide adenine dinucleotide" compartment="c" charge="-1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C21H26N7O14P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" compartment="c">
+      <species id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C21H27N7O14P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pep_c" name="Phosphoenolpyruvate" compartment="c">
+      <species id="M_pep_c" name="Phosphoenolpyruvate" compartment="c" charge="-3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H2O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pi_c" name="Phosphate" compartment="c">
+      <species id="M_pi_c" name="Phosphate" compartment="c" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: HO4P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pi_e" name="Phosphate" compartment="e">
+      <species id="M_pi_e" name="Phosphate" compartment="e" charge="-2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: HO4P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pyr_c" name="Pyruvate" compartment="c">
+      <species id="M_pyr_c" name="Pyruvate" compartment="c" charge="-1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H3O3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_glc__D_e_boundary" name="D-Glucose" compartment="e" boundaryCondition="true">
+      <species id="M_glc__D_e_boundary" name="D-Glucose" compartment="e" boundaryCondition="true" charge="0">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H12O6</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_e_boundary" name="H+" compartment="e" boundaryCondition="true">
+      <species id="M_h_e_boundary" name="H+" compartment="e" boundaryCondition="true" charge="1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
       <species id="M_lac__D_e_boundary" name="D-Lactate" compartment="e" boundaryCondition="true" charge="-1">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
     </listOfSpecies>
     <listOfReactions>
       <reaction id="R_ATPM" name="ATP maintenance requirement" reversible="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
-          <speciesReference species="M_h2o_c" stoichiometry="1"/>
           <speciesReference species="M_atp_c" stoichiometry="1"/>
+          <speciesReference species="M_h2o_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_pi_c" stoichiometry="1"/>
-          <speciesReference species="M_h_c" stoichiometry="1"/>
           <speciesReference species="M_adp_c" stoichiometry="1"/>
+          <speciesReference species="M_h_c" stoichiometry="1"/>
+          <speciesReference species="M_pi_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -179,7 +225,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b3603 or b2975</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -206,15 +251,14 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2779</p>
-            <p>SUBSYSTEM: </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">
@@ -229,12 +273,6 @@
         </kineticLaw>
       </reaction>
       <reaction id="R_EX_glc__D_e" name="D-Glucose exchange" reversible="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_glc__D_e" stoichiometry="1"/>
         </listOfReactants>
@@ -254,12 +292,6 @@
         </kineticLaw>
       </reaction>
       <reaction id="R_EX_h_e" name="H+ exchange" reversible="true">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_h_e" stoichiometry="1"/>
         </listOfReactants>
@@ -279,12 +311,6 @@
         </kineticLaw>
       </reaction>
       <reaction id="R_EX_lac__D_e" name="D-lactate exchange" reversible="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_lac__D_e" stoichiometry="1"/>
         </listOfReactants>
@@ -307,15 +333,14 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b1773 or b2097 or b2925</p>
-            <p>SUBSYSTEM: </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">
@@ -333,13 +358,12 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b1779</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
+          <speciesReference species="M_g3p_c" stoichiometry="1"/>
           <speciesReference species="M_nad_c" stoichiometry="1"/>
           <speciesReference species="M_pi_c" stoichiometry="1"/>
-          <speciesReference species="M_g3p_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_h_c" stoichiometry="1"/>
@@ -362,7 +386,6 @@
         <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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -370,8 +393,8 @@
           <speciesReference species="M_glc__D_e" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_g6p_c" stoichiometry="1"/>
           <speciesReference species="M_pyr_c" stoichiometry="1"/>
+          <speciesReference species="M_g6p_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -389,7 +412,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b0875 or s0001</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -414,17 +436,16 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2133 or b1380</p>
-            <p>SUBSYSTEM: </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_h_c" stoichiometry="1"/>
           <speciesReference species="M_pyr_c" stoichiometry="1"/>
           <speciesReference species="M_nadh_c" stoichiometry="1"/>
-          <speciesReference species="M_h_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
@@ -442,7 +463,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b3916 or b1723</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -450,8 +470,8 @@
           <speciesReference species="M_atp_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_h_c" stoichiometry="1"/>
           <speciesReference species="M_adp_c" stoichiometry="1"/>
+          <speciesReference species="M_h_c" stoichiometry="1"/>
           <speciesReference species="M_fdp_c" stoichiometry="1"/>
         </listOfProducts>
         <kineticLaw>
@@ -470,7 +490,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b4025</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -495,12 +514,11 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2926</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_atp_c" stoichiometry="1"/>
           <speciesReference species="M_3pg_c" stoichiometry="1"/>
+          <speciesReference species="M_atp_c" stoichiometry="1"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_13dpg_c" stoichiometry="1"/>
@@ -522,7 +540,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b4395 or b3612 or b0755</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -547,7 +564,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2987 or b3493</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -574,7 +590,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b1854 or b1676</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -583,8 +598,8 @@
           <speciesReference species="M_pep_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">
@@ -602,7 +617,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b3919</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
diff --git a/cobra/test/data/mini_fbc1.xml b/cobra/test/data/mini_fbc1.xml
index 2a243e9..ae48616 100644
--- a/cobra/test/data/mini_fbc1.xml
+++ b/cobra/test/data/mini_fbc1.xml
@@ -122,160 +122,205 @@
       <compartment id="e" name="extracellular" spatialDimensions="3" size="1" units="volume" constant="true"/>
     </listOfCompartments>
     <listOfSpecies>
-      <species id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-4" fbc:chemicalFormula="C3H4O10P2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O10P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_2pg_c" name="D-Glycerate 2-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_2pg_c" name="D-Glycerate 2-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-3" fbc:chemicalFormula="C3H4O7P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O7P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_3pg_c" name="3-Phospho-D-glycerate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_3pg_c" name="3-Phospho-D-glycerate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-3" fbc:chemicalFormula="C3H4O7P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H4O7P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_adp_c" name="ADP" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_adp_c" name="ADP" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-3" fbc:chemicalFormula="C10H12N5O10P2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C10H12N5O10P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_atp_c" name="ATP" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_atp_c" name="ATP" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-4" fbc:chemicalFormula="C10H12N5O13P3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C10H12N5O13P3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_dhap_c" name="Dihydroxyacetone phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_dhap_c" name="Dihydroxyacetone phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="C3H5O6P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_f6p_c" name="D-Fructose 6-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_f6p_c" name="D-Fructose 6-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="C6H11O9P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H11O9P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-4" fbc:chemicalFormula="C6H10O12P2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H10O12P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_g3p_c" name="Glyceraldehyde 3-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_g3p_c" name="Glyceraldehyde 3-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="C3H5O6P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_g6p_c" name="D-Glucose 6-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_g6p_c" name="D-Glucose 6-phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="C6H11O9P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H11O9P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_glc__D_e" name="D-Glucose" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_glc__D_e" name="D-Glucose" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:chemicalFormula="C6H12O6">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H12O6</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h2o_c" name="H2O" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_h2o_c" name="H2O" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:chemicalFormula="H2O">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H2O</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h2o_e" name="H2O" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_h2o_e" name="H2O" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:chemicalFormula="H2O">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H2O</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_c" name="H+" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_h_c" name="H+" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="1" fbc:chemicalFormula="H">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_e" name="H+" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_h_e" name="H+" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="1" fbc:chemicalFormula="H">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
-      <species id="M_lac__D_c" name="D-Lactate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1">
+      <species id="M_lac__D_c" name="D-Lactate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1" fbc:chemicalFormula="C3H5O3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_lac__D_e" name="D-Lactate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1">
+      <species id="M_lac__D_e" name="D-Lactate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1" fbc:chemicalFormula="C3H5O3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_nad_c" name="Nicotinamide adenine dinucleotide" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_nad_c" name="Nicotinamide adenine dinucleotide" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1" fbc:chemicalFormula="C21H26N7O14P2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C21H26N7O14P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="C21H27N7O14P2">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C21H27N7O14P2</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pep_c" name="Phosphoenolpyruvate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_pep_c" name="Phosphoenolpyruvate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-3" fbc:chemicalFormula="C3H2O6P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H2O6P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pi_c" name="Phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_pi_c" name="Phosphate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="HO4P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: HO4P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pi_e" name="Phosphate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_pi_e" name="Phosphate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-2" fbc:chemicalFormula="HO4P">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: HO4P</p>
+          </html>
         </notes>
       </species>
-      <species id="M_pyr_c" name="Pyruvate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
+      <species id="M_pyr_c" name="Pyruvate" compartment="c" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false" fbc:charge="-1" fbc:chemicalFormula="C3H3O3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H3O3</p>
+          </html>
         </notes>
       </species>
-      <species id="M_glc__D_e_boundary" name="D-Glucose" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false">
+      <species id="M_glc__D_e_boundary" name="D-Glucose" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false" fbc:chemicalFormula="C6H12O6">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C6H12O6</p>
+          </html>
         </notes>
       </species>
-      <species id="M_h_e_boundary" name="H+" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false">
+      <species id="M_h_e_boundary" name="H+" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false" fbc:charge="1" fbc:chemicalFormula="H">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: H</p>
+          </html>
         </notes>
       </species>
-      <species id="M_lac__D_e_boundary" name="D-Lactate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false" fbc:charge="-1">
+      <species id="M_lac__D_e_boundary" name="D-Lactate" compartment="e" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="true" constant="false" fbc:charge="-1" fbc:chemicalFormula="C3H5O3">
         <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml"/>
+          <html xmlns="http://www.w3.org/1999/xhtml">
+            <p>FORMULA: C3H5O3</p>
+          </html>
         </notes>
       </species>
     </listOfSpecies>
     <listOfReactions>
       <reaction id="R_ATPM" name="ATP maintenance requirement" reversible="false" fast="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
-          <speciesReference species="M_h2o_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_h2o_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_pi_c" stoichiometry="1" constant="true"/>
-          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_adp_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_pi_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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -291,24 +336,17 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2779</p>
-            <p>SUBSYSTEM: </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">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_glc__D_e" stoichiometry="1" constant="true"/>
         </listOfReactants>
@@ -317,12 +355,6 @@
         </listOfProducts>
       </reaction>
       <reaction id="R_EX_h_e" name="H+ exchange" reversible="true" fast="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_h_e" stoichiometry="1" constant="true"/>
         </listOfReactants>
@@ -331,12 +363,6 @@
         </listOfProducts>
       </reaction>
       <reaction id="R_EX_lac__D_e" name="D-lactate exchange" reversible="false" fast="false">
-        <notes>
-          <html xmlns="http://www.w3.org/1999/xhtml">
-            <p>GENE_ASSOCIATION: </p>
-            <p>SUBSYSTEM: </p>
-          </html>
-        </notes>
         <listOfReactants>
           <speciesReference species="M_lac__D_e" stoichiometry="1" constant="true"/>
         </listOfReactants>
@@ -348,28 +374,26 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b1773 or b2097 or b2925</p>
-            <p>SUBSYSTEM: </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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
+          <speciesReference species="M_g3p_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_nad_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pi_c" stoichiometry="1" constant="true"/>
-          <speciesReference species="M_g3p_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
@@ -381,7 +405,6 @@
         <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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -389,15 +412,14 @@
           <speciesReference species="M_glc__D_e" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_g6p_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pyr_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_g6p_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -411,24 +433,22 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2133 or b1380</p>
-            <p>SUBSYSTEM: </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_h_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_pyr_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_nadh_c" stoichiometry="1" constant="true"/>
-          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
         </listOfProducts>
       </reaction>
       <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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -436,8 +456,8 @@
           <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
-          <speciesReference species="M_h_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_adp_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>
@@ -445,7 +465,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b4025</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -459,12 +478,11 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2926</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
-          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
           <speciesReference species="M_3pg_c" stoichiometry="1" constant="true"/>
+          <speciesReference species="M_atp_c" stoichiometry="1" constant="true"/>
         </listOfReactants>
         <listOfProducts>
           <speciesReference species="M_13dpg_c" stoichiometry="1" constant="true"/>
@@ -475,7 +493,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b4395 or b3612 or b0755</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -489,7 +506,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b2987 or b3493</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -505,7 +521,6 @@
         <notes>
           <html xmlns="http://www.w3.org/1999/xhtml">
             <p>GENE_ASSOCIATION: b1854 or b1676</p>
-            <p>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
@@ -514,15 +529,14 @@
           <speciesReference species="M_pep_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>SUBSYSTEM: </p>
           </html>
         </notes>
         <listOfReactants>
diff --git a/cobra/test/data/mini_fbc2.xml b/cobra/test/data/mini_fbc2.xml
index 8ee34bd..fb0f57f 100644
--- a/cobra/test/data/mini_fbc2.xml
+++ b/cobra/test/data/mini_fbc2.xml
@@ -28,7 +28,7 @@
       <compartment constant="true" id="e" name="extracellular"/>
     </listOfCompartments>
     <listOfSpecies>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" metaid="M_13dpg_c" compartment="c" fbc:chemicalFormula="C3H4O10P2">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_13dpg_c" name="3-Phospho-D-glyceroyl phosphate" metaid="M_13dpg_c" compartment="c" fbc:charge="-4" fbc:chemicalFormula="C3H4O10P2">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_13dpg_c">
@@ -53,7 +53,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_2pg_c" name="D-Glycerate 2-phosphate" metaid="M_2pg_c" compartment="c" fbc:chemicalFormula="C3H4O7P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_2pg_c" name="D-Glycerate 2-phosphate" metaid="M_2pg_c" compartment="c" fbc:charge="-3" fbc:chemicalFormula="C3H4O7P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_2pg_c">
@@ -82,7 +82,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_3pg_c" name="3-Phospho-D-glycerate" metaid="M_3pg_c" compartment="c" fbc:chemicalFormula="C3H4O7P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_3pg_c" name="3-Phospho-D-glycerate" metaid="M_3pg_c" compartment="c" fbc:charge="-3" fbc:chemicalFormula="C3H4O7P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_3pg_c">
@@ -117,7 +117,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_adp_c" name="ADP" metaid="M_adp_c" compartment="c" fbc:chemicalFormula="C10H12N5O10P2">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_adp_c" name="ADP" metaid="M_adp_c" compartment="c" fbc:charge="-3" fbc:chemicalFormula="C10H12N5O10P2">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_adp_c">
@@ -157,7 +157,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_atp_c" name="ATP" metaid="M_atp_c" compartment="c" fbc:chemicalFormula="C10H12N5O13P3">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_atp_c" name="ATP" metaid="M_atp_c" compartment="c" fbc:charge="-4" fbc:chemicalFormula="C10H12N5O13P3">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_atp_c">
@@ -196,7 +196,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_dhap_c" name="Dihydroxyacetone phosphate" metaid="M_dhap_c" compartment="c" fbc:chemicalFormula="C3H5O6P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_dhap_c" name="Dihydroxyacetone phosphate" metaid="M_dhap_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="C3H5O6P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_dhap_c">
@@ -228,7 +228,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_f6p_c" name="D-Fructose 6-phosphate" metaid="M_f6p_c" compartment="c" fbc:chemicalFormula="C6H11O9P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_f6p_c" name="D-Fructose 6-phosphate" metaid="M_f6p_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="C6H11O9P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_f6p_c">
@@ -260,7 +260,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" metaid="M_fdp_c" compartment="c" fbc:chemicalFormula="C6H10O12P2">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_fdp_c" name="D-Fructose 1,6-bisphosphate" metaid="M_fdp_c" compartment="c" fbc:charge="-4" fbc:chemicalFormula="C6H10O12P2">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_fdp_c">
@@ -293,7 +293,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_g3p_c" name="Glyceraldehyde 3-phosphate" metaid="M_g3p_c" compartment="c" fbc:chemicalFormula="C3H5O6P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_g3p_c" name="Glyceraldehyde 3-phosphate" metaid="M_g3p_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="C3H5O6P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_g3p_c">
@@ -319,7 +319,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_g6p_c" name="D-Glucose 6-phosphate" metaid="M_g6p_c" compartment="c" fbc:chemicalFormula="C6H11O9P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_g6p_c" name="D-Glucose 6-phosphate" metaid="M_g6p_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="C6H11O9P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_g6p_c">
@@ -354,7 +354,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_glc__D_e" name="D-Glucose" metaid="M_glc__D_e" compartment="e" fbc:chemicalFormula="C6H12O6">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_glc__D_e" name="D-Glucose" metaid="M_glc__D_e" compartment="e" fbc:charge="0" fbc:chemicalFormula="C6H12O6">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_glc__D_e">
@@ -371,7 +371,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_c" name="H2O" metaid="M_h2o_c" compartment="c" fbc:chemicalFormula="H2O">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_c" name="H2O" metaid="M_h2o_c" compartment="c" fbc:charge="0" fbc:chemicalFormula="H2O">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_h2o_c">
@@ -442,7 +442,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_e" name="H2O" metaid="M_h2o_e" compartment="e" fbc:chemicalFormula="H2O">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h2o_e" name="H2O" metaid="M_h2o_e" compartment="e" fbc:charge="0" fbc:chemicalFormula="H2O">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_h2o_e">
@@ -513,7 +513,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_c" name="H+" metaid="M_h_c" compartment="c" fbc:chemicalFormula="H">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_c" name="H+" metaid="M_h_c" compartment="c" fbc:charge="1" fbc:chemicalFormula="H">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_h_c">
@@ -557,7 +557,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_e" name="H+" metaid="M_h_e" compartment="e" fbc:chemicalFormula="H">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_h_e" name="H+" metaid="M_h_e" compartment="e" fbc:charge="1" fbc:chemicalFormula="H">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_h_e">
@@ -603,7 +603,7 @@
       </species>
       <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_lac__D_c" name="D-Lactate" compartment="c" fbc:charge="-1" fbc:chemicalFormula="C3H5O3"/>
       <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_lac__D_e" name="D-Lactate" compartment="e" fbc:charge="-1" fbc:chemicalFormula="C3H5O3"/>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_nad_c" name="Nicotinamide adenine dinucleotide" metaid="M_nad_c" compartment="c" fbc:chemicalFormula="C21H26N7O14P2">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_nad_c" name="Nicotinamide adenine dinucleotide" metaid="M_nad_c" compartment="c" fbc:charge="-1" fbc:chemicalFormula="C21H26N7O14P2">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_nad_c">
@@ -640,7 +640,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" metaid="M_nadh_c" compartment="c" fbc:chemicalFormula="C21H27N7O14P2">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_nadh_c" name="Nicotinamide adenine dinucleotide - reduced" metaid="M_nadh_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="C21H27N7O14P2">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_nadh_c">
@@ -672,7 +672,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pep_c" name="Phosphoenolpyruvate" metaid="M_pep_c" compartment="c" fbc:chemicalFormula="C3H2O6P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pep_c" name="Phosphoenolpyruvate" metaid="M_pep_c" compartment="c" fbc:charge="-3" fbc:chemicalFormula="C3H2O6P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_pep_c">
@@ -703,7 +703,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_c" name="Phosphate" metaid="M_pi_c" compartment="c" fbc:chemicalFormula="HO4P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_c" name="Phosphate" metaid="M_pi_c" compartment="c" fbc:charge="-2" fbc:chemicalFormula="HO4P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_pi_c">
@@ -762,7 +762,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_e" name="Phosphate" metaid="M_pi_e" compartment="e" fbc:chemicalFormula="HO4P">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pi_e" name="Phosphate" metaid="M_pi_e" compartment="e" fbc:charge="-2" fbc:chemicalFormula="HO4P">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_pi_e">
@@ -821,7 +821,7 @@
           </rdf:RDF>
         </sbml:annotation>
       </species>
-      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pyr_c" name="Pyruvate" metaid="M_pyr_c" compartment="c" fbc:chemicalFormula="C3H3O3">
+      <species boundaryCondition="false" constant="false" hasOnlySubstanceUnits="false" id="M_pyr_c" name="Pyruvate" metaid="M_pyr_c" compartment="c" fbc:charge="-1" fbc:chemicalFormula="C3H3O3">
         <sbml:annotation xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
           <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
             <rdf:Description rdf:about="#M_pyr_c">
@@ -972,7 +972,7 @@
           <speciesReference constant="true" species="M_pi_c" stoichiometry="1"/>
         </listOfProducts>
       </reaction>
-      <reaction fast="false" id="R_D_LACt2" reversible="true" name="D_LACt2" fbc:upperFluxBound="cobra_default_ub" fbc:lowerFluxBound="cobra_default_lb">
+      <reaction fast="false" id="R_D_LACt2" reversible="true" fbc:upperFluxBound="cobra_default_ub" fbc:lowerFluxBound="cobra_default_lb">
         <listOfReactants>
           <speciesReference constant="true" species="M_h_e" stoichiometry="1"/>
           <speciesReference constant="true" species="M_lac__D_e" stoichiometry="1"/>
diff --git a/cobra/test/data/mini_fbc2.xml.bz2 b/cobra/test/data/mini_fbc2.xml.bz2
index f4efa1b..b121d7d 100644
Binary files a/cobra/test/data/mini_fbc2.xml.bz2 and b/cobra/test/data/mini_fbc2.xml.bz2 differ
diff --git a/cobra/test/data/mini_fbc2.xml.gz b/cobra/test/data/mini_fbc2.xml.gz
index 677cf08..7ce91bc 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/textbook.xml.gz b/cobra/test/data/textbook.xml.gz
index b61570f..a5d1fbc 100644
Binary files a/cobra/test/data/textbook.xml.gz and b/cobra/test/data/textbook.xml.gz differ
diff --git a/cobra/test/design.py b/cobra/test/design.py
new file mode 100644
index 0000000..0359667
--- /dev/null
+++ b/cobra/test/design.py
@@ -0,0 +1,89 @@
+from unittest import TestCase, TestLoader, TextTestRunner, skipIf
+
+import sys
+
+if __name__ == "__main__":
+    sys.path.insert(0, "../..")
+    from cobra.test import create_test_model, data_directory
+    from cobra.design import *
+    from cobra.design.design_algorithms import _add_decision_variable
+    from cobra.solvers import get_solver_name
+    sys.path.pop(0)
+else:
+    from . import create_test_model, data_directory
+    from ..design import *
+    from ..design.design_algorithms import _add_decision_variable
+    from ..solvers import get_solver_name
+
+try:
+    solver = get_solver_name(mip=True)
+except:
+    no_mip_solver = True
+else:
+    no_mip_solver = False
+
+
+class TestDesignAlgorithms(TestCase):
+    """Test functions in cobra.design"""
+
+    def test_dual(self):
+        model = create_test_model("textbook")
+        self.assertAlmostEqual(model.optimize("maximize").f, 0.874, places=3)
+        dual = dual_problem(model)
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
+
+    def test_dual_integer_vars_as_lp(self):
+        model = create_test_model("textbook")
+        var = _add_decision_variable(model, "AKGDH")
+        self.assertAlmostEqual(model.optimize("maximize").f, 0.874, places=3)
+        # as lp: make integer continuous, set to 1
+        dual = dual_problem(model, "maximize", [var.id], copy=True)
+        r = dual.reactions.get_by_id(var.id)
+        r.variable_kind = "continuous"
+        r.lower_bound = r.upper_bound = 1
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
+        r.lower_bound = r.upper_bound = 0
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
+
+    @skipIf(no_mip_solver, "no MILP solver found")
+    def test_dual_integer_vars_as_mip(self):
+        # mip
+        model = create_test_model("textbook")
+        var = _add_decision_variable(model, "AKGDH")
+        dual = dual_problem(model, "maximize", [var.id], copy=True)
+        var_in_dual = dual.reactions.get_by_id(var.id)
+
+        # minimization, so the optimal value state is to turn off AKGDH
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
+
+        # turn off AKGDH in dual
+        var_in_dual.lower_bound = var_in_dual.upper_bound = 1
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
+
+        # turn on AKGDH in dual
+        var_in_dual.lower_bound = var_in_dual.upper_bound = 0
+        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
+
+    @skipIf(no_mip_solver, "no MILP solver found")
+    def test_optknock(self):
+        model = create_test_model("textbook")
+        model.reactions.get_by_id("EX_o2_e").lower_bound = 0
+        knockable_reactions = ["ACKr", "AKGDH", "ACALD", "LDH_D"]
+        optknock_problem = set_up_optknock(model, "EX_lac__D_e",
+                                           knockable_reactions, n_knockouts=2,
+                                           copy=False)
+        solution = run_optknock(optknock_problem, tolerance_integer=1e-9)
+        self.assertIn("ACKr", solution.knockouts)
+        self.assertIn("ACALD", solution.knockouts)
+        self.assertAlmostEqual(solution.f, 17.891, places=3)
+
+# make a test suite to run all of the tests
+loader = TestLoader()
+suite = loader.loadTestsFromModule(sys.modules[__name__])
+
+
+def test_all():
+    TextTestRunner(verbosity=2).run(suite)
+
+if __name__ == "__main__":
+    test_all()
diff --git a/cobra/test/flux_analysis.py b/cobra/test/flux_analysis.py
index 14a354f..4103cf0 100644
--- a/cobra/test/flux_analysis.py
+++ b/cobra/test/flux_analysis.py
@@ -11,6 +11,10 @@ try:
     import numpy
 except:
     numpy = None
+try:
+    import matplotlib
+except:
+    matplotlib = None
 
 if __name__ == "__main__":
     sys.path.insert(0, "../..")
@@ -18,7 +22,6 @@ if __name__ == "__main__":
     from cobra import Model, Reaction, Metabolite
     from cobra.manipulation import initialize_growth_medium
     from cobra.solvers import solver_dict, get_solver_name
-    from cobra.manipulation import modify, delete
     from cobra.flux_analysis import *
     sys.path.pop(0)
 else:
@@ -26,7 +29,6 @@ else:
     from .. import Model, Reaction, Metabolite
     from ..manipulation import initialize_growth_medium
     from ..solvers import solver_dict, get_solver_name
-    from ..manipulation import modify, delete
     from ..flux_analysis import *
 
 
@@ -45,113 +47,6 @@ class TestCobraFluxAnalysis(TestCase):
             self.assertAlmostEqual(model.solution.f, 0.8739, places=3)
             self.assertAlmostEqual(sum(abs_x), 518.4221, places=3)
 
-    def test_modify_reversible(self):
-        model1 = create_test_model("textbook")
-        model1.optimize()
-        model2 = create_test_model("textbook")
-        modify.convert_to_irreversible(model2)
-        model2.optimize()
-        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
-        modify.revert_to_reversible(model2)
-        model2.optimize()
-        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
-
-        # Ensure revert_to_reversible is robust to solutions generated both
-        # before and after reversibility conversion, or not solved at all.
-        model3 = create_test_model("textbook")
-        model3.optimize()
-        modify.convert_to_irreversible(model3)
-        modify.revert_to_reversible(model3)
-        self.assertAlmostEqual(model1.solution.f, model3.solution.f, places=3)
-
-        model4 = create_test_model("textbook")
-        modify.convert_to_irreversible(model4)
-        modify.revert_to_reversible(model4)
-
-    def test_escape_ids(self):
-        model = create_test_model('textbook')
-        model.reactions.PGI.gene_reaction_rule = "a.b or c"
-        self.assertIn("a.b", model.genes)
-        modify.escape_ID(model)
-        self.assertNotIn("a.b", model.genes)
-
-    def test_gene_knockout_computation(self):
-        cobra_model = create_test_model()
-
-        # helper functions for running tests
-        delete_model_genes = delete.delete_model_genes
-        find_gene_knockout_reactions = delete.find_gene_knockout_reactions
-
-        def find_gene_knockout_reactions_fast(cobra_model, gene_list):
-            compiled_rules = delete.get_compiled_gene_reaction_rules(
-                cobra_model)
-            return find_gene_knockout_reactions(
-                cobra_model, gene_list,
-                compiled_gene_reaction_rules=compiled_rules)
-
-        def get_removed(m):
-            return {x.id for x in m._trimmed_reactions}
-
-        def test_computation(m, gene_ids, expected_reaction_ids):
-            genes = [m.genes.get_by_id(i) for i in gene_ids]
-            expected_reactions = {m.reactions.get_by_id(i)
-                                  for i in expected_reaction_ids}
-            removed1 = set(find_gene_knockout_reactions(m, genes))
-            removed2 = set(find_gene_knockout_reactions_fast(m, genes))
-            self.assertEqual(removed1, expected_reactions)
-            self.assertEqual(removed2, expected_reactions)
-            delete.delete_model_genes(m, gene_ids, cumulative_deletions=False)
-            self.assertEqual(get_removed(m), expected_reaction_ids)
-            delete.undelete_model_genes(m)
-
-        gene_list = ['STM1067', 'STM0227']
-        dependent_reactions = {'3HAD121', '3HAD160', '3HAD80', '3HAD140',
-                               '3HAD180', '3HAD100', '3HAD181', '3HAD120',
-                               '3HAD60', '3HAD141', '3HAD161', 'T2DECAI',
-                               '3HAD40'}
-        test_computation(cobra_model, gene_list, dependent_reactions)
-        test_computation(cobra_model, ['STM4221'], {'PGI'})
-        test_computation(cobra_model, ['STM1746.S'], {'4PEPTabcpp'})
-        # test cumulative behavior
-        delete_model_genes(cobra_model, gene_list[:1])
-        delete_model_genes(cobra_model, gene_list[1:],
-                           cumulative_deletions=True)
-        delete_model_genes(cobra_model, ["STM4221"],
-                           cumulative_deletions=True)
-        dependent_reactions.add('PGI')
-        self.assertEqual(get_removed(cobra_model), dependent_reactions)
-        # non-cumulative following cumulative
-        delete_model_genes(cobra_model, ["STM4221"],
-                           cumulative_deletions=False)
-        self.assertEqual(get_removed(cobra_model), {'PGI'})
-        # make sure on reset that the bounds are correct
-        reset_bound = cobra_model.reactions.get_by_id("T2DECAI").upper_bound
-        self.assertEqual(reset_bound, 1000.)
-        # test computation when gene name is a subset of another
-        test_model = Model()
-        test_reaction_1 = Reaction("test1")
-        test_reaction_1.gene_reaction_rule = "eggs or (spam and eggspam)"
-        test_model.add_reaction(test_reaction_1)
-        test_computation(test_model, ["eggs"], set())
-        test_computation(test_model, ["eggs", "spam"], {'test1'})
-        # test computation with nested boolean expression
-        test_reaction_1.gene_reaction_rule = \
-            "g1 and g2 and (g3 or g4 or (g5 and g6))"
-        test_computation(test_model, ["g3"], set())
-        test_computation(test_model, ["g1"], {'test1'})
-        test_computation(test_model, ["g5"], set())
-        test_computation(test_model, ["g3", "g4", "g5"], {'test1'})
-        # test computation when gene names are python expressions
-        test_reaction_1.gene_reaction_rule = "g1 and (for or in)"
-        test_computation(test_model, ["for", "in"], {'test1'})
-        test_computation(test_model, ["for"], set())
-        test_reaction_1.gene_reaction_rule = "g1 and g2 and g2.conjugate"
-        test_computation(test_model, ["g2"], {"test1"})
-        test_computation(test_model, ["g2.conjugate"], {"test1"})
-        test_reaction_1.gene_reaction_rule = "g1 and (try:' or 'except:1)"
-        test_computation(test_model, ["try:'"], set())
-        test_computation(test_model, ["try:'", "'except:1"], {"test1"})
-
     def test_single_gene_deletion(self):
         cobra_model = create_test_model("textbook")
         # expected knockouts for textbook model
@@ -342,46 +237,7 @@ class TestCobraFluxAnalysis(TestCase):
         self.assertEqual({i[0].id for i in result},
                          {"SMILEY_EX_b", "SMILEY_EX_c"})
 
-    def test_remove_genes(self):
-        m = Model("test")
-        m.add_reactions([Reaction("r" + str(i + 1)) for i in range(8)])
-        self.assertEqual(len(m.reactions), 8)
-        rxns = m.reactions
-        rxns.r1.gene_reaction_rule = "(a and b) or (c and a)"
-        rxns.r2.gene_reaction_rule = "(a and b and d and e)"
-        rxns.r3.gene_reaction_rule = "(a and b) or (b and c)"
-        rxns.r4.gene_reaction_rule = "(f and b) or (b and c)"
-        rxns.r5.gene_reaction_rule = "x"
-        rxns.r6.gene_reaction_rule = "y"
-        rxns.r7.gene_reaction_rule = "x or     z"
-        rxns.r8.gene_reaction_rule = ""
-        self.assertIn("a", m.genes)
-        self.assertIn("x", m.genes)
-        delete.remove_genes(m, ["a"], remove_reactions=False)
-        self.assertNotIn("a", m.genes)
-        self.assertIn("x", m.genes)
-        self.assertEqual(rxns.r1.gene_reaction_rule, "")
-        self.assertEqual(rxns.r2.gene_reaction_rule, "")
-        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
-        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
-        self.assertEqual(rxns.r5.gene_reaction_rule, "x")
-        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
-        self.assertEqual(rxns.r7.genes, {m.genes.x, m.genes.z})
-        self.assertEqual(rxns.r8.gene_reaction_rule, "")
-        delete.remove_genes(m, ["x"], remove_reactions=True)
-        self.assertEqual(len(m.reactions), 7)
-        self.assertNotIn("r5", m.reactions)
-        self.assertNotIn("x", m.genes)
-        self.assertEqual(rxns.r1.gene_reaction_rule, "")
-        self.assertEqual(rxns.r2.gene_reaction_rule, "")
-        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
-        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
-        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
-        self.assertEqual(rxns.r7.gene_reaction_rule, "z")
-        self.assertEqual(rxns.r7.genes, {m.genes.z})
-        self.assertEqual(rxns.r8.gene_reaction_rule, "")
-
-    @skipIf(numpy is None, "double deletions require numpy")
+    @skipIf(numpy is None, "phase plane requires numpy")
     def test_phenotype_phase_plane(self):
         model = create_test_model("textbook")
         data = calculate_phenotype_phase_plane(
@@ -390,6 +246,9 @@ class TestCobraFluxAnalysis(TestCase):
         self.assertEqual(data.growth_rates.shape, (20, 20))
         self.assertAlmostEqual(data.growth_rates.max(), 1.20898, places=4)
         self.assertAlmostEqual(abs(data.growth_rates[0, :]).max(), 0, places=4)
+        if matplotlib is None:
+            self.skip("can't test plots without matplotlib")
+        data.plot()
 
 # make a test suite to run all of the tests
 loader = TestLoader()
diff --git a/cobra/test/io_tests.py b/cobra/test/io_tests.py
index 515d28c..05b237e 100644
--- a/cobra/test/io_tests.py
+++ b/cobra/test/io_tests.py
@@ -2,7 +2,8 @@ from warnings import catch_warnings, warn
 from tempfile import gettempdir
 from os import unlink
 from os.path import join, split
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf
+from unittest import TestCase, TestLoader, TextTestRunner, skipIf, \
+    expectedFailure
 from functools import partial
 from pickle import load, dump
 import sys
@@ -31,14 +32,15 @@ class TestCobraIO(object):
                          len(model2.reactions))
         self.assertEqual(len(model1.metabolites),
                          len(model2.metabolites))
-        for attr in ("id", "name", "lower_bound", "upper_bound"):
+        for attr in ("id", "name", "lower_bound", "upper_bound",
+                     "objective_coefficient", "gene_reaction_rule"):
             self.assertEqual(getattr(model1.reactions[0], attr),
                              getattr(model2.reactions[0], attr))
             self.assertEqual(getattr(model1.reactions[10], attr),
                              getattr(model2.reactions[10], attr))
             self.assertEqual(getattr(model1.reactions[-1], attr),
                              getattr(model2.reactions[-1], attr))
-        for attr in ("id", "name", "compartment"):
+        for attr in ("id", "name", "compartment", "formula", "charge"):
             self.assertEqual(getattr(model1.metabolites[0], attr),
                              getattr(model2.metabolites[0], attr))
             self.assertEqual(getattr(model1.metabolites[10], attr),
@@ -94,6 +96,19 @@ class TestCobraIO(object):
         self.validate(test_output_filename)
         unlink(test_output_filename)
 
+    def test_write_empty(self):
+        test_output_filename = join(gettempdir(), split(self.test_file)[-1])
+        m = self.test_model.copy()
+        m.metabolites[0].charge = None
+        m.remove_reactions(list(m.reactions))
+        self.write_function(m, test_output_filename)
+        reread_model = self.read_function(test_output_filename)
+        self.assertEqual(len(reread_model.reactions), 0)
+        self.assertEqual(len(reread_model.metabolites), len(m.metabolites))
+        # ensure empty metabolite charge is read as None
+        self.assertIs(reread_model.metabolites[0].charge, None)
+        unlink(test_output_filename)
+
     def validate(self, filename):
         # overload if a validator exists
         None
@@ -126,6 +141,14 @@ class TestCobraIOSBMLfbc2Bz2(TestCobraIOSBMLfbc2):
         self.write_function = io.write_sbml_model
 
 
+class TestCobraSBMLValidation(TestCase):
+    def test_bad_valiation(self):
+        for i in range(3):
+            filename = join(data_directory, "invalid%d.xml" % i)
+            m, errors = io.sbml3.validate_sbml_model(filename)
+            self.assertTrue(len(errors) >= 1)
+
+
 @skipIf(not libsbml, "libsbml required")
 class TestCobraIOSBMLfbc1(TestCase, TestCobraIO):
     def setUp(self):
@@ -138,6 +161,14 @@ class TestCobraIOSBMLfbc1(TestCase, TestCobraIO):
     def extra_comparisons(self, model1, model2):
         None
 
+    @expectedFailure
+    def test_read(self):
+        TestCobraIO.test_read(self)
+
+    @expectedFailure
+    def test_write(self):
+        TestCobraIO.test_write(self)
+
 
 @skipIf(not libsbml, "libsbml required")
 class TestCobraIOSBMLcobra(TestCase, TestCobraIO):
diff --git a/cobra/test/manipulation.py b/cobra/test/manipulation.py
new file mode 100644
index 0000000..3bca63d
--- /dev/null
+++ b/cobra/test/manipulation.py
@@ -0,0 +1,250 @@
+from unittest import TestCase, TestLoader, TextTestRunner
+
+import sys
+
+if __name__ == "__main__":
+    sys.path.insert(0, "../..")
+    from cobra.test import create_test_model, data_directory
+    from cobra.core import Metabolite, Model, Reaction
+    from cobra.manipulation import *
+    sys.path.pop(0)
+else:
+    from . import create_test_model, data_directory
+    from ..core import Metabolite, Model, Reaction
+    from ..manipulation import *
+
+
+class TestManipulation(TestCase):
+    """Test functions in cobra.manipulation"""
+
+    def test_canonical_form(self):
+        model = create_test_model("textbook")
+        # add G constraint to test
+        g_constr = Metabolite("SUCCt2_2__test_G_constraint")
+        g_constr._constraint_sense = "G"
+        g_constr._bound = 5.0
+        model.reactions.get_by_id("SUCCt2_2").add_metabolites({g_constr: 1})
+        self.assertAlmostEqual(model.optimize("maximize").f, 0.855, places=3)
+        # convert to canonical form
+        model = canonical_form(model)
+        self.assertAlmostEqual(model.optimize("maximize").f, 0.855, places=3)
+
+    def test_canonical_form_minimize(self):
+        model = create_test_model("textbook")
+        # make a minimization problem
+        model.reactions.get_by_id("Biomass_Ecoli_core").lower_bound = 0.5
+        for reaction in model.reactions:
+            reaction.objective_coefficient = reaction.id == "GAPD"
+        self.assertAlmostEqual(model.optimize("minimize").f, 6.27, places=3)
+        # convert to canonical form. Convert minimize to maximize
+        model = canonical_form(model, objective_sense="minimize")
+        self.assertAlmostEqual(model.optimize("maximize").f, -6.27, places=3)
+        # lower bounds should now be <= constraints
+        self.assertEqual(
+            model.reactions.get_by_id("Biomass_Ecoli_core").lower_bound, 0.0)
+
+    def test_modify_reversible(self):
+        model1 = create_test_model("textbook")
+        model1.optimize()
+        model2 = create_test_model("textbook")
+        convert_to_irreversible(model2)
+        model2.optimize()
+        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
+        revert_to_reversible(model2)
+        model2.optimize()
+        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
+
+        # Ensure revert_to_reversible is robust to solutions generated both
+        # before and after reversibility conversion, or not solved at all.
+        model3 = create_test_model("textbook")
+        model3.optimize()
+        convert_to_irreversible(model3)
+        revert_to_reversible(model3)
+        self.assertAlmostEqual(model1.solution.f, model3.solution.f, places=3)
+
+        # test reaction where both bounds are negative
+        model4 = create_test_model("textbook")
+        glc = model4.reactions.get_by_id("EX_glc__D_e")
+        glc.upper_bound = -1
+        convert_to_irreversible(model4)
+        model4.optimize()
+        self.assertAlmostEqual(model1.solution.f, model4.solution.f, places=3)
+        glc_rev = model4.reactions.get_by_id(glc.notes["reflection"])
+        self.assertEqual(glc_rev.lower_bound, 1)
+        self.assertEqual(glc.upper_bound, 0)
+        revert_to_reversible(model4)
+        self.assertEqual(glc.upper_bound, -1)
+
+    def test_escape_ids(self):
+        model = create_test_model('textbook')
+        model.reactions.PGI.gene_reaction_rule = "a.b or c"
+        self.assertIn("a.b", model.genes)
+        escape_ID(model)
+        self.assertNotIn("a.b", model.genes)
+
+    def test_gene_knockout_computation(self):
+        cobra_model = create_test_model()
+
+        def find_gene_knockout_reactions_fast(cobra_model, gene_list):
+            compiled_rules = get_compiled_gene_reaction_rules(
+                cobra_model)
+            return find_gene_knockout_reactions(
+                cobra_model, gene_list,
+                compiled_gene_reaction_rules=compiled_rules)
+
+        def get_removed(m):
+            return {x.id for x in m._trimmed_reactions}
+
+        def test_computation(m, gene_ids, expected_reaction_ids):
+            genes = [m.genes.get_by_id(i) for i in gene_ids]
+            expected_reactions = {m.reactions.get_by_id(i)
+                                  for i in expected_reaction_ids}
+            removed1 = set(find_gene_knockout_reactions(m, genes))
+            removed2 = set(find_gene_knockout_reactions_fast(m, genes))
+            self.assertEqual(removed1, expected_reactions)
+            self.assertEqual(removed2, expected_reactions)
+            delete_model_genes(m, gene_ids, cumulative_deletions=False)
+            self.assertEqual(get_removed(m), expected_reaction_ids)
+            undelete_model_genes(m)
+
+        gene_list = ['STM1067', 'STM0227']
+        dependent_reactions = {'3HAD121', '3HAD160', '3HAD80', '3HAD140',
+                               '3HAD180', '3HAD100', '3HAD181', '3HAD120',
+                               '3HAD60', '3HAD141', '3HAD161', 'T2DECAI',
+                               '3HAD40'}
+        test_computation(cobra_model, gene_list, dependent_reactions)
+        test_computation(cobra_model, ['STM4221'], {'PGI'})
+        test_computation(cobra_model, ['STM1746.S'], {'4PEPTabcpp'})
+        # test cumulative behavior
+        delete_model_genes(cobra_model, gene_list[:1])
+        delete_model_genes(cobra_model, gene_list[1:],
+                           cumulative_deletions=True)
+        delete_model_genes(cobra_model, ["STM4221"],
+                           cumulative_deletions=True)
+        dependent_reactions.add('PGI')
+        self.assertEqual(get_removed(cobra_model), dependent_reactions)
+        # non-cumulative following cumulative
+        delete_model_genes(cobra_model, ["STM4221"],
+                           cumulative_deletions=False)
+        self.assertEqual(get_removed(cobra_model), {'PGI'})
+        # make sure on reset that the bounds are correct
+        reset_bound = cobra_model.reactions.get_by_id("T2DECAI").upper_bound
+        self.assertEqual(reset_bound, 1000.)
+        # test computation when gene name is a subset of another
+        test_model = Model()
+        test_reaction_1 = Reaction("test1")
+        test_reaction_1.gene_reaction_rule = "eggs or (spam and eggspam)"
+        test_model.add_reaction(test_reaction_1)
+        test_computation(test_model, ["eggs"], set())
+        test_computation(test_model, ["eggs", "spam"], {'test1'})
+        # test computation with nested boolean expression
+        test_reaction_1.gene_reaction_rule = \
+            "g1 and g2 and (g3 or g4 or (g5 and g6))"
+        test_computation(test_model, ["g3"], set())
+        test_computation(test_model, ["g1"], {'test1'})
+        test_computation(test_model, ["g5"], set())
+        test_computation(test_model, ["g3", "g4", "g5"], {'test1'})
+        # test computation when gene names are python expressions
+        test_reaction_1.gene_reaction_rule = "g1 and (for or in)"
+        test_computation(test_model, ["for", "in"], {'test1'})
+        test_computation(test_model, ["for"], set())
+        test_reaction_1.gene_reaction_rule = "g1 and g2 and g2.conjugate"
+        test_computation(test_model, ["g2"], {"test1"})
+        test_computation(test_model, ["g2.conjugate"], {"test1"})
+        test_reaction_1.gene_reaction_rule = "g1 and (try:' or 'except:1)"
+        test_computation(test_model, ["try:'"], set())
+        test_computation(test_model, ["try:'", "'except:1"], {"test1"})
+
+    def test_remove_genes(self):
+        m = Model("test")
+        m.add_reactions([Reaction("r" + str(i + 1)) for i in range(8)])
+        self.assertEqual(len(m.reactions), 8)
+        rxns = m.reactions
+        rxns.r1.gene_reaction_rule = "(a and b) or (c and a)"
+        rxns.r2.gene_reaction_rule = "(a and b and d and e)"
+        rxns.r3.gene_reaction_rule = "(a and b) or (b and c)"
+        rxns.r4.gene_reaction_rule = "(f and b) or (b and c)"
+        rxns.r5.gene_reaction_rule = "x"
+        rxns.r6.gene_reaction_rule = "y"
+        rxns.r7.gene_reaction_rule = "x or     z"
+        rxns.r8.gene_reaction_rule = ""
+        self.assertIn("a", m.genes)
+        self.assertIn("x", m.genes)
+        remove_genes(m, ["a"], remove_reactions=False)
+        self.assertNotIn("a", m.genes)
+        self.assertIn("x", m.genes)
+        self.assertEqual(rxns.r1.gene_reaction_rule, "")
+        self.assertEqual(rxns.r2.gene_reaction_rule, "")
+        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
+        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
+        self.assertEqual(rxns.r5.gene_reaction_rule, "x")
+        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
+        self.assertEqual(rxns.r7.genes, {m.genes.x, m.genes.z})
+        self.assertEqual(rxns.r8.gene_reaction_rule, "")
+        remove_genes(m, ["x"], remove_reactions=True)
+        self.assertEqual(len(m.reactions), 7)
+        self.assertNotIn("r5", m.reactions)
+        self.assertNotIn("x", m.genes)
+        self.assertEqual(rxns.r1.gene_reaction_rule, "")
+        self.assertEqual(rxns.r2.gene_reaction_rule, "")
+        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
+        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
+        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
+        self.assertEqual(rxns.r7.gene_reaction_rule, "z")
+        self.assertEqual(rxns.r7.genes, {m.genes.z})
+        self.assertEqual(rxns.r8.gene_reaction_rule, "")
+
+    def test_SBO_annotation(self):
+        model = create_test_model("textbook")
+        rxns = model.reactions
+        rxns.EX_o2_e.annotation.clear()
+        fake_DM = Reaction("DM_h_c")
+        model.add_reaction(fake_DM)
+        fake_DM.add_metabolites({model.metabolites.get_by_id("h_c"): -1})
+        # this exchange will be set wrong. The function should not overwrite
+        # an existing SBO annotation
+        rxns.get_by_id("EX_h_e").annotation["SBO"] = "SBO:0000628"
+        add_SBO(model)
+        self.assertEqual(rxns.EX_o2_e.annotation["SBO"], "SBO:0000627")
+        self.assertEqual(rxns.DM_h_c.annotation["SBO"], "SBO:0000628")
+        self.assertEqual(rxns.EX_h_e.annotation["SBO"], "SBO:0000628")
+
+    def test_validate_reaction_bounds(self):
+        model = create_test_model("textbook")
+        model.reactions[0].lower_bound = float("-inf")
+        model.reactions[1].lower_bound = float("nan")
+        model.reactions[0].upper_bound = float("inf")
+        model.reactions[1].upper_bound = float("nan")
+        errors = check_reaction_bounds(model)
+        self.assertEqual(len(errors), 4)
+
+    def test_validate_formula_compartment(self):
+        model = create_test_model("textbook")
+        model.metabolites[1].compartment = "fake"
+        model.metabolites[1].formula = "(a*.bcde)"
+        errors = check_metabolite_compartment_formula(model)
+        self.assertEqual(len(errors), 2)
+
+    def test_validate_mass_balance(self):
+        model = create_test_model("textbook")
+        self.assertEqual(len(check_mass_balance(model)), 0)
+        # if we remove the SBO term which marks the reaction as
+        # mass balanced, then the reaction should be detected as
+        # no longer mass balanced
+        EX_rxn = model.reactions.query("EX")[0]
+        EX_rxn.annotation.pop("SBO")
+        balance = check_mass_balance(model)
+        self.assertEqual(len(balance), 1)
+        self.assertIn(EX_rxn, balance)
+
+
+# make a test suite to run all of the tests
+loader = TestLoader()
+suite = loader.loadTestsFromModule(sys.modules[__name__])
+
+
+def test_all():
+    TextTestRunner(verbosity=2).run(suite)
+
+if __name__ == "__main__":
+    test_all()
diff --git a/cobra/test/unit_tests.py b/cobra/test/unit_tests.py
index 61cdc0a..1edc147 100644
--- a/cobra/test/unit_tests.py
+++ b/cobra/test/unit_tests.py
@@ -3,6 +3,7 @@ from unittest import TestCase, TestLoader, TextTestRunner, skipIf
 from copy import copy, deepcopy
 from pickle import loads, dumps, HIGHEST_PROTOCOL
 import warnings
+import re
 
 if __name__ == "__main__":
     sys.path.insert(0, "../..")
@@ -168,12 +169,24 @@ class TestDictList(TestCase):
 
     def testQuery(self):
         obj2 = Object("test2")
+        obj2.name = "foobar1"
         self.list.append(obj2)
         result = self.list.query("test1")  # matches only test1
         self.assertEqual(len(result), 1)
         self.assertEqual(result[0], self.obj)
+        result = self.list.query("foo", "name")  # matches only test2
+        self.assertEqual(len(result), 1)
+        self.assertEqual(result[0], obj2)
         result = self.list.query("test")  # matches test1 and test2
         self.assertEqual(len(result), 2)
+        # test with a regular expression
+        result = self.list.query(re.compile("test[0-9]"))
+        self.assertEqual(len(result), 2)
+        result = self.list.query(re.compile("test[29]"))
+        self.assertEqual(len(result), 1)
+        # test query of name
+        result = self.list.query(re.compile("foobar."), "name")
+        self.assertEqual(len(result), 1)
 
     def testRemoval(self):
         obj_list = DictList(Object("test%d" % (i)) for i in range(2, 10))
@@ -211,6 +224,29 @@ class TestDictList(TestCase):
         self.assertRaises(ValueError, obj_list.__setitem__, slice(5, 7),
                           [Object("testd"), Object("testd")])
 
+    def testSortandReverse(self):
+        dl = DictList(Object("test%d" % (i)) for i in reversed(range(10)))
+        self.assertEqual(dl[0].id, "test9")
+        dl.sort()
+        self.assertEqual(len(dl), 10)
+        self.assertEqual(dl[0].id, "test0")
+        self.assertEqual(dl.index("test0"), 0)
+        dl.reverse()
+        self.assertEqual(dl[0].id, "test9")
+        self.assertEqual(dl.index("test0"), 9)
+
+    def testDir(self):
+        """makes sure tab complete will work"""
+        attrs = dir(self.list)
+        self.assertIn("test1", attrs)
+        self.assertIn("_dict", attrs)  # attribute of DictList
+
+    def testUnion(self):
+        self.list.union([Object("test1"), Object("test2")])
+        # should only add 1 element
+        self.assertEqual(len(self.list), 2)
+        self.assertEqual(self.list.index("test2"), 1)
+
 
 class CobraTestCase(TestCase):
     def setUp(self):
diff --git a/cobra/test_all.py b/cobra/test_all.py
deleted file mode 100755
index 2d7ee0b..0000000
--- a/cobra/test_all.py
+++ /dev/null
@@ -1,7 +0,0 @@
-#!/usr/bin/env python
-assert __name__ == "__main__"
-from os.path import abspath, split, join
-import sys
-sys.path.insert(0, abspath(join(split(__file__)[0], "..")))
-from cobra.test import test_all
-test_all()
diff --git a/cobra/topology/reporter_metabolites.py b/cobra/topology/reporter_metabolites.py
index e8f928e..84bd6f2 100644
--- a/cobra/topology/reporter_metabolites.py
+++ b/cobra/topology/reporter_metabolites.py
@@ -1,10 +1,12 @@
 #cobra.topology.reporter_metabolites.py: Module for topological analysis of cobra_models
 #Based on Patil et al 2005 PNAS 102:2685-9
 #TODO: Validate cobra.core compliance
+from __future__ import print_function
 from copy import deepcopy
 from numpy import array, corrcoef, mean, std, tril, where, unique, zeros
 from scipy.stats import norm, randint
 from collections import defaultdict
+from six import iteritems
 
 def identify_reporter_metabolites(cobra_model, reaction_scores_dict,
                                   number_of_randomizations=1000, number_of_layers=1,
@@ -109,7 +111,7 @@ def identify_reporter_metabolites(cobra_model, reaction_scores_dict,
             correction_dict[i] = [mean(random_score_distribution),
                                       std(random_score_distribution,ddof=1)] 
 
-    for the_metabolite, the_score in metabolite_scores.iteritems():
+    for the_metabolite, the_score in iteritems(metabolite_scores):
         number_of_connections = metabolite_connections[the_metabolite]
         if number_of_connections > 0:
             #Correct based on background distribution
diff --git a/documentation_builder/io.ipynb b/documentation_builder/io.ipynb
index e305b68..295bde1 100644
--- a/documentation_builder/io.ipynb
+++ b/documentation_builder/io.ipynb
@@ -65,7 +65,7 @@
     {
      "data": {
       "text/plain": [
-       "<Model mini_textbook at 0x7f35fc0de3d0>"
+       "<Model mini_textbook at 0x7f246d4e2e50>"
       ]
      },
      "execution_count": 2,
@@ -107,7 +107,7 @@
     {
      "data": {
       "text/plain": [
-       "<Model mini_textbook at 0x7f35c46c7990>"
+       "<Model mini_textbook at 0x7f2436c65a10>"
       ]
      },
      "execution_count": 4,
@@ -149,7 +149,7 @@
     {
      "data": {
       "text/plain": [
-       "<Model mini_textbook at 0x7f35c46dcd50>"
+       "<Model mini_textbook at 0x7f2436c7c850>"
       ]
      },
      "execution_count": 6,
@@ -169,7 +169,7 @@
    },
    "outputs": [],
    "source": [
-    "cobra.io.write_sbml_model(textbook_model, \"test.json\")"
+    "cobra.io.save_json_model(textbook_model, \"test.json\")"
    ]
   },
   {
@@ -198,7 +198,7 @@
     {
      "data": {
       "text/plain": [
-       "<Model mini_textbook at 0x7f35c46dcc90>"
+       "<Model mini_textbook at 0x7f2436c7c810>"
       ]
      },
      "execution_count": 8,
@@ -228,7 +228,7 @@
     {
      "data": {
       "text/plain": [
-       "<Model mini_textbook at 0x7f35c46dcb90>"
+       "<Model mini_textbook at 0x7f2436c65510>"
       ]
      },
      "execution_count": 9,
diff --git a/documentation_builder/io.rst b/documentation_builder/io.rst
index e2fb878..e9dd914 100644
--- a/documentation_builder/io.rst
+++ b/documentation_builder/io.rst
@@ -3,8 +3,9 @@ Reading and Writing Models
 ==========================
 
 Cobrapy supports reading and writing models in SBML (with and without
-FBC), JSON, MAT, and pickle formats. Generally, SBML with FBC is the
-preferred format for general use.
+FBC), JSON, MAT, and pickle formats. Generally, SBML with FBC version 2
+is the preferred format for general use. The JSON format may be more
+useful for cobrapy-specific functionality.
 
 The package also ships with test models in various formats for testing
 purposes.
@@ -25,7 +26,7 @@ purposes.
 .. parsed-literal::
 
     mini test files: 
-    mini.mat, mini_cobra.xml, mini.json, mini_fbc2.xml, mini_fbc1.xml, mini.pickle
+    mini.mat, mini_cobra.xml, mini.json, mini_fbc2.xml.gz, mini_fbc2.xml.bz2, mini_fbc2.xml, mini_fbc1.xml, mini.pickle
 
 
 SBML
@@ -50,7 +51,7 @@ requirements in order to generate a valid SBML file.
 
 .. parsed-literal::
 
-    <Model mini_textbook at 0x7f82fc4aec90>
+    <Model mini_textbook at 0x7f246d4e2e50>
 
 
 
@@ -77,7 +78,7 @@ flag can be used can be used.
 
 .. parsed-literal::
 
-    <Model mini_textbook at 0x7f82cc744c10>
+    <Model mini_textbook at 0x7f2436c65a10>
 
 
 
@@ -90,8 +91,7 @@ JSON
 
 cobrapy models have a `JSON <https://en.wikipedia.org/wiki/JSON>`__
 (JavaScript Object Notation) representation. This format was crated for
-interoperability with `escher <https://escher.github.io>`__. Additional
-fields, however, will not be saved.
+interoperability with `escher <https://escher.github.io>`__.
 
 .. code:: python
 
@@ -102,13 +102,13 @@ fields, however, will not be saved.
 
 .. parsed-literal::
 
-    <Model mini_textbook at 0x7f82fc4d4f50>
+    <Model mini_textbook at 0x7f2436c7c850>
 
 
 
 .. code:: python
 
-    cobra.io.write_sbml_model(textbook_model, "test.json")
+    cobra.io.save_json_model(textbook_model, "test.json")
 
 MATLAB
 ------
@@ -133,7 +133,7 @@ reading function:
 
 .. parsed-literal::
 
-    <Model mini_textbook at 0x7f82cc97b110>
+    <Model mini_textbook at 0x7f2436c7c810>
 
 
 
@@ -149,7 +149,7 @@ variable to read from, and the variable\_name paramter is unnecessary.
 
 .. parsed-literal::
 
-    <Model mini_textbook at 0x7f82fc4d4450>
+    <Model mini_textbook at 0x7f2436c65510>
 
 
 
@@ -163,19 +163,8 @@ Pickle
 ------
 
 Cobra models can be serialized using the python serialization format,
-`pickle <https://docs.python.org/2/library/pickle.html>`__. While this
-will save any extra fields which may have been created, it does not work
-with any other tools and can break between cobrapy major versions. JSON,
-SBML, and MAT are generally the preferred format.
+`pickle <https://docs.python.org/2/library/pickle.html>`__.
 
-.. code:: python
-
-    from pickle import load, dump
-    
-    # read in the test models
-    with open(os.path.join(cobra.test.data_directory, "mini.pickle"), "rb") as infile:
-        mini_model = load(infile)
-    
-    # output to a file
-    with open("test.pickle", "wb") as outfile:
-        dump(textbook_model, outfile)
+Please note that use of the pickle format is generally not recommended
+for most use cases. JSON, SBML, and MAT are generally the preferred
+formats.
diff --git a/documentation_builder/loopless.ipynb b/documentation_builder/loopless.ipynb
index 7015d42..1c8da93 100644
--- a/documentation_builder/loopless.ipynb
+++ b/documentation_builder/loopless.ipynb
@@ -18,8 +18,8 @@
     "$$ 0 \\le lb \\le v \\le ub \\le \\max(ub) $$\n",
     "\n",
     "We will add in boolean indicators as well, such that\n",
-    "$$ \\max(ub) i \\ge v $$\n",
-    "$$ i \\in {0, 1} $$\n",
+    "$$ \\max(ub) \\cdot i \\ge v $$\n",
+    "$$ i \\in \\{0, 1\\} $$\n",
     "\n",
     "We also want to ensure that an entry in the row space of S also exists with negative values wherever v is nonzero. In this expression, $1-i$ acts as a not to indicate inactivity of a reaction.\n",
     "\n",
diff --git a/documentation_builder/loopless.rst b/documentation_builder/loopless.rst
index 38f5dcc..ef613cc 100644
--- a/documentation_builder/loopless.rst
+++ b/documentation_builder/loopless.rst
@@ -26,9 +26,9 @@ We can make the model irreversible, so that all reactions will satisfy
 
 We will add in boolean indicators as well, such that
 
-.. math::  \max(ub) i \ge v 
+.. math::  \max(ub) \cdot i \ge v 
 
-.. math::  i \in {0, 1} 
+.. math::  i \in \{0, 1\} 
 
 We also want to ensure that an entry in the row space of S also exists
 with negative values wherever v is nonzero. In this expression,
diff --git a/documentation_builder/phenotype_phase_plane.ipynb b/documentation_builder/phenotype_phase_plane.ipynb
index c5eb486..77d39a9 100644
--- a/documentation_builder/phenotype_phase_plane.ipynb
+++ b/documentation_builder/phenotype_phase_plane.ipynb
@@ -22,7 +22,6 @@
     "%matplotlib inline\n",
     "from time import time\n",
     "\n",
-    "\n",
     "import cobra.test\n",
     "from cobra.flux_analysis import calculate_phenotype_phase_plane\n",
     "\n",
@@ -45,9 +44,9 @@
    "outputs": [
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXeYHOd15vv7KnVO0z0JGAQCnAGYwQzmnKNkURQVKK1s\nS6Ys2V7JkmV7bTquvbvXe6/WsnVtK1kiRZGEmASRFJMokmIAEwiQAJEBTujpiZ1Dxf3jq+YMIWQM\nBqDY7/PMg5mq6kJ9PT2n3jrnPe8RnufRQgsttNDC7EA53BfQQgsttPBBQivottBCCy3MIlpBt4UW\nWmhhFtEKui200EILs4hW0G2hhRZamEW0gm4LLbTQwiyiFXRbaKGFFmYRraDbQgsttDCLaAXdFlpo\noYVZRCvottBCCy3MIlpBt4UWWmhhFtEKui200EILs4hW0G2hhRZamEW0gm4LLbTQwiyiFXRbaKGF\nFmYR [...]
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXewHHV69/vpNNOT8wk6ygEhIVgEkhBBBIksYFmWDez6\nXe/i9etQ1/dulf2vq9Z/2eWyb9W133rtd9fZC7sLLDlnJKIIEghEkFA8acKZHDveP1o9jA7nSEfS\nCRLMp4oq0WdCd0/PM99+fs/zfQTbtunSpUuXLrODONc70KVLly5fJ7pBt0uXLl1mkW7Q7dKlS5dZ\npBt0u3Tp0mUW6QbdLl26dJlFukG3S5cuXWYR+QR/79aTdenSpcvJI0z2h67S7dKlS5dZpBt0u3Tp\n0mUW6QbdLl26dJlFukG3S5cuXWaRbtDt0qVLl1mkG3S7dOnSZRbpBt0uXbp0mUW6QbdLly5dZpFu\n0O3S [...]
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x7fe334311210>"
+       "<matplotlib.figure.Figure at 0x7fd5f40df810>"
       ]
      },
      "metadata": {},
@@ -63,7 +62,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If [brewer2mpl](https://pypi.python.org/pypi/brewer2mpl/) is installed, other color schemes can be used as well"
+    "If [palettable](https://github.com/jiffyclub/palettable) is installed, other color schemes can be used as well"
    ]
   },
   {
@@ -75,9 +74,9 @@
    "outputs": [
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdYU4f6xz/nnAwghA2CouAA3HuPaq3XrbV111FXq7Va\nW7un3eN6a+u2w2rVWq2tVWsdddQ9cG8FFREwKIIIJIQkJ+f3x0ku1OsW0P6az/PwtHKS4zkY3nzz\nju8rKIqCBw8ePHgoHcT7fQEePHjw8E/CE3Q9ePDgoRTxBF0PHjx4KEU8QdeDBw8eShFP0PXgwYOH\nUsQTdD148OChFPEEXQ8ePHgoRTxB14MHDx5KEU/Q9eDBg4dSxBN0PXjw4KEU8QRdDx48eChFPEHX\ngwcPHkoRT9D14MGDh1LEE3Q9ePDgoRTxBF0PHjx4KEU8QdeDBw8eShFP0PXgwYOHUsQTdD148OCh\nFNHc [...]
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWdYlHfWxn9TmBkGZuhNUewNFbuosfdewN5LLLGlbt5s\nNtlkd7MxMcbYW4zG3nuNLSo2bKBiQwVUFAWp02Da+4E8syMBRUXQ3fldl19kGOYpc577f/7n3Edk\ntVpx4MCBAwfFg7ikP4ADBw4c/C/hCLoOHDhwUIw4gq4DBw4cFCOOoOvAgQMHxYgj6Dpw4MBBMeII\nug4cOHBQjEif83NHPZkDBw4cvDiign7gULoOHDhwUIw4gq4DBw4cFCOOoOvAgQMHxYgj6Dpw4MBB\nMeIIug4cOHBQjDiCrgMHDhwUI46g68CBAwfFiCPoOnDgwEEx4gi6Dhw4cFCMOIKuAwcOHBQjjqDr\nwIED [...]
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x7fe301c440d0>"
+       "<matplotlib.figure.Figure at 0x7fd5c1a56e10>"
       ]
      },
      "metadata": {},
@@ -85,9 +84,9 @@
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXecFPX5x9/faduvH9zR69GbiHSwF6wxscaoscQSjYlJ\nNOUX0ajxpzGJJYm/mJjEkqgxUWMDpIrYQAHp3AEHHHBw/W77zs7M74/vrHcSOseBcd+v1730dnaH\nmb3dZz7zfJ/n8wjHcciSJUuWLB2DcrQPIEuWLFm+TGSDbpYsWbJ0INmgmyVLliwdSDboZsmSJUsH\nkg26WbJkydKBZINulixZsnQg2aCbJUuWLB1INuhmyZIlSweSDbpZsmTJ0oFkg26WLFmydCDZoJsl\nS5YsHUg26GbJkiVLB5INulmyZMnSgWSDbpYsWbJ0INmgmyVLliwdSDboZsmSJUsHkg26WbJkydKB\nZINu [...]
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWd4HPW1h98p27tW1bLce+8d29iYahtCh9ADISEhkISb\n3Nw0ICQEAgklIQTSICEQQiDgjnvvvfdepJW02t5n7gd5Nmsh2bItSzbM+zz+4NXu7M7s7JnfnP85\nvyOoqoqOjo6OTvMgtvQH0NHR0fkioQddHR0dnWZED7o6Ojo6zYgedHV0dHSaET3o6ujo6DQjetDV\n0dHRaUbkM/xdryfT0dHROXuEhv6gK10dHR2dZkQPujo6OjrNiB50dXR0dJoRPejq6OjoNCN60NXR\n0dFpRvSgq6Ojo9OM6EFXR0dHpxnRg66Ojo5OM6IHXR0dHZ1mRA+6Ojo6Os2IHnR1dHR0mpEzeS/o\n6Fz0 [...]
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x7fe301e8f650>"
+       "<matplotlib.figure.Figure at 0x7fd5c1737390>"
       ]
      },
      "metadata": {},
@@ -115,9 +114,9 @@
    "outputs": [
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXeYHNWV/v+5VZ3D5KgMQoEkBBIiCBA5KCDA2BgHbGOb\ntb1417v79dfrdeC3Oa+/Xq/Xuw5gY2NsAzYgiYxBIJIAIQRC0owEitMTemLn7qq6vz9u1UzPTPcE\naYIw/T6PHlBX3epbrZlTp8857/sKKSUllFBCCSVMDbTp3kAJJZRQwgcJpaBbQgkllDCFKAXdEkoo\noYQpRCnollBCCSVMIUpBt4QSSihhClEKuiWUUEIJU4hS0C2hhBJKmEKUgm4JJZRQwhSiFHRLKKGE\nEqYQpaBbQgkllDCFKAXdEkoooYQpRCnollBCCSVMIUpBt4QSSihhClEKuiWUUEIJU4hS0C2hhBJK\nmEKU [...]
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvWmUHPWV5v2L3NeqzMrMWlRakIQQCEloBwHa0AJaSqy2\ncdtAQ4OBadvted/u4zPvjH1m3N0zbY9nusftHnvMoe22jRmMDWiXEAiEFrSgBUksQmhXrZlZue+Z\nEe+HqMjKvUql2mTyOUcfVJn/yIjIyBs37r3P8wiSJFFFFVVUUcXwQDXSO1BFFVVU8UVCNehWUUUV\nVQwjqkG3iiqqqGIYUQ26VVRRRRXDiGrQraKKKqoYRlSDbhVVVFHFMELTx+vVebIqqqiiiquHUO6F\naqZbRRVVVDGMqAbdKqqooophRDXoVlFFFVUMI6pBt4oqqqhiGFENulVUUUUVw4hq0K2iiiqqGEZU\ng24V [...]
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x7fe301e44910>"
+       "<matplotlib.figure.Figure at 0x7fd5c163ff50>"
       ]
      },
      "metadata": {},
@@ -148,8 +147,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "took 0.40 seconds with 1 process\n",
-      "took 0.21 seconds with 4 process\n"
+      "took 0.41 seconds with 1 process\n",
+      "took 0.29 seconds with 4 process\n"
      ]
     }
    ],
diff --git a/documentation_builder/phenotype_phase_plane.rst b/documentation_builder/phenotype_phase_plane.rst
index ce7af26..783279d 100644
--- a/documentation_builder/phenotype_phase_plane.rst
+++ b/documentation_builder/phenotype_phase_plane.rst
@@ -15,7 +15,6 @@ Here, we will make one for the "textbook" *E. coli* core model.
     %matplotlib inline
     from time import time
     
-    
     import cobra.test
     from cobra.flux_analysis import calculate_phenotype_phase_plane
     
@@ -34,7 +33,7 @@ and Oxygen.
 .. image:: phenotype_phase_plane_files/phenotype_phase_plane_3_0.png
 
 
-If `brewer2mpl <https://pypi.python.org/pypi/brewer2mpl/>`__ is
+If `palettable <https://github.com/jiffyclub/palettable>`__ is
 installed, other color schemes can be used as well
 
 .. code:: python
@@ -81,6 +80,6 @@ The code can also use multiple processes to speed up calculations
 
 .. parsed-literal::
 
-    took 0.40 seconds with 1 process
-    took 0.21 seconds with 4 process
+    took 0.41 seconds with 1 process
+    took 0.29 seconds with 4 process
 
diff --git a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_3_0.png b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_3_0.png
index b3f96fd..bc6dac4 100644
Binary files a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_3_0.png and b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_3_0.png differ
diff --git a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_0.png b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_0.png
index 46c5994..8413a31 100644
Binary files a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_0.png and b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_0.png differ
diff --git a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_1.png b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_1.png
index 0c2c6c4..ff8c9d4 100644
Binary files a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_1.png and b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_5_1.png differ
diff --git a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_7_0.png b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_7_0.png
index a49bce8..47ce17c 100644
Binary files a/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_7_0.png and b/documentation_builder/phenotype_phase_plane_files/phenotype_phase_plane_7_0.png differ
diff --git a/documentation_builder/simulating.ipynb b/documentation_builder/simulating.ipynb
index d83a1c1..00c6b54 100644
--- a/documentation_builder/simulating.ipynb
+++ b/documentation_builder/simulating.ipynb
@@ -11,7 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 1,
    "metadata": {
     "collapsed": false
    },
@@ -41,7 +41,7 @@
     {
      "data": {
       "text/plain": [
-       "<Solution 0.87 at 0x7fb3f849ca10>"
+       "<Solution 0.87 at 0x7fe558058b50>"
       ]
      },
      "execution_count": 2,
@@ -137,7 +137,7 @@
     {
      "data": {
       "text/plain": [
-       "{<Reaction Biomass_Ecoli_core at 0x7fb3c899be90>: 1.0}"
+       "{<Reaction Biomass_Ecoli_core at 0x7fe526516490>: 1.0}"
       ]
      },
      "execution_count": 5,
@@ -166,7 +166,7 @@
     {
      "data": {
       "text/plain": [
-       "{<Reaction ATPM at 0x7fb3c899bbd0>: 1}"
+       "{<Reaction ATPM at 0x7fe526516210>: 1}"
       ]
      },
      "execution_count": 6,
@@ -192,7 +192,7 @@
     {
      "data": {
       "text/plain": [
-       "<Solution 175.00 at 0x7fb3c895de50>"
+       "174.99999999999997"
       ]
      },
      "execution_count": 7,
@@ -201,7 +201,7 @@
     }
    ],
    "source": [
-    "model.optimize()"
+    "model.optimize().f"
    ]
   },
   {
@@ -213,7 +213,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 8,
    "metadata": {
     "collapsed": false
    },
@@ -221,10 +221,10 @@
     {
      "data": {
       "text/plain": [
-       "{<Reaction Biomass_Ecoli_core at 0x7fb3c899be90>: 1.0}"
+       "{<Reaction Biomass_Ecoli_core at 0x7fe526516490>: 1.0}"
       ]
      },
-     "execution_count": 9,
+     "execution_count": 8,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -246,7 +246,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -392,7 +392,7 @@
        "EX_ac_e            -3.933509e-15  0.000000e+00"
       ]
      },
-     "execution_count": 35,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -411,7 +411,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 10,
    "metadata": {
     "collapsed": false
    },
@@ -481,7 +481,7 @@
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ATPM</th>\n",
-       "      <td>8.390000e+00</td>\n",
+       "      <td>2.555100e+01</td>\n",
        "      <td>8.390000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
@@ -545,7 +545,7 @@
        "AKGDH               8.045934e+00   0.000000\n",
        "AKGt2r              0.000000e+00  -1.430083\n",
        "ALCD2x              0.000000e+00  -2.214323\n",
-       "ATPM                8.390000e+00   8.390000\n",
+       "ATPM                2.555100e+01   8.390000\n",
        "ATPS4r              5.938106e+01  34.825618\n",
        "Biomass_Ecoli_core  8.739215e-01   0.786529\n",
        "CO2t               -1.520653e+01 -26.528850\n",
@@ -557,7 +557,7 @@
        "EX_ac_e             3.813556e+00   0.000000"
       ]
      },
-     "execution_count": 36,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -566,6 +566,56 @@
     "fva_result = cobra.flux_analysis.flux_variability_analysis(model, model.reactions[:20], fraction_of_optimum=0.9)\n",
     "pandas.DataFrame.from_dict(fva_result).T"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Running pFBA\n",
+    "\n",
+    "Parsimonious FBA (often written pFBA) finds a flux distribution which gives the optimal growth rate, but minimizes the total sum of flux. This involves solving two sequential linear programs, but is handled transparently by cobrapy. For more details on pFBA, please see [Lewis et al. (2010)](http://dx.doi.org/10.1038/msb.2010.47)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "FBA_solution = model.optimize()\n",
+    "pFBA_solution = cobra.flux_analysis.optimize_minimal_flux(model)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "These functions should give approximately the same objective value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1.1102230246251565e-16"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "abs(FBA_solution.f - pFBA_solution.f)"
+   ]
   }
  ],
  "metadata": {
diff --git a/documentation_builder/simulating.rst b/documentation_builder/simulating.rst
index 6dd67bc..d007044 100644
--- a/documentation_builder/simulating.rst
+++ b/documentation_builder/simulating.rst
@@ -26,7 +26,7 @@ Running FBA
 
 .. parsed-literal::
 
-    <Solution 0.87 at 0x7fb3f849ca10>
+    <Solution 0.87 at 0x7fe558058b50>
 
 
 
@@ -88,7 +88,7 @@ only one objective reaction, with an objective coefficient of 1.
 
 .. parsed-literal::
 
-    {<Reaction Biomass_Ecoli_core at 0x7fb3c899be90>: 1.0}
+    {<Reaction Biomass_Ecoli_core at 0x7fe526516490>: 1.0}
 
 
 
@@ -109,20 +109,20 @@ which can be a reaction object (or just it's name), or a dict of
 
 .. parsed-literal::
 
-    {<Reaction ATPM at 0x7fb3c899bbd0>: 1}
+    {<Reaction ATPM at 0x7fe526516210>: 1}
 
 
 
 .. code:: python
 
-    model.optimize()
+    model.optimize().f
 
 
 
 
 .. parsed-literal::
 
-    <Solution 175.00 at 0x7fb3c895de50>
+    174.99999999999997
 
 
 
@@ -140,7 +140,7 @@ Reaction.objective\_coefficient directly.
 
 .. parsed-literal::
 
-    {<Reaction Biomass_Ecoli_core at 0x7fb3c899be90>: 1.0}
+    {<Reaction Biomass_Ecoli_core at 0x7fe526516490>: 1.0}
 
 
 
@@ -352,7 +352,7 @@ for reactions at 90% optimality.
         </tr>
         <tr>
           <th>ATPM</th>
-          <td>8.390000e+00</td>
+          <td>2.555100e+01</td>
           <td>8.390000</td>
         </tr>
         <tr>
@@ -405,3 +405,32 @@ for reactions at 90% optimality.
     </div>
 
 
+
+Running pFBA
+------------
+
+Parsimonious FBA (often written pFBA) finds a flux distribution which
+gives the optimal growth rate, but minimizes the total sum of flux. This
+involves solving two sequential linear programs, but is handled
+transparently by cobrapy. For more details on pFBA, please see `Lewis et
+al. (2010) <http://dx.doi.org/10.1038/msb.2010.47>`__.
+
+.. code:: python
+
+    FBA_solution = model.optimize()
+    pFBA_solution = cobra.flux_analysis.optimize_minimal_flux(model)
+
+These functions should give approximately the same objective value
+
+.. code:: python
+
+    abs(FBA_solution.f - pFBA_solution.f)
+
+
+
+
+.. parsed-literal::
+
+    1.1102230246251565e-16
+
+
diff --git a/documentation_builder/solvers.ipynb b/documentation_builder/solvers.ipynb
index 86c355f..e39688c 100644
--- a/documentation_builder/solvers.ipynb
+++ b/documentation_builder/solvers.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#Solver Interface\n",
+    "# Solver Interface\n",
     "\n",
     "Each cobrapy solver must expose the following API. The solvers all will have their own distinct LP object types, but each can be manipulated by these functions. This API can be used directly when implementing algorithms efficiently on linear programs because it has 2 primary benefits:\n",
     "\n",
diff --git a/setup.py b/setup.py
index 4282ba0..23bb89b 100644
--- a/setup.py
+++ b/setup.py
@@ -102,15 +102,19 @@ try:
         library_dirs.append(abspath("."))
     if isfile("glpk.h"):
         include_dirs.append(abspath("."))
-    if name == "posix":
+    # if the glpk files are not in the current directory attempt to
+    # auto-detect their location by finding the location of the glpsol
+    # command
+    if name == "posix" and len(include_dirs) == 0 and len(library_dirs) == 0:
         from subprocess import check_output
         try:
-            glpksol_path = check_output(["which", "glpsol"]).strip()
+            glpksol_path = check_output(["which", "glpsol"],
+                                        universal_newlines=True).strip()
             glpk_path = abspath(join(dirname(glpksol_path), ".."))
             include_dirs.append(join(glpk_path, "include"))
             library_dirs.append(join(glpk_path, "lib"))
-        except:
-            None
+        except Exception as e:
+            print('Could not autodetect include and library dirs: ' + str(e))
     if len(include_dirs) > 0:
         build_args["include_dirs"] = include_dirs
     if len(library_dirs) > 0:
@@ -124,14 +128,15 @@ try:
     else:
         ext_modules = [Extension("cobra.solvers.cglpk",
                                  ["cobra/solvers/cglpk.c"], **build_args)]
-except:
+except Exception as e:
+    print('Could not build CGLPK: {}'.format(e))
     ext_modules = None
 
 extras = {
     'matlab': ["pymatbridge"],
     'sbml': ["python-libsbml", "lxml"],
     'array': ["numpy>=1.6", "scipy>=11.0"],
-    'display': ["matplotlib", "brewer2mpl", "pandas"]
+    'display': ["matplotlib", "palettable", "pandas"]
 }
 
 all_extras = {'Cython>=0.21'}
@@ -175,7 +180,7 @@ setup(
 
     author="Daniel Robert Hyduke <danielhyduke at gmail.com>, "
     "Ali Ebrahim <aebrahim at ucsd.edu>",
-    author_email="danielhyduke at gmail.com",
+    author_email="aebrahim at ucsd.edu",
     description="COBRApy is a package for constraints-based modeling of "
     "biological networks",
     license="LGPL/GPL v2+",
@@ -195,6 +200,7 @@ setup(
         'Operating System :: OS Independent',
         'Programming Language :: Python :: 2.7',
         'Programming Language :: Python :: 3.4',
+        'Programming Language :: Python :: 3.5',
         'Programming Language :: Cython',
         'Programming Language :: Python :: Implementation :: CPython',
         'Topic :: Scientific/Engineering',

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-cobra.git



More information about the debian-med-commit mailing list