A first-principles self-consistent model that couples plasma and neutral physics suitable for the simulation of turbulent plasma behavior in the tokamak SOL is presented. While the plasma is modeled by the drift-reduced two fluid Braginskii equations, a kinetic model for the neutrals is developed, valid in short and in long mean free path scenarios. The model includes ionization, charge-exchange, recombination, and elastic collisional processes. The solution of the neutral kinetic equation is implemented within the GBS plasma turbulence code (Ricci et al 2012 Plasma Phys. Control. Fusion54 124047) and it is performed by using the method of characteristics. The details of the numerical implementation are discussed. Finally, we show initial results of the first self-consistent simulations of plasma turbulence and neutral dynamics.