In this work, we study the asymptotic behaviour of solutions to the heat equation in exterior domains, i.e., domains which are the complement of a smooth compact set in $\mathbb{R}^N$. Different homogeneous boundary conditions are considered, including Dirichlet, Robin, and Neumann conditions for integrable initial data in $L^1(\Omega)$. After taking into account the loss of mass of the solution through the boundary, depending on the boundary conditions, we describe the asymptotic spatial distribution of the remaining mass in terms of the Gaussian and of a suitable asymptotic profile function. We show that our results have optimal time rates.