Extensive Molecular Dynamics simulations have been performed to investigate the structural and thermal properties of Si-doped fullerenes containing one and two silicon atoms. Both, a many-body potential and ab initio Density Functional Theory (DFT) have been used to investigate the structural features of the heterofullerenes. The competition between the exohedral and the substitutional types of doping, as a function of fullerene size (both even and odd heterofullerenes have been considered) and Si concentration, is analyzed. The DFT calculations confirm the main structural trends obtained with the many-body potential. The thermal stability and the structural transformations of the heterofullerenes have been also studied as a function of temperature (T=0–5000 K). The structural transformations include, local rearrangement of atoms, isomerization transitions, diffusion of atoms, eventual destruction of the cage, and sublimation of atoms. The isomerization transition between exohedral and substitutional isomers has been observed experimentally.