Understanding the physical mechanisms at play in the interaction between turbulent plasma and neutral particles is a crucial issue that we approach in this Thesis by using a first-principles self-consistent model of the tokamak periphery implemented in the GBS code. 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 neutral kinetic equation is solved by using the method of characteristics. We identify the key elements determining the interaction between neutrals and the turbulent plasma focusing on a tokamak with a toroidal rail limiter on the high-field side equatorial midplane. For this purpose, we simulate the dynamics of the plasma and the neutrals in a domain that includes both the confined edge region and the scrape-off layer (SOL). It turns out that, in the considered plasma conditions, neither the fluctuations of the neutral moments, nor the friction between neutrals and the plasma impact the time-averaged plasma profiles significantly. Thanks to this study, we derive a simple model for the neutral-plasma interaction, which is helpful to identify and understand the principal physical processes at play in the tokamak periphery. By studying the dynamics of the neutral-plasma interplay along the magnetic field lines in the SOL, we derive a refined two-point model from the drift-reduced Braginskii equations that balances the parallel and perpendicular transport of plasma and heat, and takes into account the plasma-neutral interaction. The model estimates the electron temperature drop along a field line, from a region far from the limiter to the limiter plates. The refined two-point model is shown to be in very good agreement with the simulation results. Finally, we self-consistently simulate a diagnostic neutral gas puff, which is often used experimentally as a tool to learn about the turbulence properties in the tokamak periphery. In particular, we investigate the impact of neutral density fluctuations on the D-α light emission, finding that at a radial distance from the gas puff smaller than the neutral mean free path, neutral density fluctuations are anti-correlated with plasma density, electron temperature, and D-α fluctuations, while at distances from the gas puff larger than the neutral mean free path, a non-local shadowing effect influences the neutrals, and the D-α fluctuations are correlated with the neutral density fluctuations.