In this paper, we derive a mathematical formulation of the higher order adjoint NEM-M2B2 equations by preconditioning the nodal interface neutron currents equations of the forward equations system, and by using the Lagrangian Multipliers analysis method. In the NEM-M2B2 system of equations, the quadratic transverse leakage approximation is used to model the leakage of neutrons between each node in the system. The solution of the adjoint equation can be used to perform adjoint-based predictive sensitivity/perturbation analysis. As an example, we use the mathematical adjoint solution as sensitivity weighting for predicting the response of the IAEA-3D benchmark’s eigenvalue to a perturbation in the independent parameters of the system (i.e., cross-sections). We also derive perturbation equations associated with the particular NEM-M2B2 model we are using. These perturbation-equations are used in predicting the model eigenvalue change without resorting to recalculating the forward NEM-M2B2 system of equations again (labeled as exact calculations). They also enabled construction of a reactivity sensitivity map showing the importance of each calculation node of the benchmark depending on its spatial and spectral coordinates. Perturbations were imposed on both the absorption cross-sections (fast and thermal) and the scattering cross-section of the IAEA-3D benchmark problem. Several verification steps were taken to ensure that the developed mathematical adjoint solver is adequate for adjoint analysis (e.g., commutativity checks, and comparison against exact calculations).