diff --git a/src/cobra/core/metabolite.py b/src/cobra/core/metabolite.py
index ae68da9b1..ef918b210 100644
--- a/src/cobra/core/metabolite.py
+++ b/src/cobra/core/metabolite.py
@@ -73,8 +73,8 @@ def constraint(self):
 
     @property
     def elements(self):
-        """ Dictionary of elements as keys and their count in the metabolite
-        as integer. When set, the `formula` property is update accordingly """
+        """Dictionary of elements as keys and their count in the metabolite
+        as integer. When set, the `formula` property is update accordingly"""
         tmp_formula = self.formula
         if tmp_formula is None:
             return {}
@@ -246,7 +246,10 @@ def summary(self, solution=None, fva=None):
         from cobra.summary import MetaboliteSummary
 
         return MetaboliteSummary(
-            metabolite=self, model=self._model, solution=solution, fva=fva,
+            metabolite=self,
+            model=self._model,
+            solution=solution,
+            fva=fva,
         )
 
     def _repr_html_(self):
diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index 95b89c0aa..d92412246 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -73,8 +73,7 @@ class Model(Object):
     """
 
     def __setstate__(self, state):
-        """Make sure all cobra.Objects in the model point to the model.
-        """
+        """Make sure all cobra.Objects in the model point to the model."""
         self.__dict__.update(state)
         for y in ["reactions", "genes", "metabolites"]:
             for x in getattr(self, y):
diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index b08609dd0..8fe9e131a 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -130,7 +130,7 @@ def flux_expression(self):
             The expression representing the the forward flux (if associated
             with model), otherwise None. Representing the net flux if
             model.reversible_encoding == 'unsplit' or None if reaction is
-            not associated with a model """
+            not associated with a model"""
         if self.model is not None:
             return 1.0 * self.forward_variable - 1.0 * self.reverse_variable
         else:
@@ -169,7 +169,7 @@ def reverse_variable(self):
 
     @property
     def objective_coefficient(self):
-        """ Get the coefficient for this reaction in a linear
+        """Get the coefficient for this reaction in a linear
         objective (float)
 
         Assuming that the objective of the associated model is summation of
@@ -294,7 +294,7 @@ def upper_bound(self, value):
 
     @property
     def bounds(self):
-        """ Get or set the bounds directly from a tuple
+        """Get or set the bounds directly from a tuple
 
         Convenience method for setting upper and lower bounds in one line
         using a tuple of lower and upper bound. Invalid bounds will raise an
@@ -1174,7 +1174,10 @@ def summary(self, solution=None, fva=None):
         from cobra.summary import ReactionSummary
 
         return ReactionSummary(
-            reaction=self, model=self._model, solution=solution, fva=fva,
+            reaction=self,
+            model=self._model,
+            solution=solution,
+            fva=fva,
         )
 
     def __str__(self):
diff --git a/src/cobra/flux_analysis/deletion.py b/src/cobra/flux_analysis/deletion.py
index f691ed9c0..195bdaae4 100644
--- a/src/cobra/flux_analysis/deletion.py
+++ b/src/cobra/flux_analysis/deletion.py
@@ -132,7 +132,14 @@ def _multi_deletion(
 
         def extract_knockout_results(result_iter):
             result = pd.DataFrame(
-                [(set(ids), growth, status,) for (ids, growth, status) in result_iter],
+                [
+                    (
+                        set(ids),
+                        growth,
+                        status,
+                    )
+                    for (ids, growth, status) in result_iter
+                ],
                 columns=["ids", "growth", "status"],
             )
             return result
diff --git a/src/cobra/flux_analysis/gapfilling.py b/src/cobra/flux_analysis/gapfilling.py
index bc1fa1de2..3e5c0a414 100644
--- a/src/cobra/flux_analysis/gapfilling.py
+++ b/src/cobra/flux_analysis/gapfilling.py
@@ -80,7 +80,7 @@ class GapFiller(object):
        e1004808. doi:10.1371/journal.pcbi.1004808.
 
        [5] Diener, Christian https://github.com/cdiener/corda
-     """
+    """
 
     def __init__(
         self,
@@ -178,8 +178,7 @@ def update_costs(self):
         self.model.objective.set_linear_coefficients(self.costs)
 
     def add_switches_and_objective(self):
-        """ Update gapfilling model with switches and the indicator objective.
-        """
+        """Update gapfilling model with switches and the indicator objective."""
         constraints = list()
         big_m = max(max(abs(b) for b in r.bounds) for r in self.model.reactions)
         prob = self.model.problem
@@ -291,54 +290,54 @@ def gapfill(
 ):
     """Perform gapfilling on a model.
 
