We have investigated the origin of the S1-T1 energy levels inversion for heptazine, and other N-doped π-conjugated hydrocarbons, leading thus to an unusually negative singlet-triplet energy gap ((Formula presented.)). Since this inversion might rely on substantial doubly-excited configurations to the S1 and/or T1 wavefunctions, we have systematically applied multi-configurational SA-CASSCF and SC-NEVPT2 methods, SCS-corrected CC2 and ADC(2) approaches, and linear-response TD-DFT, to analyze if the latter method could also face this challenging issue. We have also extended the study to B-doped π-conjugated systems, to see the effect of chemical composition on the results. For all the systems studied, an intricate interplay between the singlet-triplet exchange interaction, the influence of doubly-excited configurations, and the impact of dynamic correlation effects, serves to explain the (Formula presented.) values found for most of the compounds, which is not predicted by TD-DFT.