We apply mean-field theory and Hirsch-Fye quantum Monte Carlo method to study the spin-spin interaction in the bulk of three-dimensional topological insulators. We find that the spin-spin interaction has three different components: the longitudinal, the transverse and the transverse Dzyaloshinskii-Moriya-like terms. When the Fermi energy is located in the bulk gap of topological insulators, the spin-spin interaction decays exponentially due to Bloembergen-Rowland interaction. The longitudinal correlation is antiferromagnetic and the transverse correlations are ferromagnetic. When the chemical potential is in the conduction or valence band, the spin-spin interaction follows power law decay, and isotropic ferromagnetic interaction dominates in short separation limit.