-    See documentation for the class GapFiller.
+        See documentation for the class GapFiller.
 
-    Parameters
-    ----------
-    model : cobra.Model
-        The model to perform gap filling on.
-    universal : cobra.Model, None
-        A universal model with reactions that can be used to complete the
-        model. Only gapfill considering demand and exchange reactions if
-        left missing.
-    lower_bound : float
-        The minimally accepted flux for the objective in the filled model.
-    penalties : dict, None
-        A dictionary with keys being 'universal' (all reactions included in
-        the universal model), 'exchange' and 'demand' (all additionally
-        added exchange and demand reactions) for the three reaction types.
-        Can also have reaction identifiers for reaction specific costs.
-        Defaults are 1, 100 and 1 respectively.
-    iterations : int
-        The number of rounds of gapfilling to perform. For every iteration,
-        the penalty for every used reaction increases linearly. This way,
-        the algorithm is encouraged to search for alternative solutions
-        which may include previously used reactions. I.e., with enough
-        iterations pathways including 10 steps will eventually be reported
-        even if the shortest pathway is a single reaction.
-    exchange_reactions : bool
-        Consider adding exchange (uptake) reactions for all metabolites
-        in the model.
-    demand_reactions : bool
-        Consider adding demand reactions for all metabolites.
+        Parameters
+        ----------
+        model : cobra.Model
+            The model to perform gap filling on.
+        universal : cobra.Model, None
+            A universal model with reactions that can be used to complete the
+            model. Only gapfill considering demand and exchange reactions if
+            left missing.
+        lower_bound : float
+            The minimally accepted flux for the objective in the filled model.
+        penalties : dict, None
+            A dictionary with keys being 'universal' (all reactions included in
+            the universal model), 'exchange' and 'demand' (all additionally
+            added exchange and demand reactions) for the three reaction types.
+            Can also have reaction identifiers for reaction specific costs.
+            Defaults are 1, 100 and 1 respectively.
+        iterations : int
+            The number of rounds of gapfilling to perform. For every iteration,
+            the penalty for every used reaction increases linearly. This way,
+            the algorithm is encouraged to search for alternative solutions
+            which may include previously used reactions. I.e., with enough
+            iterations pathways including 10 steps will eventually be reported
+            even if the shortest pathway is a single reaction.
+        exchange_reactions : bool
+            Consider adding exchange (uptake) reactions for all metabolites
+            in the model.
+        demand_reactions : bool
+            Consider adding demand reactions for all metabolites.
 
