diff --git a/festim/boundary_conditions/flux_bc.py b/festim/boundary_conditions/flux_bc.py
index 6d356eca6..4699c7391 100644
--- a/festim/boundary_conditions/flux_bc.py
+++ b/festim/boundary_conditions/flux_bc.py
@@ -141,6 +141,8 @@ class ParticleFluxBC(FluxBCBase):
             is applied
         value (float, fem.Constant, callable): the value of the particle flux
         species (festim.Species): the species to which the flux is applied
+        species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species}
+            where "name" is the variable name in the callable value and species is a festim.Species object.
 
     Attributes:
         subdomain (festim.SurfaceSubdomain): the surface subdomain where the particle flux
@@ -151,6 +153,8 @@ class ParticleFluxBC(FluxBCBase):
             fenics format
         bc_expr (fem.Expression): the expression of the particle flux that is used to
             update the value_fenics
+        species_dependent_value (dict): a dictionary containing the species that the value. Example: {"name": species}
+            where "name" is the variable name in the callable value and species is a festim.Species object.
 
 
     Usage:
@@ -160,11 +164,56 @@ class ParticleFluxBC(FluxBCBase):
         >>> ParticleFluxBC(subdomain=my_subdomain, value=lambda t: 1 + t, species="H")
         >>> ParticleFluxBC(subdomain=my_subdomain, value=lambda T: 1 + T, species="H")
         >>> ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + t, species="H")
+        >>> ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2, species="H", species_dependent_value={"c1": species1})
     """
 
-    def __init__(self, subdomain, value, species):
+    def __init__(self, subdomain, value, species, species_dependent_value={}):
         super().__init__(subdomain=subdomain, value=value)
         self.species = species
+        self.species_dependent_value = species_dependent_value
+
+    def create_value_fenics(self, mesh, temperature, t: fem.Constant):
+        """Creates the value of the boundary condition as a fenics object and sets it to
+        self.value_fenics.
+        If the value is a constant, it is converted to a fenics.Constant.
+        If the value is a function of t, it is converted to a fenics.Constant.
+        Otherwise, it is converted to a ufl Expression
+
+        Args:
+            mesh (dolfinx.mesh.Mesh) : the mesh
+            temperature (float): the temperature
+            t (dolfinx.fem.Constant): the time
+        """
+        x = ufl.SpatialCoordinate(mesh)
+
+        if isinstance(self.value, (int, float)):
+            self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.value)
+
+        elif callable(self.value):
+            arguments = self.value.__code__.co_varnames
+
+            if "t" in arguments and "x" not in arguments and "T" not in arguments:
+                # only t is an argument
+                if not isinstance(self.value(t=float(t)), (float, int)):
+                    raise ValueError(
+                        f"self.value should return a float or an int, not {type(self.value(t=float(t)))} "
+                    )
+                self.value_fenics = F.as_fenics_constant(
+                    mesh=mesh, value=self.value(t=float(t))
+                )
+            else:
+                kwargs = {}
+                if "t" in arguments:
+                    kwargs["t"] = t
+                if "x" in arguments:
+                    kwargs["x"] = x
+                if "T" in arguments:
+                    kwargs["T"] = temperature
+
+                for name, species in self.species_dependent_value.items():
+                    kwargs[name] = species.concentration
+
+                self.value_fenics = self.value(**kwargs)
 
 
 class HeatFluxBC(FluxBCBase):
diff --git a/test/test_flux_bc.py b/test/test_flux_bc.py
index fe2545fac..5aa05c8bc 100644
--- a/test/test_flux_bc.py
+++ b/test/test_flux_bc.py
@@ -88,6 +88,29 @@ def test_create_value_fenics_value(value, expected_value):
         assert np.isclose(bc.value_fenics.value, expected_value)
 
 
+def test_create_value_fenics_dependent_conc():
+    """Test that the value_fenics of ParticleFluxBC is set correctly when the value is dependent on the concentration"""
+    # BUILD
+    left = F.SurfaceSubdomain1D(1, x=0)
+    my_species = F.Species("test")
+    my_species.solution = F.as_fenics_constant(12, mesh)
+    T = F.as_fenics_constant(1, mesh)
+    t = F.as_fenics_constant(0, mesh)
+    bc = F.ParticleFluxBC(
+        subdomain=left,
+        value=lambda c: 1.0 + c,
+        species=my_species,
+        species_dependent_value={"c": my_species},
+    )
+
+    # RUN
+    bc.create_value_fenics(mesh, T, t)
+
+    # TEST
+    assert isinstance(bc.value_fenics, ufl.core.expr.Expr)
+    assert bc.value_fenics == 1.0 + my_species.solution
+
+
 def test_value_fenics_setter_error():
     left = F.SurfaceSubdomain1D(1, x=0)
     my_species = F.Species("test")
@@ -237,7 +260,7 @@ def test_create_value_fenics_type_HeatFluxBC(value, expected_type):
         (lambda t: 50.0 if t < 1 else 0.0, 50),
     ],
 )
-def test_create_value_fenics_value(value, expected_value):
+def test_create_value_fenics_value_HeatFluxBC(value, expected_value):
     """Test that"""
     # BUILD
     left = F.SurfaceSubdomain1D(1, x=0)