Doping Bi$_2$Se$_3$ by magnetic ions represents an interesting problem since it may break the time reversal symmetry needed to maintain the topological insulator character. Mn dopants in Bi$_2$Se$_3$ represent one of the most studied examples here. However, there is a lot of open questions regarding their magnetic ordering. In the experimental literature different Curie temperatures or no ferromagnetic order at all are reported for comparable Mn concentrations. This suggests that magnetic ordering phenomena are complex and highly susceptible to different growth parameters, which are known to affect material defect concentrations. So far theory focused on Mn dopants in one possible position, and neglected relaxation effects as well as native defects. We have used ab initio methods to calculate the Bi$_2$Se$_3$ electronic structure influenced by magnetic Mn dopants, and exchange interactions between them. We have considered two possible Mn positions, the substitutional and interstitial one, and also native defects. We have found a sizable relaxation of atoms around Mn, which affects significantly magnetic interactions. Surprisingly, very strong interactions correspond to a specific position of Mn atoms separated by van der Waals gap. Based on the calculated data we performed spin dynamics simulations to examine systematically the resulting magnetic order for various defect contents. We have found under which conditions the experimentally measured Curie temperatures ${T_{rm{C}}}$ can be reproduced, noticing that interstitial Mn atoms appear to be important here. Our theory predicts the change of ${T_{rm{C}}}$ with a shift of Fermi level, which opens the way to tune the system magnetic properties by selective doping.