diff --git a/cara/models.py b/cara/models.py
index 521c40c2..8f2a8c42 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -644,20 +644,19 @@ class ConcentrationModel:
         k = (vg * 3600) / h
 
         return k + self.virus.decay_constant + self.ventilation.air_exchange(
-                                    self.room, self._next_state_change(time))
+                                                    self.room, time)
 
     @cached()
     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.
+        value reached by the concentration when time goes to infinity,
+        if all parameters were to stay time-independent.
         """
         V = self.room.volume
         IVRR = self.infectious_virus_removal_rate(time)
 
-        return (self.infected.emission_rate(self._next_state_change(time))) / (IVRR * V)
+        return (self.infected.emission_rate(time)) / (IVRR * V)
 
     @cached()
     def state_change_times(self):
@@ -694,19 +693,32 @@ class ConcentrationModel:
 
     @cached()
     def concentration(self, time: float) -> _VectorisedFloat:
-        # Note that time is not vectorised. You can only pass a single float
-        # to this method.
+        """
+        Virus quanta concentration, as a function of time.
+        The formulas used here assume that all parameters (ventilation,
+        emission rate) are constant between two state changes - only
+        the value of these parameters at the next state change, are used.
+
+        Note that time is not vectorised. You can only pass a single float
+        to this method.
+        """
 
         if time == 0:
             return 0.0
-        IVRR = self.infectious_virus_removal_rate(time)
+        try:
+            IVRR = self.infectious_virus_removal_rate(self._next_state_change(time))
+            concentration_limit = self._concentration_limit(self._next_state_change(time))
+        except ValueError:
+            raise ValueError("Concentration cannot be computed at a time"
+                             " larger than last state change (parameters"
+                             " are not defined)")
 
         t_last_state_change = self.last_state_change(time)
         concentration_at_last_state_change = self.concentration(t_last_state_change)
 
         delta_time = time - t_last_state_change
         fac = np.exp(-IVRR * delta_time)
-        return self._concentration_limit(time) * (1 - fac) + concentration_at_last_state_change * fac
+        return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
 
 
 @dataclass(frozen=True)