  We build on the simplified spectral deferred corrections (SDC) coupling of
hydrodynamics and reactions to handle the case of nuclear statistical
equilibrium (NSE) and electron/positron captures/decays in the cores of massive
stars. Our approach blends a traditional reaction network on the grid with a
tabulated NSE state from a very large, $\mathcal{O}(100)$ nuclei, network. We
demonstrate how to achieve second-order accuracy in the simplified-SDC
framework when coupling NSE to hydrodynamics, with the ability to evolve the
star on the hydrodynamics timestep. We discuss the application of this method
to convection in massive stars leading up to core-collapse. We also show how to
initialize the initial convective state from a 1D model in a self-consistent
fashion. All of these developments are done in the publicly available Castro
simulation code and the entire simulation methodology is fully GPU accelerated.
