Numerical approximation of Turbulent flows poses significant computational challenges due to its high degrees of freedom and intricate statistical features\cite{fri}. In this work we introduce a novel approach to address these challenges by implementing a turbulence subgrid closure model for the Lattice Boltzmann Method (LBM) using Artificial Neural Networks (ANNs). The LBM is a broad class of computational methods, originating from the kinetic theory of gases, able to accurately model the dynamics of fluid flows at the mesoscopic level. The fluid flow is described by the dynamics of a set of discrete particle distribution functions (populations) following the stream and collide paradigm,  in which at each time step populations hop from lattice-site to lattice-site  and then incoming populations collide among one another. Leveraging on Direct Numerical Simulation (DNS) data, we learn a correction term for the BGK collisional operator as a physics-constrained Artificial Neural Network\cite{LBM-ML} taking as input the pre-collisional populations. Remarkably, our approach is fully local both in space and time, leading to reduced computational costs, enhanced generalization capabilities and interpretability. We show that the model offers superior performance compared to traditional methods, such as the Smagorinsky model, closely capturing both energy spectra and higher-order statistics of filtered DNS data. Moreover, we explore potential insights into the learned model's subgrid scale energy fluxes and possible connections with Multiple Relaxation Time BGK.