How is the volume kept constant in an NVE MD simulation?
The answer to this is linked to how you define the “volume” of your system, and the boundary conditions. The simplest case is simply to consider an isolated systems, like a molecule in gas phase. This corresponds to $N$ atoms (or molecules or particles) with no restrictions of their coordinates and no boundary conditions. This is essentially an infinite box. Its volume naturally does not and cannot vary.

Note that the volume of the system is not defined as “the volume actually occupied by the atoms”, such as for example the Connolly volume of your system. Instead, the volume of the system is the boundary of space where the atoms can go.

If you take a finite volume for your simulation, for example of cube of size $L$ and volume $V=L^3$. You have to put periodic boundary conditions, for example reflexive (your box is surrounded by a wall, and particles bounce on it) or periodic (look up for periodic boundary conditions in a textbook). Then from the initial positions and velocities of the atoms, Newton's equation allows you to propagate those quantities and determine their evolution in time. But nothing in Newton's law affects the system size $L$ (and thus $V=L^3$) or the boundaries. They will stay fixed and the volume will be constant.

Also, in all cases, because of conservation of matter, $N$ is constant. And if the forces acting on the system are conservative (as are all interatomic forces), then the total energy in the system $E$ is conserved.* Thus, the system samples the $(N, V, E)$ thermodynamic ensemble.



* the momentum of the system is also conserved… and it is not as trivial as it may sound :)