-    Returns
-    -------
-    iterable
-        list of lists with on set of reactions that completes the model per
-        requested iteration.
-
-    Examples
-    --------
-
-import cobra as ct
-    >>> from cobra import Model
-    >>> from cobra.flux_analysis import gapfill
-    >>> model = ct.create_test_model("salmonella")
-    >>> universal = Model('universal')
-    >>> universal.add_reactions(model.reactions.GF6PTA.copy())
-    >>> model.remove_reactions([model.reactions.GF6PTA])
-    >>> gapfill(model, universal)
+        Returns
+        -------
+        iterable
+            list of lists with on set of reactions that completes the model per
+            requested iteration.
+
+        Examples
+        --------
+
+    import cobra as ct
+        >>> from cobra import Model
+        >>> from cobra.flux_analysis import gapfill
+        >>> model = ct.create_test_model("salmonella")
+        >>> universal = Model('universal')
+        >>> universal.add_reactions(model.reactions.GF6PTA.copy())
+        >>> model.remove_reactions([model.reactions.GF6PTA])
+        >>> gapfill(model, universal)
     """
     gapfiller = GapFiller(
         model,
diff --git a/src/cobra/flux_analysis/loopless.py b/src/cobra/flux_analysis/loopless.py
index 36c8dfb2b..203135d55 100644
--- a/src/cobra/flux_analysis/loopless.py
+++ b/src/cobra/flux_analysis/loopless.py
@@ -169,7 +169,9 @@ def loopless_solution(model, fluxes=None):
         prob = model.problem
         # Needs one fixed bound for cplex...
         loopless_obj_constraint = prob.Constraint(
-            model.objective.expression, lb=-1e32, name="loopless_obj_constraint",
+            model.objective.expression,
+            lb=-1e32,
+            name="loopless_obj_constraint",
         )
         model.add_cons_vars([loopless_obj_constraint])
         _add_cycle_free(model, fluxes)
diff --git a/src/cobra/flux_analysis/variability.py b/src/cobra/flux_analysis/variability.py
index 3b5c52275..471fc7004 100644
--- a/src/cobra/flux_analysis/variability.py
+++ b/src/cobra/flux_analysis/variability.py
@@ -186,7 +186,8 @@ def flux_variability_analysis(
         if pfba_factor is not None:
             if pfba_factor < 1.0:
                 warn(
-                    "The 'pfba_factor' should be larger or equal to 1.", UserWarning,
+                    "The 'pfba_factor' should be larger or equal to 1.",
+                    UserWarning,
                 )
             with model:
                 add_pfba(model, fraction_of_optimum=0)
@@ -227,7 +228,11 @@ def flux_variability_analysis(
 
 
 def find_blocked_reactions(
-    model, reaction_list=None, zero_cutoff=None, open_exchanges=False, processes=None,
+    model,
+    reaction_list=None,
+    zero_cutoff=None,
+    open_exchanges=False,
+    processes=None,
 ):
     """
     Find reactions that cannot carry any flux.
diff --git a/src/cobra/io/sbml.py b/src/cobra/io/sbml.py
index 3351f5a1c..1fea29157 100644
--- a/src/cobra/io/sbml.py
+++ b/src/cobra/io/sbml.py
@@ -111,7 +111,7 @@ class CobraSBMLError(Exception):
 
 def _escape_non_alphanum(nonASCII):
     """converts a non alphanumeric character to a string representation of
-    its ascii number """
+    its ascii number"""
     return "__" + str(ord(nonASCII.group())) + "__"
 
 
@@ -567,7 +567,7 @@ def _sbml_to_model(
 
     # GPR rules
     def process_association(ass):
-        """ Recursively convert gpr association to a gpr string.
+        """Recursively convert gpr association to a gpr string.
         Defined as inline functions to not pass the replacement dict around.
         """
         if ass.isFbcOr():
@@ -1406,7 +1406,7 @@ def _check(value, message):
 # Notes
 # -----------------------------------------------------------------------------
 def _parse_notes_dict(sbase):
-    """ Creates dictionary of COBRA notes.
+    """Creates dictionary of COBRA notes.
 
     Parameters
     ----------
diff --git a/src/cobra/sampling/achr.py b/src/cobra/sampling/achr.py
index 6baa7c00b..616a81863 100644
--- a/src/cobra/sampling/achr.py
+++ b/src/cobra/sampling/achr.py
@@ -154,7 +154,8 @@ def sample(self, n, fluxes=True):
             names = [r.id for r in self.model.reactions]
 
             return pandas.DataFrame(
-                samples[:, self.fwd_idx] - samples[:, self.rev_idx], columns=names,
+                samples[:, self.fwd_idx] - samples[:, self.rev_idx],
+                columns=names,
             )
         else:
             names = [v.name for v in self.model.variables]
