The interaction between neutral particles and the plasma plays a key role in determining the dynamics of the tokamak boundary that, in turn, significantly impact the overall performance of the device. Leveraging the work in Wersal and Ricci [Nucl. Fusion 55, 123014 (2015)], the present thesis describes the development and implementation in the GBS code of a mass-conserving multi-component self-consistent model to simulate the interplay between neutrals and plasma in the tokamak boundary. The simulation results are shown to be a useful tool to disentangle the physics at play in the tokamak boundary, and in particular for the interpretation of gas puff imaging (GPI). Developed in the last decade, the GBS code [Ricci et al, Plasma Phys. Control. Fusion 54, 124047 (2012)] allows for self-consistent three-dimensional numerical simulations of the turbulent plasma and neutral dynamics in the tokamak boundary. In GBS, a set of Braginskii equations in the drift limit describes the plasma time evolution, and the neutrals are modelled by solving a kinetic equation using the method of characteristics. While GBS enables simulations in arbitrary magnetic configurations, here we focus on limited plasmas. We first describe the geometrical operators and proper boundary conditions to ensure that mass conservation is satisfied. In comparison to the previous non-mass-conserving model of GBS, the mass-conserving simulations capture more accurately the sharp transition of the plasma and neutral quantities between the edge and scrape-off layer regions. In addition, we show that mass conserving simulations allow for reliable quantitative studies of particle fluxes in the tokamak boundary. A multi-component model of the neutral-plasma interaction is then developed by extending the single-component model to the description of a deuterium plasma that includes electrons, D+ ions, D atoms, D2 molecules and D2+ ions. The molecular dynamics is introduced through a set of drift-reduced Braginskii equations for the D2+ species and considering a kinetic equation for D2 molecules, in addition to the kinetic equation for D atoms, thus resulting in a coupled system of kinetic equations for the atomic and molecular neutral distribution functions. The first multi-component GBS simulations show that, in the sheath-limited regime under consideration, most of the D2 molecules cross the last closed flux surface (LCFS) and are dissociated or ionized in the edge region, thus giving rise to sources of D atoms inside the LCFS. This leads to an inward radial shift of the peak of the plasma source due to ionization of D atoms with respect to the single-component simulations. The multi-component model is applied to the simulation of GPI diagnostics, where the presence of a molecular gas puff is simulated self-consistently with the plasma and neutral dynamics of the tokamak boundary. The injected molecules interact with the boundary plasma, resulting in the emission of light in the D alpha wavelength that can be measured to infer the turbulent properties of the plasma. The simulated mechanisms underlying the light emission, which include the excitation of D atoms and dissociation of both D2 and D2+, provide a reliable tool for the interpretation of GPI experimental measurements. The impact of neutral fluctuations on the D alpha emission rate is investigated, as well as the correlation between the D alpha emission and the plasma and neutral quantities.