The focus of this work is on subgrid-scale models that leverage machine learning methods. In a conventional direct approach, training of the machine learning model attempts to reduce the mismatch between the subgrid-scale model and an accepted model (e.g., direct numerical simulation). We present an approach to expose the training process to the governing equations by using the adjoint equations. Specifically, we embed an untrained feedforward deep neural network into the governing equations and train by solving the adjoints of the governing equations, yielding end-to-end sensitivities of the flow variables with respect to the neural network parameters. This leverages the physics encoded in the governing equations, which provides greater robustness to the learned model. To illustrate, models based on the direct approach may result in unstable flow, while the embedded models are stable and capable of extrapolating on the out-of-sample scenarios that the model has not been exposed to in training [1, 2]. Moreover, the training target is distinct from the closure, thereby allowing training for any flow observable. The approach has been demonstrated for the simulation of isotropic turbulence [1] and free-shear-flow turbulence [2]; we extend it to the greater coupled-physics challenge of a turbulent reacting flow. Whereas standard large-eddy simulations of reacting flows model the turbulence and the combustion separately, our embedded approach considers them as a whole, leading to a model that takes into account their interaction. We demonstrate the approach on a planar premixed flame propagating through isotropic turbulence within the “thin reaction zone” regime, in which the well-known subgrid-scale turbulence-combustion coupling is a significant challenge to established models.