The determination of the Higgs boson trilinear self-coupling $\lambda_{3H}$ is a key goal of the LHC physics programme. Its precise measurement provides insight into the shape of the scalar potential and the mechanism of electroweak symmetry breaking. Higgs boson pair production in the process $gg \to HH$, particularly in the final state $HH \to b\bar{b}\gamma\gamma$, offers direct sensitivity to $\lambda_{3H}$. We present the first implementation of the Matrix Element Method at next-to-leading order (MEM@NLO) for this channel. Extending the method to NLO introduces significant methodological challenges, which we address through a new formalism integrated into the MoMEMta framework. Studies with simulated pseudo-datasets show that the method provides strong discriminating power and can extract the coupling modifier $\kappa_\lambda = \lambda_{3H} / \lambda^{\mathrm{SM}}_{3H}$ with high precision.