We present an efficient integral equation approach to solve the heat equation in a two-dimensional, multiply connected domain, and with Dirichlet boundary conditions. Instead of using integral equations based on the heat kernel, we take the approach of discretizing in time, first. This leads to a non-homogeneous modified Helmholtz equation that is solved at each time step. The solution to this equation is formulated as a volume potential plus a double layer potential. The volume potential is evaluated using a fast multipole-accelerated solver. The boundary conditions are then satisfied by solving an integral equation for the homogeneous modified Helmholtz equation. The integral equation solver is also accelerated by the fast multipole method (FMM). For a total of N points in the discretization of the boundary and the domain, the total computational cost per time step is O(N) or O(N log N).