The stability of toroidal magnetic fields in radiative stellar interiors remains a major open issue to advance our understanding of the rotational and chemical evolution of low-mass stars. We perform 3D direct numerical simulations in a spherical geometry to examine the Tayler instability, a kink-type instability of purely toroidal fields expected to occur in stably stratified stellar interiors. The simulations are novel in that they consider a consistent background state derived from magnetohydrostatic equilibrium and explore the combined effect of gravity and thermal diffusion, as well as of fluid viscosity and magnetic resistivity. We trace the entire evolution of the instability from the linear to the nonlinear phase. Our simulations show that stable stratification and magnetic diffusivity can inhibit unstable modes, in agreement with linear stability analysis predictions by Bonanno & Urpin (2012). This suggests that toroidal fields in radiative stellar interiors may be only partially affected by Tayler instability and that the associated turbulent transport is lower than expected. These results may have implications for explaining the solid body rotation of the solar radiative core and the origin of the magnetic fields recently observed in red giant cores.