We study the dynamical evolution of the M87 globular cluster system (GCS) with a number of numerical simulations. We explore a range of different initial conditions for the GCS mass function (GCMF), for the GCS spatial distribution, and for the GCS velocity distribution. Our simulations include the effects of two-body relaxation, dynamical friction, and mass loss due to stellar evolution. We first confirm that an initial power-law GCMF such as that observed in young cluster systems can be readily transformed through dynamical processes into a bell-shaped GCMF. However, only models with initial velocity distributions characterized by a strong radial anisotropy increasing with the galactocentric distance are able to reproduce the observed constancy of the GCMF at all radii. We show that such strongly radial orbital distributions are inconsistent with the observed kinematics of the M87 GCS. The evolution of models with a bell-shaped GCMF with a turnover similar to that currently observed in old GCSs is also investigated. We show that models with this initial GCMF can satisfy all the observational constraints currently available on the GCS spatial distribution, the GCS velocity distribution, and on the GCMF properties. In particular, these models successfully reproduce both the lack of a radial gradient of the GCS mean mass recently found in an analysis of Hubble Space Telescope images of M87 at multiple locations and the observed kinematics of the M87 GCS. Our simulations also show that evolutionary processes significantly affect the initial GCS properties by leading to the disruption of many clusters and changing the masses of those that survive. The preferential disruption of inner clusters flattens the initial GCS number density profile, and it can explain the rising specific frequency with radius; we show that the inner flattening observed in the M87 GCS spatial distribution can be the result of the effects of dynamical evolution on an initially steep density profile.