The presence of a particulate phase defines the effective response of a suspension to changes of the average heat flux. Using the random-point approximation we show that within the first order in the concentration, one needs to solve the problem for the temperature field created by a single inclusion in a matrix subject to an unsteady temperature gradient at infinity. We solve this problem by means of Laplace transform and use the solution as the first-order kernel in the functional expansion. From this kernel, we find the statistical average for the heat flux, which turns out to be a memory integral of the spatially averaged time-dependent temperature gradient. Thus, we discover that the constructive relationship between the average flux and averaged temperature gradient is not local in time, but rather involves a convolution integral that represents the memory due to the heterogeneity of the system. This is a novel result, which inter alia gives a rigorous justification to the usage of generalizations of the heat conduction law involving fractional time derivatives. The decay of the kernel is very close to t −1/2 for dimensionless times lesser than one and abruptly changes to t −3/2 for larger times.