diff --git a/festim/__init__.py b/festim/__init__.py
index 2af946b92..853092014 100644
--- a/festim/__init__.py
+++ b/festim/__init__.py
@@ -28,10 +28,12 @@
 
 from .hydrogen_transport_problem import HydrogenTransportProblem
 
-from .species import Species, Trap
+from .species import Species, Trap, ImplicitSpecies
 
 from .subdomain.surface_subdomain import SurfaceSubdomain1D
 from .subdomain.volume_subdomain import VolumeSubdomain1D
 
 from .exports.vtx import VTXExport
 from .exports.xdmf import XDMFExport
+
+from .reaction import Reaction
diff --git a/festim/reaction.py b/festim/reaction.py
new file mode 100644
index 000000000..104cebbb7
--- /dev/null
+++ b/festim/reaction.py
@@ -0,0 +1,84 @@
+import festim as F
+from typing import Union
+
+from ufl import exp
+
+
+class Reaction:
+    """A reaction between two species, with a forward and backward rate.
+
+    Arguments:
+        reactant1 (Union[F.Species, F.ImplicitSpecies]): The first reactant.
+        reactant2 (Union[F.Species, F.ImplicitSpecies]): The second reactant.
+        product (F.Species): The product.
+        k_0 (float): The forward rate constant pre-exponential factor.
+        E_k (float): The forward rate constant activation energy.
+        p_0 (float): The backward rate constant pre-exponential factor.
+        E_p (float): The backward rate constant activation energy.
+
+    Attributes:
+        reactant1 (Union[F.Species, F.ImplicitSpecies]): The first reactant.
+        reactant2 (Union[F.Species, F.ImplicitSpecies]): The second reactant.
+        product (F.Species): The product.
+        k_0 (float): The forward rate constant pre-exponential factor.
+        E_k (float): The forward rate constant activation energy.
+        p_0 (float): The backward rate constant pre-exponential factor.
+        E_p (float): The backward rate constant activation energy.
+
+    Usage:
+        >>> # create two species
+        >>> species1 = F.Species("A")
+        >>> species2 = F.Species("B")
+
+        >>> # create a product species
+        >>> product = F.Species("C")
+
+        >>> # create a reaction between the two species
+        >>> reaction = Reaction(species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3)
+        >>> print(reaction)
+        A + B <--> C
+
+        >>> # compute the reaction term at a given temperature
+        >>> temperature = 300.0
+        >>> reaction_term = reaction.reaction_term(temperature)
+
+    """
+
+    def __init__(
+        self,
+        reactant1: Union[F.Species, F.ImplicitSpecies],
+        reactant2: Union[F.Species, F.ImplicitSpecies],
+        product: F.Species,
+        k_0: float,
+        E_k: float,
+        p_0: float,
+        E_p: float,
+    ) -> None:
+        self.reactant1 = reactant1
+        self.reactant2 = reactant2
+        self.product = product
+        self.k_0 = k_0
+        self.E_k = E_k
+        self.p_0 = p_0
+        self.E_p = E_p
+
+    def __repr__(self) -> str:
+        return f"Reaction({self.reactant1} + {self.reactant2} <--> {self.product}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})"
+
+    def __str__(self) -> str:
+        return f"{self.reactant1} + {self.reactant2} <--> {self.product}"
+
+    def reaction_term(self, temperature):
+        """Compute the reaction term at a given temperature.
+
+        Arguments:
+            temperature (): The temperature at which the reaction term is computed.
+        """
+        k = self.k_0 * exp(-self.E_k / (F.k_B * temperature))
+        p = self.p_0 * exp(-self.E_p / (F.k_B * temperature))
+
+        c_A = self.reactant1.concentration
+        c_B = self.reactant2.concentration
+
+        c_C = self.product.solution
+        return k * c_A * c_B - p * c_C
diff --git a/festim/species.py b/festim/species.py
index eda76f2e1..a3f45eacd 100644
--- a/festim/species.py
+++ b/festim/species.py
@@ -1,3 +1,6 @@
+from typing import List
+
+
 class Species:
     """
     Hydrogen species class for H transport simulation.
@@ -10,6 +13,7 @@ class Species:
         solution (dolfinx.fem.Function or ...): the solution for the current timestep
         prev_solution (dolfinx.fem.Function or ...): the solution for the previous timestep
         test_function (ufl.Argument or ...): the testfunction associated with this species
+        concentration (dolfinx.fem.Function): the concentration of the species
 
     Usage:
         >>> from festim import Species, HTransportProblem
@@ -28,6 +32,16 @@ def __init__(self, name: str = None) -> None:
         self.prev_solution = None
         self.test_function = None
 
+    def __repr__(self) -> str:
+        return f"Species({self.name})"
+
+    def __str__(self) -> str:
+        return f"{self.name}"
+
+    @property
+    def concentration(self):
+        return self.solution
+
 
 class Trap(Species):
     """Trap species class for H transport simulation.
@@ -51,3 +65,49 @@ class Trap(Species):
 
     def __init__(self, name: str = None) -> None:
         super().__init__(name)
