diff --git a/cara/models.py b/cara/models.py
index 609da9a3..521c40c2 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -647,7 +647,13 @@ class ConcentrationModel:
                                     self.room, self._next_state_change(time))
 
     @cached()
-    def concentration_limit(self, time: float) -> _VectorisedFloat:
+    def _concentration_limit(self, time: float) -> _VectorisedFloat:
+        """
+        Provides a constant that represents the theoretical asymptotic 
+        value reached by the concentration when time goes to infinity.
+        It's a piecewise constant function - at any time it takes the
+        value it would have at the end of the current interval.
+        """
         V = self.room.volume
         IVRR = self.infectious_virus_removal_rate(time)
 
@@ -700,8 +706,7 @@ class ConcentrationModel:
 
         delta_time = time - t_last_state_change
         fac = np.exp(-IVRR * delta_time)
-        concentration_limit = self.concentration_limit(time)
-        return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
+        return self._concentration_limit(time) * (1 - fac) + concentration_at_last_state_change * fac
 
 
 @dataclass(frozen=True)
diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py
index 10cf6a8c..f14735d8 100644
--- a/cara/tests/models/test_concentration_model.py
+++ b/cara/tests/models/test_concentration_model.py
@@ -156,7 +156,7 @@ def test_concentration_model_constant_parameters():
         )
     )
     times = [0.1, 10, 20, 24]
-    concentration_limits = np.array([c_model.concentration_limit(t)
+    concentration_limits = np.array([c_model._concentration_limit(t)
                                      for t in times])
     IVRRs = np.array([c_model.infectious_virus_removal_rate(t)
                       for t in times])