First ditch all the object orientation and encapsulation and stuff. Your data is a 2xN Numpy array. Your visualization is a scatter plot in Matplotlib. Voila, 80% of the code is gone.
For the position updates, you either use a repulsive potential to approximate the hard spheres and do molecular dynamics, showing how to integrate Newton's second law and the Verlet scheme and ergodicity and the whole shebang. Or you do Monte Carlo for the positional updates and keep the exact hard spheres. You discuss statistical mechanics concepts like ensembles and thermostats and stuff.
Then you produce results like the pair correlation function and compare it with the Carnahan-Starling equation, dig into the really cool stuff. Compute velocity autocorrelation functions, test what happens when you change density and temperature, talk about phase diagrams, etc.
This is actually an amazingly deep subject, yet very accessible and intuitive, that sits on the border between physics and chemistry. Sad to see it treated like this. Would suggest that people have a look at the book by Allen and Tildesley which is much much better. They have both Python and Fortran example code on Github.