Flip-flop processes due to magnetic dipole-dipole interaction between neighboring ions in crystals doped by rare-earth ions are one of the mechanisms of relaxation between hyperfine levels. Modeling of this mechanism has so far been macroscopic, characterized by an average rate describing the relaxation of all ions. Here, however, we present a microscopic model of flip-flop interactions between individual nuclear spins of dopant ions. Every ion is situated in a unique local environment in the crystal, where each ion has different distances and a unique orientation relative to its nearest neighbors, as determined by the lattice structure. Thus each ion has a unique flip-flop rate and the collective relaxation dynamics of all ions in a bulk crystal is a sum of many exponential decays, giving rise to a distribution of rates rather than a single average decay rate. We employ this model to calculate flip-flop rates in ${\mathrm{Pr}}^{3+}$:${\mathrm{Y}}_{2}{\mathrm{SiO}}_{5}$ and show experimental measurements of population decay of the ground state hyperfine levels at $\ensuremath{\sim}2$ K. We also present a method to measure rates of individual transitions from hole burning spectra that requires significantly fewer fitting parameters in theoretical rate equations compared to earlier work. Furthermore, we measure the effect of external magnetic field on the flip-flop rates and observe that the rates slow down by two orders of magnitude in a field of 5--10 mT.