The overall performance of fusion devices, such as tokamaks, is strongly correlated to the phenomena that occur in the boundary region, the outer plasma region that faces the wall of the device. The boundary plays a crucial role in regulating the heat and particle fluxes on the machine wall, impurity concentration, helium ash removal, and neutral dynamics. Modeling the turbulent plasma dynamics in the boundary is challenging because of its wide range of plasma parameters, particularly the plasma collisionality. Indeed, the presence of regions of low collisionality can limit the use of fluid models, such that a gyrokinetic (GK) description is necessary. At the same time, the use of first-principles GK simulations is challenged by their high computational cost, the complexity of implementing collision operators, and the ordering parameters of the boundary. The present thesis aims to develop a GK boundary model to simulate plasma dynamics in the boundary region efficiently. By considering a proper ordering for the boundary region, a GK model is first derived to describe electromagnetic fluctuations and collisional effects using a first-principles GK Coulomb (Fokker-Planck) collision operator. Then, with the goal of simulating the boundary efficiently, the full gyrocenter (full-F) distribution function is projected onto a Hermite and Laguerre velocity-space polynomial basis, which yields a hierarchy of fluid-like equations for the expansion coefficients, referred to as the gyro-moments (GMs). The GM hierarchy equation is valid at an arbitrary level of collisionality. To assess the properties of the model developed and the role of collisional effects, a linearized version of the GK model, which includes a first-principles GK Coulomb collision operator, is derived and numerically implemented for the first time, removing the use of ad-hoc approximations often considered by previous collision operator models. The GM approach is also applied to other advanced collision operator models. A comprehensive comparison of collision operators is performed and shows that simplified collision operator models deviate significantly from the GK Coulomb operator in predicting the linear growth rates of microinstabilities and their nonlinear saturation. These results question the accuracy of previously-used simplified collision operators and their applications to the boundary region. The ability of the GM approach to describe the high and low collisionality regimes is also assessed. Supported by analytical results and comparisons with well-established models and numerical codes, the fluid predictions are recovered with a low number of GMs. On the other hand, increasing the number of GMs allows the description of microinstabilities in the collisionless limit. Finally, the first nonlinear and full-F simulations of the GM approach are presented for the cases of one-dimensional electrostatic plasma waves and a linear plasma device. This thesis demonstrates the potential of the GM approach to enable comprehensive simulations of the boundary with a lower computational cost than GK simulations, particularly at high collisionality. At the same time, the GM approach provides an improved fluid description over the drift-reduced Braginskii-like fluid model in the low collisionality limit.