The theoretical description of the ultrafast optical response of nitrogen-vacancy centers has traditionally been developed using the nonlinear response function formalism based on a displaced harmonic oscillator model. However, restrictive defect geometries that lead to rich phonon mode distributions can limit the validity of this model. Here, we extend this framework by incorporating state-dependent oscillator modes characterized by different frequencies. We demonstrate that, as a consequence, the nonlinear response function is modified by a straightforward correction term that can be evaluated using operator algebraic techniques. These modifications directly impact two-dimensional correlation spectroscopy, which serves as a vital characterization tool for color centers by resolving ultrafast correlated electron-vibrational dynamics.