From f291aec74482cbd75fdc9fe23a394fd79fd91c4c Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Thu, 28 Jul 2022 21:42:26 -0400
Subject: [PATCH 1/6] solving the problem in 673 and adding test

---
 src/cobra/core/model.py       | 38 +++++++++++++++++++++--
 src/cobra/core/reaction.py    |  2 +-
 src/cobra/core/species.py     |  9 ++++--
 tests/test_core/test_model.py | 58 +++++++++++++++++++++++++++++++++++
 4 files changed, 102 insertions(+), 5 deletions(-)

diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index 597c3802b..fd934f0f1 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -500,6 +500,9 @@ def add_metabolites(self, metabolite_list: Union[List, Metabolite]) -> None:
 
         # First check whether the metabolites exist in the model
         metabolite_list = [x for x in metabolite_list if x.id not in self.metabolites]
+        metabolite_list = [
+            x.copy() if x.model and x.model != self else x for x in metabolite_list
+        ]
 
         bad_ids = [
             m for m in metabolite_list if not isinstance(m.id, str) or len(m.id) < 1
@@ -718,8 +721,39 @@ def existing_filter(rxn: Reaction) -> bool:
                 return False
             return True
 
-        # First check whether the reactions exist in the model.
-        pruned = DictList(filter(existing_filter, reaction_list))
+        # First check whether the reactions exist in the model. Copy reactions that
+        # belong to another model, in order not to remove them (see
+        # https://github.com/opencobra/cobrapy/issues/673#issuecomment-371361476).
+        exist_but_different_stoich = set()
+        other_model = ""
+        for reaction in filter(lambda x: not existing_filter(x), reaction_list):
+            met_id_dict = {
+                met.id: stoich for met, stoich in reaction.metabolites.items()
+            }
+            existing_met_id_dict = {
+                met.id: stoich
+                for met, stoich in self.reactions.get_by_id(
+                    reaction.id
+                ).metabolites.items()
+            }
+            if met_id_dict != existing_met_id_dict:
+                exist_but_different_stoich.add(reaction.id)
+                other_model = reaction.model
+        if len(exist_but_different_stoich):
+            logger.warning(
+                f"The reactions {', '.join(exist_but_different_stoich)} exist"
+                f" in both {self} and {other_model}, but have different "
+                f"stoichiometries in the two mdoels. Please check to see "
+                f"which stoichiometry you'd like to add."
+            )
+        pruned = DictList(
+            [
+                reaction
+                if not reaction._model or reaction._model == self
+                else reaction.copy()
+                for reaction in filter(existing_filter, reaction_list)
+            ]
+        )
 
         context = get_context(self)
 
diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index 0623ae170..a703d5d45 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -1232,7 +1232,7 @@ def add_metabolites(
             # Make sure metabolites being added belong to the same model, or
             # else copy them.
             if isinstance(metabolite, Metabolite):
-                if (metabolite.model is not None) and (
+                if (metabolite.model is not None) and self.model and (
                     metabolite.model is not self._model
                 ):
                     metabolite = metabolite.copy()
diff --git a/src/cobra/core/species.py b/src/cobra/core/species.py
index 39749fcd0..849c7bfd2 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -2,13 +2,14 @@
 
 
 from copy import deepcopy
+from functools import partial
 from typing import TYPE_CHECKING, FrozenSet, Optional
 
 from ..core.object import Object
-
+from ..util import HistoryManager
 
 if TYPE_CHECKING:
-    from .. import Model
+    from .. import Model, Reaction
 
 
 class Species(Object):
@@ -55,6 +56,10 @@ def reactions(self) -> FrozenSet:
         FrozenSet
             A frozenset that includes the reactions of the species.
         """
+        if self.model and self.__class__.__name__ == "Gene":
+            return self.model.reactions.query(lambda x: self in x, 'genes')
+        elif self.model and self.__class__.__name__ == "Metabolite":
+            return self.model.reactions.query(lambda x: self in x, 'metabolites')
         return frozenset(self._reaction)
 
     def __getstate__(self) -> dict:
diff --git a/tests/test_core/test_model.py b/tests/test_core/test_model.py
index e7e606ae1..264d28972 100644
--- a/tests/test_core/test_model.py
+++ b/tests/test_core/test_model.py
@@ -4,6 +4,7 @@
 import warnings
 from copy import copy, deepcopy
 from math import isnan
+from pathlib import Path
 from typing import List, Tuple
 
 import numpy as np
@@ -14,6 +15,7 @@
 from cobra import Solution
 from cobra.core import Group, Metabolite, Model, Reaction
 from cobra.exceptions import OptimizationError
+from cobra.io import read_sbml_model
 from cobra.manipulation.delete import remove_genes
 from cobra.util import solver as su
 from cobra.util.solver import SolverNotFound, set_objective, solvers
@@ -982,6 +984,62 @@ def test_add_reactions(model: Model) -> None:
     assert coefficients_dict[model.reactions.r2.reverse_variable] == -3.0
 
 
+def test_add_reactions_with_context(model: Model) -> None:
+    """Test add_reactions() function to add reactions to model in a context.
+
+    Parameters
+    ----------
+    model: cobra.Model
+    """
+    r1 = Reaction("r1")
+    r1.add_metabolites({Metabolite("A"): -1, Metabolite("B"): 1})
+    r1.lower_bound, r1.upper_bound = -999999.0, 999999.0
+    r2 = Reaction("r2")
+    r2.add_metabolites({Metabolite("A"): -1, Metabolite("C"): 1, Metabolite("D"): 1})
+    r2.lower_bound, r2.upper_bound = 0.0, 999999.0
+    with model:
+        model.add_reactions([r1, r2])
+        r2.objective_coefficient = 3.0
+        assert r2.objective_coefficient == 3.0
+        assert model.reactions[-2] == r1
+        assert model.reactions[-1] == r2
+        assert isinstance(model.reactions[-2].reverse_variable, model.problem.Variable)
+        coefficients_dict = model.objective.expression.as_coefficients_dict()
+        biomass_r = model.reactions.get_by_id("Biomass_Ecoli_core")
+        assert coefficients_dict[biomass_r.forward_variable] == 1.0
+        assert coefficients_dict[biomass_r.reverse_variable] == -1.0
+        assert coefficients_dict[model.reactions.r2.forward_variable] == 3.0
+        assert coefficients_dict[model.reactions.r2.reverse_variable] == -3.0
+    assert model.reactions[-2] != r1
+    assert model.reactions[-1] != r2
+    assert isinstance(model.reactions[-2].reverse_variable, model.problem.Variable)
+    coefficients_dict = model.objective.expression.as_coefficients_dict()
+    biomass_r = model.reactions.get_by_id("Biomass_Ecoli_core")
+    assert coefficients_dict[biomass_r.forward_variable] == 1.0
+    assert coefficients_dict[biomass_r.reverse_variable] == -1.0
+    assert r2.forward_variable not in coefficients_dict
+    assert r2.reverse_variable not in coefficients_dict
+
+
+def test_add_reactions_from_another_model(model: Model, data_directory: Path) -> None:
+    """Test adding reaction from an existing model.
+
+    This will test that adding a reaction with an existing model does not add
+    reactions of the metabolites, and does not remove the reaction from the existing
+    model.
+
+    Parameters
+    ----------
+    model: corbra.Model
+    data_directory: Path
+        Directory where mini_cobra.xml is found.
+    """
+    mini = read_sbml_model(data_directory / 'mini_cobra.xml')
+    mini.add_reactions([model.reactions.get_by_id('ACALD')])
+    assert len(mini.metabolites.get_by_id('acald_c').reactions) == 1
+
+
+
 def test_add_reactions_single_existing(model: Model) -> None:
     """Test adding a reaction already present to a model.
 

From 3ed87af625a6ee6e87e3b3e58764c568d6363c9f Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Thu, 28 Jul 2022 21:56:25 -0400
Subject: [PATCH 2/6] minimize usage of private attributes added species.py
 functions to set ._reaction field made metabolites a settable property, and
 used that to simplify bits of model (including copy and add_reaction)

---
 src/cobra/core/model.py    | 74 ++++++++++++++++++--------------------
 src/cobra/core/reaction.py | 38 ++++++++++++++------
 src/cobra/core/species.py  | 34 ++++++++++++++++++
 3 files changed, 96 insertions(+), 50 deletions(-)

diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index fd934f0f1..84524c321 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -428,10 +428,11 @@ def copy(self) -> "Model":
             new_reaction._model = new
             new.reactions.append(new_reaction)
             # update awareness
-            for metabolite, stoic in reaction._metabolites.items():
-                new_met = new.metabolites.get_by_id(metabolite.id)
-                new_reaction._metabolites[new_met] = stoic
-                new_met._reaction.add(new_reaction)
+            new_metabolites = {
+                new.metabolites.get_by_id(metabolite.id): stoic
+                for metabolite, stoic in reaction.metabolites.items()
+            }
+            new_reaction.metabolites = new_metabolites
             new_reaction.update_genes_from_gpr()
 
         new.groups = DictList()
@@ -561,12 +562,12 @@ def remove_metabolites(
                 group.remove_members(x)
 
             if not destructive:
-                for the_reaction in list(x._reaction):  # noqa W0212
+                for the_reaction in x.reactions:
                     the_coefficient = the_reaction._metabolites[x]  # noqa W0212
                     the_reaction.subtract_metabolites({x: the_coefficient})
 
             else:
-                for x2 in list(x._reaction):  # noqa W0212
+                for x2 in x.reactions:
                     x2.remove_from_model()
 
         self.metabolites -= metabolite_list
@@ -762,23 +763,22 @@ def existing_filter(rxn: Reaction) -> bool:
             reaction._model = self
             if context:
                 context(partial(setattr, reaction, "_model", None))
-            # Build a `list()` because the dict will be modified in the loop.
-            for metabolite in list(reaction.metabolites):
-                # TODO: Maybe this can happen with
-                #  Reaction.add_metabolites(combine=False)
-                # TODO: Should we add a copy of the metabolite instead?
-                if metabolite not in self.metabolites:
-                    self.add_metabolites(metabolite)
-                # A copy of the metabolite exists in the model, the reaction
-                # needs to point to the metabolite in the model.
-                else:
-                    # FIXME: Modifying 'private' attributes is horrible.
-                    stoichiometry = reaction._metabolites.pop(metabolite)
-                    model_metabolite = self.metabolites.get_by_id(metabolite.id)
-                    reaction._metabolites[model_metabolite] = stoichiometry
-                    model_metabolite._reaction.add(reaction)
-                    if context:
-                        context(partial(model_metabolite._reaction.remove, reaction))
+            met_id_dict = {
+                met.id: stoich for met, stoich in reaction.metabolites.items()
+            }
+            met_dict = {met.id: met for met in reaction.metabolites.keys()}
+            mets_to_add = [
+                met_dict[met_id]
+                for met_id in set(met_id_dict.keys()).difference(
+                    self.metabolites.list_attr("id")
+                )
+            ]
+            self.add_metabolites(mets_to_add)
+            new_stoich = {
+                self.metabolites.get_by_id(met_id): stoich
+                for met_id, stoich in met_id_dict.items()
+            }
+            reaction.metabolites = new_stoich
             reaction.update_genes_from_gpr()
 
         self.reactions += pruned
@@ -846,21 +846,17 @@ def remove_reactions(
                 self.reactions.remove(reaction)
                 reaction._model = None
 
-                for met in reaction._metabolites:
-                    if reaction in met._reaction:
-                        met._reaction.remove(reaction)
-                        if context:
-                            context(partial(met._reaction.add, reaction))
-                        if remove_orphans and len(met._reaction) == 0:
+                for met in reaction.metabolites:
+                    if reaction in met.reactions:
+                        met.reaction_remove(reaction, context)
+                        if remove_orphans and len(met.reactions) == 0:
                             self.remove_metabolites(met)
 
-                for gene in reaction._genes:
-                    if reaction in gene._reaction:
-                        gene._reaction.remove(reaction)
-                        if context:
-                            context(partial(gene._reaction.add, reaction))
+                for gene in reaction.genes:
+                    if reaction in gene.reactions:
+                        gene.reaction_remove(reaction, context)
 
-                        if remove_orphans and len(gene._reaction) == 0:
+                        if remove_orphans and len(gene.reactions) == 0:
                             self.genes.remove(gene)
                             if context:
                                 context(partial(self.genes.add, gene))
@@ -1291,13 +1287,13 @@ def repair(
             self.groups._generate_index()
         if rebuild_relationships:
             for met in self.metabolites:
-                met._reaction.clear()
+                met.reaction_clear()
             for gene in self.genes:
-                gene._reaction.clear()
+                gene.reaction_clear()
             for rxn in self.reactions:
                 rxn.update_genes_from_gpr()
-                for met in rxn._metabolites:
-                    met._reaction.add(rxn)
+                for met in rxn.metabolites:
+                    met.reaction_add(rxn)
 
         # point _model to self
         for dict_list in (self.reactions, self.genes, self.metabolites, self.groups):
diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index a703d5d45..f234d545f 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -579,6 +579,23 @@ def metabolites(self) -> Dict[Metabolite, float]:
         """
         return self._metabolites.copy()
 
+    @metabolites.setter
+    @resettable
+    def metabolites(self, value:  Dict[Metabolite, float]) -> None:
+        """Set metabolites to a dictionary of metabolites and coefficients.
+
+        Parameters
+        -------
+        metabolites: Dict[Metabolite, float]
+            A dictionary of cobra.Metabolite for keys and floats for coeffecieints.
+            Positive coefficient means the reaction produces this metabolite, while
+            negative coefficient means the reaction consumes this metabolite.
+        """
+        self._metabolites = value
+        context = get_context(self)
+        for met in value.keys():
+            met.reaction_add(self, context)
+
     @property
     def genes(self) -> FrozenSet:
         """Return the genes of the reaction.
@@ -837,11 +854,13 @@ def _update_awareness(self) -> None:
         Make sure all metabolites and genes that are associated with
         this reaction are aware of it.
 
+        .. deprecated:: 0.26
+            This function is deprecated and shouldn't be used.
         """
         for x in self._metabolites:
-            x._reaction.add(self)
+            x.reaction_add(self)
         for x in self._genes:
-            x._reaction.add(self)
+            x.reaction_add(self)
 
     def remove_from_model(self, remove_orphans: bool = False) -> None:
         """Remove the reaction from a model.
@@ -930,12 +949,9 @@ def __setstate__(self, state: Dict) -> None:
             state["_gpr"] = GPR.from_string(state["_gpr"])
 
         self.__dict__.update(state)
-        for x in state["_metabolites"]:
-            x._model = self._model
-            x._reaction.add(self)
-        for x in state["_genes"]:
+        for x in set(state["_metabolites"]).union(state["_genes"]):
             x._model = self._model
-            x._reaction.add(self)
+            x.reaction_add(self)
 
     def copy(self) -> "Reaction":
         """Copy a reaction.
@@ -1268,7 +1284,7 @@ def add_metabolites(
                 self._metabolites[metabolite] = coefficient
                 # make the metabolite aware that it is involved in this
                 # reaction
-                metabolite._reaction.add(self)
+                metabolite.reaction_add(self)
 
         # from cameo ...
         model = self.model
@@ -1287,7 +1303,7 @@ def add_metabolites(
             if the_coefficient == 0:
                 # make the metabolite aware that it no longer participates
                 # in this reaction
-                metabolite._reaction.remove(self)
+                metabolite.reaction_remove(self)
                 self._metabolites.pop(metabolite)
 
         context = get_context(self)
@@ -1487,7 +1503,7 @@ def _associate_gene(self, cobra_gene: Gene) -> None:
 
         """
         self._genes.add(cobra_gene)
-        cobra_gene._reaction.add(self)
+        cobra_gene.reaction_add(self)
         cobra_gene._model = self._model
 
     def _dissociate_gene(self, cobra_gene: Gene) -> None:
@@ -1499,7 +1515,7 @@ def _dissociate_gene(self, cobra_gene: Gene) -> None:
 
         """
         self._genes.discard(cobra_gene)
-        cobra_gene._reaction.discard(self)
+        cobra_gene.reaction_remove(self)
 
     def knock_out(self) -> None:
         """Knockout reaction by setting its bounds to zero."""
diff --git a/src/cobra/core/species.py b/src/cobra/core/species.py
index 849c7bfd2..ec0b0f4fc 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -62,6 +62,40 @@ def reactions(self) -> FrozenSet:
             return self.model.reactions.query(lambda x: self in x, 'metabolites')
         return frozenset(self._reaction)
 
+    def reaction_add(self, reaction: Reaction, context: Optional[HistoryManager] = None) -> None:
+        """Add reaction to .reaction field, with context.
+
+        If this is called with a context, will be reversed when exiting the context.
+
+        Parmeters
+        ---------
+        reaction: cobra.Reaction
+        context: HistoryManager, optional
+            context this action is in, optional (defualt None).
+        """
+        self._reaction.add(reaction)
+        if context:
+            context(partial(self._reaction.remove, reaction))
+
+    def reaction_remove(self, reaction: Reaction, context: Optional[HistoryManager] = None) -> None:
+        """Remove reaction from .reaction field, with context.
+
+        If this is called with a context, will be reversed when exiting the context.
+
+        Parmeters
+        ---------
+        reaction: cobra.Reaction
+        context: HistoryManager, optional
+            context this action is in, optional (defualt None).
+        """
+        self._reaction.remove(reaction)
+        if context:
+            context(partial(self._reaction.add, reaction))
+
+    def reaction_clear(self) -> None:
+        """Clear the reaction field."""
+        self._reaction.clear()
+
     def __getstate__(self) -> dict:
         """Return the state of the species.
 

From 86f1c28f33ed03a25793aa66d3925671a463a71e Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Thu, 28 Jul 2022 21:58:52 -0400
Subject: [PATCH 3/6] ran black

---
 src/cobra/core/reaction.py    |  8 +++++---
 src/cobra/core/species.py     | 13 +++++++++----
 tests/test_core/test_model.py |  7 +++----
 3 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index f234d545f..2fbc2c8b1 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -581,7 +581,7 @@ def metabolites(self) -> Dict[Metabolite, float]:
 
     @metabolites.setter
     @resettable
-    def metabolites(self, value:  Dict[Metabolite, float]) -> None:
+    def metabolites(self, value: Dict[Metabolite, float]) -> None:
         """Set metabolites to a dictionary of metabolites and coefficients.
 
         Parameters
@@ -1248,8 +1248,10 @@ def add_metabolites(
             # Make sure metabolites being added belong to the same model, or
             # else copy them.
             if isinstance(metabolite, Metabolite):
-                if (metabolite.model is not None) and self.model and (
-                    metabolite.model is not self._model
+                if (
+                    (metabolite.model is not None)
+                    and self.model
+                    and (metabolite.model is not self._model)
                 ):
                     metabolite = metabolite.copy()
 
diff --git a/src/cobra/core/species.py b/src/cobra/core/species.py
index ec0b0f4fc..3ea9cef63 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -8,6 +8,7 @@
 from ..core.object import Object
 from ..util import HistoryManager
 
+
 if TYPE_CHECKING:
     from .. import Model, Reaction
 
@@ -57,12 +58,14 @@ def reactions(self) -> FrozenSet:
             A frozenset that includes the reactions of the species.
         """
         if self.model and self.__class__.__name__ == "Gene":
-            return self.model.reactions.query(lambda x: self in x, 'genes')
+            return self.model.reactions.query(lambda x: self in x, "genes")
         elif self.model and self.__class__.__name__ == "Metabolite":
-            return self.model.reactions.query(lambda x: self in x, 'metabolites')
+            return self.model.reactions.query(lambda x: self in x, "metabolites")
         return frozenset(self._reaction)
 
-    def reaction_add(self, reaction: Reaction, context: Optional[HistoryManager] = None) -> None:
+    def reaction_add(
+        self, reaction: Reaction, context: Optional[HistoryManager] = None
+    ) -> None:
         """Add reaction to .reaction field, with context.
 
         If this is called with a context, will be reversed when exiting the context.
@@ -77,7 +80,9 @@ def reaction_add(self, reaction: Reaction, context: Optional[HistoryManager] = N
         if context:
             context(partial(self._reaction.remove, reaction))
 
-    def reaction_remove(self, reaction: Reaction, context: Optional[HistoryManager] = None) -> None:
+    def reaction_remove(
+        self, reaction: Reaction, context: Optional[HistoryManager] = None
+    ) -> None:
         """Remove reaction from .reaction field, with context.
 
         If this is called with a context, will be reversed when exiting the context.
diff --git a/tests/test_core/test_model.py b/tests/test_core/test_model.py
index 264d28972..3cf981bf4 100644
--- a/tests/test_core/test_model.py
+++ b/tests/test_core/test_model.py
@@ -1034,10 +1034,9 @@ def test_add_reactions_from_another_model(model: Model, data_directory: Path) ->
     data_directory: Path
         Directory where mini_cobra.xml is found.
     """
-    mini = read_sbml_model(data_directory / 'mini_cobra.xml')
-    mini.add_reactions([model.reactions.get_by_id('ACALD')])
-    assert len(mini.metabolites.get_by_id('acald_c').reactions) == 1
-
+    mini = read_sbml_model(data_directory / "mini_cobra.xml")
+    mini.add_reactions([model.reactions.get_by_id("ACALD")])
+    assert len(mini.metabolites.get_by_id("acald_c").reactions) == 1
 
 
 def test_add_reactions_single_existing(model: Model) -> None:

From 9baeaa649e15536291adeb0eed17e2c983bb3411 Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Sun, 7 Aug 2022 19:56:02 -0400
Subject: [PATCH 4/6] fixed some bugs I introduced

---
 release-notes/next-release.md |  1 +
 src/cobra/core/model.py       |  4 ++--
 src/cobra/core/species.py     | 10 ++++++----
 3 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/release-notes/next-release.md b/release-notes/next-release.md
index 790c94af5..d2586332a 100644
--- a/release-notes/next-release.md
+++ b/release-notes/next-release.md
@@ -10,6 +10,7 @@
 of SBO:0000633 (see https://sourceforge.net/p/sbo/term-request/113/)
 * Correct reading and writing of subsystem in mat.
 * General cleanup of code in mat.py
+* Fixed issue 673 (problem when copying reaction from a different model)
 
 ## Other
 
diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index 84524c321..efa815461 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -847,13 +847,13 @@ def remove_reactions(
                 reaction._model = None
 
                 for met in reaction.metabolites:
-                    if reaction in met.reactions:
+                    if reaction in met._reaction:
                         met.reaction_remove(reaction, context)
                         if remove_orphans and len(met.reactions) == 0:
                             self.remove_metabolites(met)
 
                 for gene in reaction.genes:
-                    if reaction in gene.reactions:
+                    if reaction in gene._reaction:
                         gene.reaction_remove(reaction, context)
 
                         if remove_orphans and len(gene.reactions) == 0:
diff --git a/src/cobra/core/species.py b/src/cobra/core/species.py
index 3ea9cef63..658b5d55a 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -58,13 +58,15 @@ def reactions(self) -> FrozenSet:
             A frozenset that includes the reactions of the species.
         """
         if self.model and self.__class__.__name__ == "Gene":
-            return self.model.reactions.query(lambda x: self in x, "genes")
+            return frozenset(self.model.reactions.query(lambda x: self in x, "genes"))
         elif self.model and self.__class__.__name__ == "Metabolite":
-            return self.model.reactions.query(lambda x: self in x, "metabolites")
+            return frozenset(
+                self.model.reactions.query(lambda x: self in x, "metabolites")
+            )
         return frozenset(self._reaction)
 
     def reaction_add(
-        self, reaction: Reaction, context: Optional[HistoryManager] = None
+        self, reaction: "Reaction", context: Optional[HistoryManager] = None
     ) -> None:
         """Add reaction to .reaction field, with context.
 
@@ -81,7 +83,7 @@ def reaction_add(
             context(partial(self._reaction.remove, reaction))
 
     def reaction_remove(
-        self, reaction: Reaction, context: Optional[HistoryManager] = None
+        self, reaction: "Reaction", context: Optional[HistoryManager] = None
     ) -> None:
         """Remove reaction from .reaction field, with context.
 

From 7ea50ad4c903c37e078e23d77c2d42ee815f249b Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Sun, 7 Aug 2022 19:59:01 -0400
Subject: [PATCH 5/6] minor fixes

---
 src/cobra/core/model.py    |  2 +-
 src/cobra/core/reaction.py | 19 +++++++++----------
 2 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index efa815461..ab0d40d8a 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -899,7 +899,7 @@ def existing_filter(new_group: Group) -> bool:
             """
             if new_group.id in self.groups:
                 logger.warning(
-                    f"Ignoring group '{new_group.id}'" f" since it already exists."
+                    f"Ignoring group '{new_group.id}' since it already exists."
                 )
                 return False
             return True
diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index 2fbc2c8b1..dc2c9a2c4 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -585,8 +585,8 @@ def metabolites(self, value: Dict[Metabolite, float]) -> None:
         """Set metabolites to a dictionary of metabolites and coefficients.
 
         Parameters
-        -------
-        metabolites: Dict[Metabolite, float]
+        ----------
+        value: Dict[Metabolite, float]
             A dictionary of cobra.Metabolite for keys and floats for coeffecieints.
             Positive coefficient means the reaction produces this metabolite, while
             negative coefficient means the reaction consumes this metabolite.
@@ -1034,7 +1034,7 @@ def __iadd__(self, other: "Reaction") -> "Reaction":
         rule2 = other.gene_reaction_rule.strip()
         if rule1 != "" and rule2 != "":
             self.gene_reaction_rule = (
-                f"({self.gene_reaction_rule}) and " f"({other.gene_reaction_rule})"
+                f"({self.gene_reaction_rule}) and ({other.gene_reaction_rule})"
             )
         elif rule1 != "" and rule2 == "":
             self.gene_reaction_rule = rule1
@@ -1247,13 +1247,12 @@ def add_metabolites(
 
             # Make sure metabolites being added belong to the same model, or
             # else copy them.
-            if isinstance(metabolite, Metabolite):
-                if (
-                    (metabolite.model is not None)
-                    and self.model
-                    and (metabolite.model is not self._model)
-                ):
-                    metabolite = metabolite.copy()
+            if isinstance(metabolite, Metabolite) and (
+                (metabolite.model is not None)
+                and self.model
+                and (metabolite.model is not self._model)
+            ):
+                metabolite = metabolite.copy()
 
             met_id = str(metabolite)
             # If a metabolite already exists in the reaction then

From 63ce0577f85cedac8531b79865d73dd5977d1f12 Mon Sep 17 00:00:00 2001
From: "uri.akavia" <uri.akavia@mcgill.ca>
Date: Thu, 11 Aug 2022 13:07:55 -0400
Subject: [PATCH 6/6] fixed according to comments

---
 src/cobra/core/gene.py       | 26 +++++++++++++--
 src/cobra/core/metabolite.py | 24 ++++++++++++--
 src/cobra/core/model.py      | 59 +++++++++++++--------------------
 src/cobra/core/reaction.py   | 64 +++++++++++++++++++++++++++++++-----
 src/cobra/core/species.py    | 30 +++--------------
 5 files changed, 128 insertions(+), 75 deletions(-)

diff --git a/src/cobra/core/gene.py b/src/cobra/core/gene.py
index 0d1f86f29..5f4212dee 100644
--- a/src/cobra/core/gene.py
+++ b/src/cobra/core/gene.py
@@ -18,7 +18,7 @@
 from ast import parse as ast_parse
 from copy import deepcopy
 from keyword import kwlist
-from typing import FrozenSet, Iterable, Optional, Set, Tuple, Union
+from typing import TYPE_CHECKING, FrozenSet, Iterable, Optional, Set, Tuple, Union
 from warnings import warn
 
 import sympy.logic.boolalg as spl
@@ -30,6 +30,10 @@
 from cobra.util.util import format_long_string
 
 
+if TYPE_CHECKING:
+    from cobra.core import Reaction
+
+
 # TODO - When https://github.com/symengine/symengine.py/issues/334 is resolved,
 #  change sympy.Symbol (above in imports) to optlang.symbolics.Symbol
 
@@ -246,10 +250,28 @@ def knock_out(self) -> None:
         context.
         """
         self.functional = False
-        for reaction in self.reactions:
+        for reaction in self.reactions:  # type: ignore
             if not reaction.functional:
                 reaction.bounds = (0, 0)
 
+    @property
+    def reactions(self) -> FrozenSet["Reaction"]:
+        """Return a frozenset of reactions.
+
+        If the gene is present in a model, calculates this set by querying
+        the model reactions for this gene.
+        If not, returns the _reaction field, see cobra.Species.
+
+        Returns
+        -------
+        FrozenSet
+            A frozenset that includes the reactions of the gene.
+        """
+        if self.model:
+            return frozenset(self.model.reactions.query(lambda x: self in x, "genes"))
+        else:
+            return frozenset(self._reaction)
+
     def _repr_html_(self):
         return f"""
         <table>
diff --git a/src/cobra/core/metabolite.py b/src/cobra/core/metabolite.py
index 6b3d31dc9..69be504c2 100644
--- a/src/cobra/core/metabolite.py
+++ b/src/cobra/core/metabolite.py
@@ -1,7 +1,7 @@
 """Define the Metabolite class."""
 
 import re
-from typing import TYPE_CHECKING, Dict, Optional, Union
+from typing import TYPE_CHECKING, Dict, FrozenSet, Optional, Union
 from warnings import warn
 
 from ..exceptions import OptimizationError
@@ -15,7 +15,7 @@
     from optlang.interface import Container
     from pandas import DataFrame
 
-    from cobra.core import Solution
+    from cobra.core import Reaction, Solution
     from cobra.summary.metabolite_summary import MetaboliteSummary
 
 # Numbers are not required because of the |(?=[A-Z])? block. See the
@@ -276,6 +276,26 @@ def remove_from_model(self, destructive: bool = False) -> None:
         """
         self._model.remove_metabolites(self, destructive)
 
+    @property
+    def reactions(self) -> FrozenSet["Reaction"]:
+        """Return a frozenset of reactions.
+
+        If the metabolite is present in a model, calculates this set by querying
+        the model reactions for this metabolite.
+        If not, returns the _reaction field, see cobra.Species.
+
+        Returns
+        -------
+        FrozenSet
+            A frozenset that includes the reactions of the metabolite.
+        """
+        if self.model:
+            return frozenset(
+                self.model.reactions.query(lambda x: self in x, "metabolites")
+            )
+        else:
+            return frozenset(self._reaction)
+
     def summary(
         self,
         solution: Optional["Solution"] = None,
diff --git a/src/cobra/core/model.py b/src/cobra/core/model.py
index ab0d40d8a..22851606f 100644
--- a/src/cobra/core/model.py
+++ b/src/cobra/core/model.py
@@ -723,34 +723,24 @@ def existing_filter(rxn: Reaction) -> bool:
             return True
 
         # First check whether the reactions exist in the model. Copy reactions that
-        # belong to another model, in order not to remove them (see
-        # https://github.com/opencobra/cobrapy/issues/673#issuecomment-371361476).
-        exist_but_different_stoich = set()
+        # belong to another model, in order not to remove them from the original model.
+        exist_but_different = set()
         other_model = ""
         for reaction in filter(lambda x: not existing_filter(x), reaction_list):
-            met_id_dict = {
-                met.id: stoich for met, stoich in reaction.metabolites.items()
-            }
-            existing_met_id_dict = {
-                met.id: stoich
-                for met, stoich in self.reactions.get_by_id(
-                    reaction.id
-                ).metabolites.items()
-            }
-            if met_id_dict != existing_met_id_dict:
-                exist_but_different_stoich.add(reaction.id)
-                other_model = reaction.model
-        if len(exist_but_different_stoich):
+            if reaction != self.reactions.get_by_id(reaction.id):
+                exist_but_different.add(reaction.id)
+                other_model = other_model or reaction.model
+        if len(exist_but_different):
             logger.warning(
-                f"The reactions {', '.join(exist_but_different_stoich)} exist"
+                f"The reactions {', '.join(exist_but_different)} exist"
                 f" in both {self} and {other_model}, but have different "
-                f"stoichiometries in the two mdoels. Please check to see "
-                f"which stoichiometry you'd like to add."
+                f"properties in the two mdoels. Please check to see "
+                f"which properties are the correct ones you'd like to add."
             )
         pruned = DictList(
             [
                 reaction
-                if not reaction._model or reaction._model == self
+                if not reaction.model or reaction.model == self
                 else reaction.copy()
                 for reaction in filter(existing_filter, reaction_list)
             ]
@@ -763,20 +753,15 @@ def existing_filter(rxn: Reaction) -> bool:
             reaction._model = self
             if context:
                 context(partial(setattr, reaction, "_model", None))
-            met_id_dict = {
-                met.id: stoich for met, stoich in reaction.metabolites.items()
-            }
-            met_dict = {met.id: met for met in reaction.metabolites.keys()}
             mets_to_add = [
-                met_dict[met_id]
-                for met_id in set(met_id_dict.keys()).difference(
-                    self.metabolites.list_attr("id")
-                )
+                met
+                for met in reaction.metabolites
+                if met.id not in self.metabolites.list_attr("id")
             ]
             self.add_metabolites(mets_to_add)
             new_stoich = {
-                self.metabolites.get_by_id(met_id): stoich
-                for met_id, stoich in met_id_dict.items()
+                self.metabolites.get_by_id(met.id): stoich
+                for met, stoich in reaction.metabolites.items()
             }
             reaction.metabolites = new_stoich
             reaction.update_genes_from_gpr()
@@ -846,15 +831,15 @@ def remove_reactions(
                 self.reactions.remove(reaction)
                 reaction._model = None
 
-                for met in reaction.metabolites:
+                for met in reaction._metabolites:
                     if reaction in met._reaction:
-                        met.reaction_remove(reaction, context)
+                        met.remove_reaction(reaction)
                         if remove_orphans and len(met.reactions) == 0:
                             self.remove_metabolites(met)
 
-                for gene in reaction.genes:
+                for gene in reaction._genes:
                     if reaction in gene._reaction:
-                        gene.reaction_remove(reaction, context)
+                        gene.remove_reaction(reaction)
 
                         if remove_orphans and len(gene.reactions) == 0:
                             self.genes.remove(gene)
@@ -1287,13 +1272,13 @@ def repair(
             self.groups._generate_index()
         if rebuild_relationships:
             for met in self.metabolites:
-                met.reaction_clear()
+                met.clear_reaction()
             for gene in self.genes:
-                gene.reaction_clear()
+                gene.clear_reaction()
             for rxn in self.reactions:
                 rxn.update_genes_from_gpr()
                 for met in rxn.metabolites:
-                    met.reaction_add(rxn)
+                    met.add_reaction(rxn)
 
         # point _model to self
         for dict_list in (self.reactions, self.genes, self.metabolites, self.groups):
diff --git a/src/cobra/core/reaction.py b/src/cobra/core/reaction.py
index dc2c9a2c4..1b463bae7 100644
--- a/src/cobra/core/reaction.py
+++ b/src/cobra/core/reaction.py
@@ -1,6 +1,7 @@
 """Define the Reaction class."""
 
 import hashlib
+import math
 import re
 from collections import defaultdict
 from copy import copy, deepcopy
@@ -594,7 +595,7 @@ def metabolites(self, value: Dict[Metabolite, float]) -> None:
         self._metabolites = value
         context = get_context(self)
         for met in value.keys():
-            met.reaction_add(self, context)
+            met.add_reaction(self)
 
     @property
     def genes(self) -> FrozenSet:
@@ -858,9 +859,9 @@ def _update_awareness(self) -> None:
             This function is deprecated and shouldn't be used.
         """
         for x in self._metabolites:
-            x.reaction_add(self)
+            x.add_reaction(self)
         for x in self._genes:
-            x.reaction_add(self)
+            x.add_reaction(self)
 
     def remove_from_model(self, remove_orphans: bool = False) -> None:
         """Remove the reaction from a model.
@@ -951,7 +952,7 @@ def __setstate__(self, state: Dict) -> None:
         self.__dict__.update(state)
         for x in set(state["_metabolites"]).union(state["_genes"]):
             x._model = self._model
-            x.reaction_add(self)
+            x.add_reaction(self)
 
     def copy(self) -> "Reaction":
         """Copy a reaction.
@@ -1001,7 +1002,7 @@ def __add__(self, other: "Reaction") -> "Reaction":
 
         """
         new_reaction = self.copy()
-        if other == 0:
+        if isinstance(other, int) and other == 0:
             return new_reaction
         else:
             new_reaction += other
@@ -1285,7 +1286,7 @@ def add_metabolites(
                 self._metabolites[metabolite] = coefficient
                 # make the metabolite aware that it is involved in this
                 # reaction
-                metabolite.reaction_add(self)
+                metabolite.add_reaction(self)
 
         # from cameo ...
         model = self.model
@@ -1304,7 +1305,7 @@ def add_metabolites(
             if the_coefficient == 0:
                 # make the metabolite aware that it no longer participates
                 # in this reaction
-                metabolite.reaction_remove(self)
+                metabolite.remove_reaction(self)
                 self._metabolites.pop(metabolite)
 
         context = get_context(self)
@@ -1504,7 +1505,7 @@ def _associate_gene(self, cobra_gene: Gene) -> None:
 
         """
         self._genes.add(cobra_gene)
-        cobra_gene.reaction_add(self)
+        cobra_gene.add_reaction(self)
         cobra_gene._model = self._model
 
     def _dissociate_gene(self, cobra_gene: Gene) -> None:
@@ -1516,7 +1517,7 @@ def _dissociate_gene(self, cobra_gene: Gene) -> None:
 
         """
         self._genes.discard(cobra_gene)
-        cobra_gene.reaction_remove(self)
+        cobra_gene.remove_reaction(self)
 
     def knock_out(self) -> None:
         """Knockout reaction by setting its bounds to zero."""
@@ -1673,6 +1674,51 @@ def summary(
             fva=fva,
         )
 
+    def __eq__(self, other) -> bool:
+        self_state = self.__getstate__()
+        self_state["_metabolites"] = {
+            met.id: stoic for met, stoic in self.metabolites.items()
+        }
+        self_state["_genes"] = {gene.id for gene in self.genes}
+        self_state.pop("_model")
+        other_state = other.__getstate__()
+        other_state["_metabolites"] = {
+            met.id: stoic for met, stoic in other.metabolites.items()
+        }
+        other_state["_genes"] = {gene.id for gene in other.genes}
+        other_state.pop("_model")
+
+        self_state.pop("_lower_bound")
+        self_state.pop("_upper_bound")
+        other_state.pop("_lower_bound")
+        other_state.pop("_upper_bound")
+
+        _tolerance = 1e-7
+        if self.model:
+            _tolerance = self.model.tolerance
+        elif other.model:
+            _tolerance = other.model.tolerance
+        elif self.model and other.model:
+            _tolerance = min(self.model.tolerance, other.model.tolerance)
+        if not math.isclose(self.lower_bound, other.lower_bound, abs_tol=_tolerance):
+            return False
+        elif not math.isclose(self.upper_bound, other.upper_bound, abs_tol=_tolerance):
+            return False
+        else:
+            return self_state == other_state
+
+    def __hash__(self) -> int:
+        """Define __hash__ explicitly to allow usage of reaction in sets.
+
+        Python objects that define __eq__ MUST define __hash__ explicitly, or they
+        are unhashable. __hash__ is set to the __hash__ function of cobra.Object.
+
+        Returns
+        -------
+        int
+        """
+        return Object.__hash__(self)
+
     def __str__(self) -> str:
         """Return reaction id and reaction as str.
 
diff --git a/src/cobra/core/species.py b/src/cobra/core/species.py
index 658b5d55a..162e0361c 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -2,11 +2,9 @@
 
 
 from copy import deepcopy
-from functools import partial
 from typing import TYPE_CHECKING, FrozenSet, Optional
 
 from ..core.object import Object
-from ..util import HistoryManager
 
 
 if TYPE_CHECKING:
@@ -57,49 +55,31 @@ def reactions(self) -> FrozenSet:
         FrozenSet
             A frozenset that includes the reactions of the species.
         """
-        if self.model and self.__class__.__name__ == "Gene":
-            return frozenset(self.model.reactions.query(lambda x: self in x, "genes"))
-        elif self.model and self.__class__.__name__ == "Metabolite":
-            return frozenset(
-                self.model.reactions.query(lambda x: self in x, "metabolites")
-            )
         return frozenset(self._reaction)
 
-    def reaction_add(
-        self, reaction: "Reaction", context: Optional[HistoryManager] = None
-    ) -> None:
+    def add_reaction(self, reaction: "Reaction") -> None:
         """Add reaction to .reaction field, with context.
 
-        If this is called with a context, will be reversed when exiting the context.
+        If this is called within a context, will be reversed when exiting the context.
 
         Parmeters
         ---------
         reaction: cobra.Reaction
-        context: HistoryManager, optional
-            context this action is in, optional (defualt None).
         """
         self._reaction.add(reaction)
-        if context:
-            context(partial(self._reaction.remove, reaction))
 
-    def reaction_remove(
-        self, reaction: "Reaction", context: Optional[HistoryManager] = None
-    ) -> None:
+    def remove_reaction(self, reaction: "Reaction") -> None:
         """Remove reaction from .reaction field, with context.
 
-        If this is called with a context, will be reversed when exiting the context.
+        If this is called within a context, will be reversed when exiting the context.
 
         Parmeters
         ---------
         reaction: cobra.Reaction
-        context: HistoryManager, optional
-            context this action is in, optional (defualt None).
         """
         self._reaction.remove(reaction)
-        if context:
-            context(partial(self._reaction.add, reaction))
 
-    def reaction_clear(self) -> None:
+    def clear_reaction(self) -> None:
         """Clear the reaction field."""
         self._reaction.clear()