Nebular line emission diagnostics have been instrumental in studying the galaxy properties over the years. As we move into the era of next generation telescopes giving us a wider access to the high redshift universe the applicability of these diagnostics are put into question. We are working on developing a self consistent model of nebular emission using state of the art cosmological hydrodynamic galaxy formation simulation (SIMBA) and photoionization code (CLOUDY). We have adopted a novel approach where we have developed the first-ever model for calculating the nebular emission line spectrum on a particle-by-particle basis for cosmological galaxy formation simulations which is a significant step forward from more commonly-used lookup-table methods. We have already developed a basic framework and shown in a pilot study that it can reproduce the observed SDSS-DR8 BPT curve. In this talk, I will be going over the modeling details and results from the pilot study as well as some preliminary diagnostic calibration for high redshift galaxies.