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/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 597c3802b..22851606f 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()
@@ -500,6 +501,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
@@ -558,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
@@ -718,8 +722,29 @@ 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 from the original model.
+        exist_but_different = set()
+        other_model = ""
+        for reaction in filter(lambda x: not existing_filter(x), reaction_list):
+            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)} exist"
+                f" in both {self} and {other_model}, but have different "
+                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
+                else reaction.copy()
+                for reaction in filter(existing_filter, reaction_list)
+            ]
+        )
 
         context = get_context(self)
 
@@ -728,23 +753,17 @@ 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))
+            mets_to_add = [
+                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, stoich in reaction.metabolites.items()
+            }
+            reaction.metabolites = new_stoich
             reaction.update_genes_from_gpr()
 
         self.reactions += pruned
@@ -814,19 +833,15 @@ def remove_reactions(
 
                 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:
+                        met.remove_reaction(reaction)
+                        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))
+                        gene.remove_reaction(reaction)
 
-                        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))
@@ -869,7 +884,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
@@ -1257,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)
+                for met in rxn.metabolites:
+                    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 0623ae170..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
@@ -579,6 +580,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
+        ----------
+        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.
+        """
+        self._metabolites = value
+        context = get_context(self)
+        for met in value.keys():
+            met.add_reaction(self)
+
     @property
     def genes(self) -> FrozenSet:
         """Return the genes of the reaction.
@@ -837,11 +855,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.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.
@@ -930,12 +950,9 @@ def __setstate__(self, state: Dict) -> None:
             state["_gpr"] = GPR.from_string(state["_gpr"])
 
         self.__dict__.update(state)
-        for x in state["_metabolites"]:
+        for x in set(state["_metabolites"]).union(state["_genes"]):
             x._model = self._model
-            x._reaction.add(self)
-        for x in state["_genes"]:
-            x._model = self._model
-            x._reaction.add(self)
+            x.add_reaction(self)
 
     def copy(self) -> "Reaction":
         """Copy a reaction.
@@ -985,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
@@ -1018,7 +1035,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
@@ -1231,11 +1248,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 (
-                    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
@@ -1268,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
@@ -1287,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)
@@ -1487,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:
@@ -1499,7 +1517,7 @@ def _dissociate_gene(self, cobra_gene: Gene) -> None:
 
         """
         self._genes.discard(cobra_gene)
-        cobra_gene._reaction.discard(self)
+        cobra_gene.remove_reaction(self)
 
     def knock_out(self) -> None:
         """Knockout reaction by setting its bounds to zero."""
@@ -1656,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 39749fcd0..162e0361c 100644
--- a/src/cobra/core/species.py
+++ b/src/cobra/core/species.py
@@ -8,7 +8,7 @@
 
 
 if TYPE_CHECKING:
-    from .. import Model
+    from .. import Model, Reaction
 
 
 class Species(Object):
@@ -57,6 +57,32 @@ def reactions(self) -> FrozenSet:
         """
         return frozenset(self._reaction)
 
+    def add_reaction(self, reaction: "Reaction") -> None:
+        """Add reaction to .reaction field, with context.
+
+        If this is called within a context, will be reversed when exiting the context.
+
+        Parmeters
+        ---------
+        reaction: cobra.Reaction
+        """
+        self._reaction.add(reaction)
+
+    def remove_reaction(self, reaction: "Reaction") -> None:
+        """Remove reaction from .reaction field, with context.
+
+        If this is called within a context, will be reversed when exiting the context.
+
+        Parmeters
+        ---------
+        reaction: cobra.Reaction
+        """
+        self._reaction.remove(reaction)
+
+    def clear_reaction(self) -> None:
+        """Clear the reaction field."""
+        self._reaction.clear()
+
     def __getstate__(self) -> dict:
         """Return the state of the species.
 
diff --git a/tests/test_core/test_model.py b/tests/test_core/test_model.py
index e7e606ae1..3cf981bf4 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,61 @@ 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.