diff --git a/src/cobra/sampling/hr_sampler.py b/src/cobra/sampling/hr_sampler.py
index cd2c1e181..70d930650 100644
--- a/src/cobra/sampling/hr_sampler.py
+++ b/src/cobra/sampling/hr_sampler.py
@@ -285,7 +285,9 @@ def generate_fva_warmup(self):
 
                 primals = self.model.solver.primal_values
                 sol = [primals[v.name] for v in self.model.variables]
-                self.warmup[self.n_warmup,] = sol
+                self.warmup[
+                    self.n_warmup,
+                ] = sol
                 self.n_warmup += 1
 
                 # Reset objective
@@ -391,13 +393,33 @@ def _bounds_dist(self, p):
         """Get the lower and upper bound distances. Negative is bad."""
 
         prob = self.problem
-        lb_dist = (p - prob.variable_bounds[0,]).min()
-        ub_dist = (prob.variable_bounds[1,] - p).min()
+        lb_dist = (
+            p
+            - prob.variable_bounds[
+                0,
+            ]
+        ).min()
+        ub_dist = (
+            prob.variable_bounds[
+                1,
+            ]
+            - p
+        ).min()
 
         if prob.bounds.shape[0] > 0:
             const = prob.inequalities.dot(p)
-            const_lb_dist = (const - prob.bounds[0,]).min()
-            const_ub_dist = (prob.bounds[1,] - const).min()
+            const_lb_dist = (
+                const
+                - prob.bounds[
+                    0,
+                ]
+            ).min()
+            const_ub_dist = (
+                prob.bounds[
+                    1,
+                ]
+                - const
+            ).min()
             lb_dist = min(lb_dist, const_lb_dist)
             ub_dist = min(ub_dist, const_ub_dist)
 
@@ -486,13 +508,39 @@ def validate(self, samples):
             )
 
         feasibility = np.abs(S.dot(samples.T).T - b).max(axis=1)
-        lb_error = (samples - bounds[0,]).min(axis=1)
-        ub_error = (bounds[1,] - samples).min(axis=1)
+        lb_error = (
+            samples
+            - bounds[
+                0,
+            ]
+        ).min(axis=1)
+        ub_error = (
+            bounds[
+                1,
+            ]
+            - samples
+        ).min(axis=1)
 
         if samples.shape[1] == len(self.model.variables) and prob.inequalities.shape[0]:
             consts = prob.inequalities.dot(samples.T)
