A novel first-principles self-consistent model that couples plasma and neutral atom physics suitable for the simulation of turbulent plasma behaviour in the tokamak edge region has been developed and implemented in the GBS code. While the plasma is modelled by the drift-reduced two fluid Braginskii equations, a kinetic model is used for the neutrals, valid in short and in long mean free path scenarios. The model includes ionization, charge-exchange, recombination, and elastic collisional processes. The model was used to study the transition form the sheath to the conduction limited regime, to include gas puffs in the simulations, and to investigate the interplay between neutral atoms and plasma turbulence.