+
+
+class ImplicitSpecies:
+    """Implicit species class for H transport simulation.
+    c = n - others
+
+    Args:
+        n (float): the total concentration of the species
+        others (List[Species]): the list of species from which the implicit
+            species concentration is computed (c = n - others)
+        name (str, optional): a name given to the species. Defaults to None.
+
+    Attributes:
+        name (str): a name given to the species.
+        n (float): the total concentration of the species
+        others (List[Species]): the list of species from which the implicit
+            species concentration is computed (c = n - others)
+        concentration (form): the concentration of the species
+
+    """
+
+    def __init__(
+        self,
+        n: float,
+        others: List[Species] = None,
+        name: str = None,
+    ) -> None:
+        self.name = name
+        self.n = n
+        self.others = others
+
+    def __repr__(self) -> str:
+        return f"ImplicitSpecies({self.name}, {self.n}, {self.others})"
+
+    def __str__(self) -> str:
+        return f"{self.name}"
+
+    @property
+    def concentration(self):
+        if len(self.others) > 0:
+            for other in self.others:
+                if other.solution is None:
+                    raise ValueError(
+                        f"Cannot compute concentration of {self.name} because {other.name} has no solution"
+                    )
+        return self.n - sum([other.solution for other in self.others])
diff --git a/test/test_reaction.py b/test/test_reaction.py
new file mode 100644
index 000000000..58fcfbcb7
--- /dev/null
+++ b/test/test_reaction.py
@@ -0,0 +1,99 @@
+import pytest
+import festim as F
+from dolfinx.fem import FunctionSpace, Function
+from dolfinx.mesh import create_unit_cube
+from mpi4py import MPI
+from ufl import exp
+
+
+def test_reaction_init():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # create a reaction between the two species
+    reaction = F.Reaction(
+        species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3
+    )
+
+    # check that the attributes are set correctly
+    assert reaction.reactant1 == species1
+    assert reaction.reactant2 == species2
+    assert reaction.product == product
+    assert reaction.k_0 == 1.0
+    assert reaction.E_k == 0.2
+    assert reaction.p_0 == 0.1
+    assert reaction.E_p == 0.3
+
+
+def test_reaction_repr():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # create a reaction between the two species
+    reaction = F.Reaction(
+        species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3
+    )
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = "Reaction(A + B <--> C, 1.0, 0.2, 0.1, 0.3)"
+    assert repr(reaction) == expected_repr
+
+
+def test_reaction_str():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # create a reaction between the two species
+    reaction = F.Reaction(
+        species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3
+    )
+
+    # check that the __str__ method returns the expected string
+    expected_str = "A + B <--> C"
+    assert str(reaction) == expected_str
+
+
+@pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0])
+def test_reaction_reaction_term(temperature):
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # create a reaction between the two species
+    reaction = F.Reaction(
+        species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3
+    )
+
+    # set the concentrations of the species
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+    species2.solution = Function(V)
+    product.solution = Function(V)
+
+    # test the reaction term at a given temperature
+    def arrhenius(pre, act, T):
+        return pre * exp(-act / (F.k_B * T))
+
+    k = arrhenius(reaction.k_0, reaction.E_k, temperature)
+    p = arrhenius(reaction.p_0, reaction.E_p, temperature)
+    expected_reaction_term = (
+        k * species1.solution * species2.solution - p * product.solution
+    )
+
+    assert reaction.reaction_term(temperature) == expected_reaction_term
diff --git a/test/test_species.py b/test/test_species.py
index d4f4f3b6b..2eafd9844 100644
--- a/test/test_species.py
+++ b/test/test_species.py
@@ -2,6 +2,10 @@
 import dolfinx
 import ufl
 import numpy as np
+import pytest
+from dolfinx.fem import FunctionSpace, Function
+from dolfinx.mesh import create_unit_cube
+from mpi4py import MPI
 
 
 def test_assign_functions_to_species():
@@ -29,3 +33,85 @@ def test_assign_functions_to_species():
         assert isinstance(spe.solution, dolfinx.fem.Function)
         assert isinstance(spe.prev_solution, dolfinx.fem.Function)
         assert isinstance(spe.test_function, ufl.Argument)
+
+
+def test_species_repr_and_str():
+    """Test that the __repr__ and __str__ methods of the Species class returns the
+    expected string.
+    """
+    # create a species
+    species = F.Species("A")
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = "Species(A)"
+    assert repr(species) == expected_repr
+
+    # check that the __str__ method returns the expected string
+    expected_str = "A"
+    assert str(species) == expected_str
+
+
+def test_implicit_species_repr_and_str():
+    """Test that the __repr__ and __str__ methods of the ImplicitSpecies class
+    returns the expected string.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = f"ImplicitSpecies(C, 3.0, {[species1, species2]})"
+    assert repr(implicit_species) == expected_repr
+
+    # check that the __str__ method returns the expected string
+    expected_str = "C"
+    assert str(implicit_species) == expected_str
+
+
+def test_implicit_species_concentration():
+    """Test that the concentration of an implicit species is computed
+    correctly.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # set the solutions of the two species
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+    species2.solution = Function(V)
+
+    # test the concentration of the implicit species
+    expected_concentration = implicit_species.n - (
+        species1.solution + species2.solution
+    )
+    assert implicit_species.concentration == expected_concentration
+
+
+def test_implicit_species_concentration_with_no_solution():
+    """Test that a ValueError is raised when on of the 'others' species
+    has no solution and the concentration of the implicit species is
+    requested.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # set the solution of the first species
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+
+    # test that a ValueError is raised when the second species has no solution
+    with pytest.raises(ValueError):
+        implicit_species.concentration