A nonlinear gyrokinetic simulation model incorporating equilibrium current has been formulated for studying kinetic magnetohydrodynamic processes in magnetized plasmas. This complete formulation enables gyrokinetic simulation of both pressure-gradient-driven and current-driven instabilities as well as their nonlinear interactions in multi-scale simulations. The gyrokinetic simulation model recovers the ideal magnetohydrodynamic theory in the linear long wavelength regime including ideal and kinetic ballooning modes, kink modes and shear Alfvén waves. The implementation of this model in the global gyrokinetic particle code has been verified for the simulation of the effects of equilibrium current on the reversed shear Alfvén eigenmode in tokamaks.