-            lb_error = np.minimum(lb_error, (consts - prob.bounds[0,]).min(axis=1))
-            ub_error = np.minimum(ub_error, (prob.bounds[1,] - consts).min(axis=1))
+            lb_error = np.minimum(
+                lb_error,
+                (
+                    consts
+                    - prob.bounds[
+                        0,
+                    ]
+                ).min(axis=1),
+            )
+            ub_error = np.minimum(
+                ub_error,
+                (
+                    prob.bounds[
+                        1,
+                    ]
+                    - consts
+                ).min(axis=1),
+            )
 
         valid = (
             (feasibility < self.feasibility_tol)
diff --git a/src/cobra/sampling/optgp.py b/src/cobra/sampling/optgp.py
index bd5a07e11..b34450a93 100644
--- a/src/cobra/sampling/optgp.py
+++ b/src/cobra/sampling/optgp.py
@@ -234,7 +234,8 @@ def sample(self, n, fluxes=True):
             names = [r.id for r in self.model.reactions]
 
             return pandas.DataFrame(
-                chains[:, self.fwd_idx] - chains[:, self.rev_idx], columns=names,
+                chains[:, self.fwd_idx] - chains[:, self.rev_idx],
+                columns=names,
             )
         else:
             names = [v.name for v in self.model.variables]
diff --git a/src/cobra/summary/metabolite_summary.py b/src/cobra/summary/metabolite_summary.py
index 4bc30db1c..f9fdef4ac 100644
--- a/src/cobra/summary/metabolite_summary.py
+++ b/src/cobra/summary/metabolite_summary.py
@@ -128,7 +128,11 @@ def _generate(
         # Create the basic flux table.
         flux = pd.DataFrame(
             data=[
-                (r.id, solution[r.id], r.get_coefficient(self._metabolite.id),)
+                (
+                    r.id,
+                    solution[r.id],
+                    r.get_coefficient(self._metabolite.id),
+                )
                 for r in self._reactions
             ],
             columns=["reaction", "flux", "factor"],
@@ -394,11 +398,13 @@ def to_html(
             metabolite = self._metabolite.id
 
         production = self._html_table(
-            self._display_flux(self.producing_flux, names, threshold), float_format,
+            self._display_flux(self.producing_flux, names, threshold),
+            float_format,
         )
 
         consumption = self._html_table(
-            self._display_flux(self.consuming_flux, names, threshold), float_format,
+            self._display_flux(self.consuming_flux, names, threshold),
+            float_format,
         )
 
         return (
diff --git a/src/cobra/summary/model_summary.py b/src/cobra/summary/model_summary.py
index c35b72320..832fe6f15 100644
--- a/src/cobra/summary/model_summary.py
+++ b/src/cobra/summary/model_summary.py
@@ -118,7 +118,9 @@ def _generate(
         if isinstance(fva, float):
             logger.info("Performing flux variability analysis.")
             fva = flux_variability_analysis(
-                model=model, reaction_list=model.boundary, fraction_of_optimum=fva,
+                model=model,
+                reaction_list=model.boundary,
+                fraction_of_optimum=fva,
             )
         if coefficients:
             self._objective: Dict["Reaction", float] = {
diff --git a/src/cobra/summary/reaction_summary.py b/src/cobra/summary/reaction_summary.py
index 45a01507d..a8e7f20ba 100644
--- a/src/cobra/summary/reaction_summary.py
+++ b/src/cobra/summary/reaction_summary.py
@@ -107,17 +107,24 @@ def _generate(
         if isinstance(fva, float):
             logger.info("Performing flux variability analysis.")
             fva = flux_variability_analysis(
-                model, reaction_list=[self._reaction], fraction_of_optimum=fva,
+                model,
+                reaction_list=[self._reaction],
+                fraction_of_optimum=fva,
             )
 
         # Create the basic flux table.
         self._flux = pd.DataFrame(
-            data={"flux": [solution[self._reaction.id]]}, index=[self._reaction.id],
+            data={"flux": [solution[self._reaction.id]]},
+            index=[self._reaction.id],
         )
         if fva is not None:
             self._flux = self._flux.join(fva)
 
-    def _string_flux(self, threshold: float, float_format: str,) -> str:
+    def _string_flux(
+        self,
+        threshold: float,
+        float_format: str,
+    ) -> str:
         """
         Transform a flux data frame to a string.
 
diff --git a/src/cobra/summary/summary.py b/src/cobra/summary/summary.py
index f975a9873..e7403b051 100644
--- a/src/cobra/summary/summary.py
+++ b/src/cobra/summary/summary.py
@@ -27,7 +27,10 @@ class Summary(ABC):
 
     """
 
-    def __init__(self, **kwargs,) -> None:
+    def __init__(
+        self,
+        **kwargs,
+    ) -> None:
         """
         Initialize a summary.
 
diff --git a/src/cobra/test/__init__.py b/src/cobra/test/__init__.py
index 4ae886e20..a3247c446 100644
--- a/src/cobra/test/__init__.py
+++ b/src/cobra/test/__init__.py
@@ -46,8 +46,7 @@ def create_test_model(model_name="salmonella"):
 
 
 def test_all(args=None):
-    """ alias for running all unit-tests on installed cobra
-    """
+    """alias for running all unit-tests on installed cobra"""
     if pytest:
         args = args if args else []
 
diff --git a/src/cobra/test/test_summary/test_metabolite_summary.py b/src/cobra/test/test_summary/test_metabolite_summary.py
index 5b49005ac..1f543bfd0 100644
--- a/src/cobra/test/test_summary/test_metabolite_summary.py
+++ b/src/cobra/test/test_summary/test_metabolite_summary.py
@@ -12,13 +12,18 @@ def test_metabolite_summary_interface(model, opt_solver):
     model.solver = opt_solver
     metabolite = model.metabolites.get_by_id("q8_c")
     MetaboliteSummary(
-        metabolite=metabolite, model=model,
+        metabolite=metabolite,
+        model=model,
     )
     MetaboliteSummary(
-        metabolite=metabolite, model=model, solution=pfba(model),
+        metabolite=metabolite,
+        model=model,
+        solution=pfba(model),
     )
     MetaboliteSummary(
-        metabolite=metabolite, model=model, fva=0.95,
+        metabolite=metabolite,
+        model=model,
+        fva=0.95,
     )
     MetaboliteSummary(
         metabolite=metabolite,
diff --git a/src/cobra/test/test_summary/test_model_summary.py b/src/cobra/test/test_summary/test_model_summary.py
index 7300642a3..c05e8ab00 100644
--- a/src/cobra/test/test_summary/test_model_summary.py
+++ b/src/cobra/test/test_summary/test_model_summary.py
@@ -10,12 +10,16 @@
 def test_model_summary_interface(model, opt_solver):
     """Test that a summary can be created successfully."""
     model.solver = opt_solver
-    ModelSummary(model=model,)
     ModelSummary(
-        model=model, solution=pfba(model),
+        model=model,
     )
     ModelSummary(
-        model=model, fva=0.95,
+        model=model,
+        solution=pfba(model),
+    )
+    ModelSummary(
+        model=model,
+        fva=0.95,
     )
     ModelSummary(
         model=model,
diff --git a/src/cobra/test/test_summary/test_reaction_summary.py b/src/cobra/test/test_summary/test_reaction_summary.py
index 7f11d3219..8926285a2 100644
--- a/src/cobra/test/test_summary/test_reaction_summary.py
+++ b/src/cobra/test/test_summary/test_reaction_summary.py
@@ -12,13 +12,18 @@ def test_reaction_summary_interface(model, opt_solver):
     model.solver = opt_solver
     reaction = model.reactions.get_by_id("FUM")
     ReactionSummary(
-        reaction=reaction, model=model,
+        reaction=reaction,
+        model=model,
     )
     ReactionSummary(
-        reaction=reaction, model=model, solution=pfba(model),
+        reaction=reaction,
+        model=model,
+        solution=pfba(model),
     )
     ReactionSummary(
-        reaction=reaction, model=model, fva=0.95,
+        reaction=reaction,
+        model=model,
+        fva=0.95,
     )
     ReactionSummary(
         reaction=reaction,
diff --git a/src/cobra/util/array.py b/src/cobra/util/array.py
index 79972fb19..b14cb2cc8 100644
--- a/src/cobra/util/array.py
+++ b/src/cobra/util/array.py
@@ -136,7 +136,9 @@ def nullspace(A: np.ndarray, atol: float = 1e-13, rtol: float = 0.0) -> np.ndarr
 
 
 def constraint_matrices(
-    model: "Model", array_type: str = "dense", zero_tol: float = 1e-6,
+    model: "Model",
+    array_type: str = "dense",
+    zero_tol: float = 1e-6,
 ) -> NamedTuple:
     """Create a matrix representation of the problem.
 
diff --git a/src/cobra/util/util.py b/src/cobra/util/util.py
index 9754d8f5e..154552353 100644
--- a/src/cobra/util/util.py
+++ b/src/cobra/util/util.py
@@ -12,7 +12,7 @@ def format_long_string(string, max_length=50):
 
 class AutoVivification(dict):
     """Implementation of perl's autovivification feature. Checkout
-    http://stackoverflow.com/a/652284/280182 """
+    http://stackoverflow.com/a/652284/280182"""
 
     def __getitem__(self, item):
         try: