Numerical simulation of gravitational N-body systems is an important tool for studying the dynamic behaviour of stellar systems, and in some cases is the only option available given the extremely large time scales involved. The direct summation approach, which evaluates the force between each pair of particles at each time step, produces the most accurate results. However despite many algorithmic advances this method remains a computationally challenging problem owing to its O(N²) scaling characteristics. The desire to model increasingly larger systems has spurred the adoption of parallel computation techniques, but unfortunately many of the strategies used to accelerate sequential direct N-body simulations hinder their efficient parallelization. This paper investigates the use of parallel discrete event simulation as an alternative to the usual iterative time-stepping approach. By decomposing typical operations into finer-grained events, it is shown that there exists considerable potential for exploiting the model's inherent concurrency. In addition, it is demonstrated how certain optimizations that are normally difficult to parallelize are incorporated naturally into the parallel discrete event paradigm.