A numerical simulation has been developed to study the formation in the intensities of ≲0.1-to~10.0-MeV protons of the long-lived energetic storm particle (ESP) events observed in association with the passage of flare-produced fast mode MHD shocks. Protons are traced numerically backward in time from a given observation point (i.e., particle detector), through their adiabatic scatter-free motion along the assumed laminar spiral interplanetary magnetic field (IMF), and finally through the shock discontinuity to their pre-shock interaction state. An ambient proton can undergo at most one shock encounter, during which, as seen in the shock frame, the proton is accelerated by its effective 'grad-B drift' along the induced electric field each time it crosses the shock. The time-reversed calculation reveals which of the observed protons interacted with the shock, where the shock interaction occurred, and by how much each interacting proton's kinetic energy was increased at the shock. Time profiles of the omnidirectional and unidirectional protons fluxes are formed from the trajectories by invoking Liouville's theorem and an ambient proton unidirectional flux j~T-&ggr;, where T is the kinetic energy and &ggr; the spectral exponent. Simulation predictions are shown for different observed proton energies under various assumptions concerning the shock speed and strength, the slope of the ambient proton energy spectrum, the variation of the angle &bgr;1 between the shock surface and the upstream IMF vector, and the heliocentric radial position r of the observation point. The simulation model correctly recovers the gross features of the intensity variations, flux anisotropies, and energy spectra time evolutions of many observed ESP events. Representative simulation results predict shock-induced enhancements in low-energy proton fluxes that are characterized by (1) long, steady preshock or upstream rises to peaks before the shock arrival and abrupt declines following the shock passage; (2) peak enhancement amplitudes that increase for lower-energy protons, stronger and faster shocks, softer ambient proton spectra, decreased &bgr;1, and increased r, (3) large, highly field-aligned preshock flux anisotropies directed directed away from the sun along the IMF, (4) small postshock flux anisotropies nearly perpendicular to the magnetic field and directed towards the sun along the IMF, and (5) steadily softening preshock energy that are steepest at the peak flux enhancement and slightly softer than the ambient spectra after the shock passage.