The stability of boron and carbon dopants (D) at vacant octahedral sites in rutile has been investigated using density functional theory quantum mechanical (QM) modeling. Three different types of sites were considered: vacant Ti lattice site (Dvac), vacant Ti site + Ti interstitial (DFrenkel), and interstitial site (Dint). The defect formation energies, Ed, at different temperatures and gas partial pressures were calculated from the QM total energies of the relaxed defect and host structures, combined with chemical potentials of the reservoir components that were corrected for both temperature and gas partial pressure variations. The contribution of vibrational free energy to Ed as a function of temperature was evaluated for the BFrenkel model using ab initio phonon density of states calculations. The calculated vibrational and configurational free energy terms were of opposite sign and partially cancelled, giving a relatively small (0.3 eV at 700 K) combined contribution to Ed. Under strongly reducing conditions at 1500 K, boron incorporation at interstitial and Frenkel sites is favored, with Ed values of 0.53 and 0.8 eV respectively, whereas the Ed values for carbon doping were high (>4 eV) for all three models under high-temperature reducing conditions. Under oxygen-rich conditions relevant to sol-gel processing, the Dvac model was favored for both boron and carbon. Further stabilization of the Dvac model for boron was obtained at a protonated vacancy site, giving Ed = 1.16 eV for (B+H)vac at 700 K.