Low energy protons (< 100 − 300 keV) in the Van Allen belt and the outer regions can enter the field of view of X-ray focusing telescopes, interact with the Wolter-I optics, and reach the focal plane. The funneling of soft protons was discovered after the damaging of the Chandra/ACIS Front-Illuminated CCDs in September 1999 after the first passages through the radiation belt. The use of special filters protects the XMM-Newton focal plane below an altitude of 70000 km, but above this limit the effect of soft protons is still present in the form of sudden flares in the count rate of the EPIC instruments that can last from hundreds of seconds to hours and can hardly be disentangled from X-ray photons, causing the loss of large amounts of observing time. The accurate characterization of (i) the distribution of the soft proton population, (ii) the physics interaction at play, and (iii) the effect on the focal plane, are mandatory to evaluate the background and design the proton magnetic diverter on board future X-ray focusing telescopes (e.g. ATHENA). Several solutions have been proposed so far for the primary population and the physics interaction, however the difficulty in precise angle and energy measurements in laboratory makes the smoking gun still unclear. Since the only real data available is the XMM-Newton spectrum of soft proton flares in orbit, we try to characterize the input proton population and the physics interaction by simulating, using the BoGEMMS framework, the proton interaction with a simplified model of the X-ray mirror module and the focal plane, and comparing the result with a real observation. The analysis of ten orbits of observations of the EPIC/pn instrument show that the detection of flares in regions far outside the radiation belt is largely influenced by the different orientation of the Earth's magnetosphere respect with XMM-Newton's orbit, confirming the solar origin of the soft proton population. The Equator-S proton spectrum at 70000 km altitude is used for the proton population entering the optics, where a combined multiple and Firsov scattering is used as physics interaction. If the thick filter is used, the soft protons in the 30-70 keV energy range are the main contributors to the simulated spectrum below 10 keV. We are able to reproduce the proton vignetting observed in real data-sets, with a ∼ 50% decrease from the inner to the outer region, but a maximum flux of ∼ 0.01 counts cm −2 s −1 keV −1 is obtained below 10 keV, about 5 times lower than the EPIC/MOS detection and 100 times lower than the EPIC/pn one. Given the high variability of the flare intensity, we conclude that an average spectrum, based on the analysis of a full season of soft proton events is required to compare Monte Carlo simulations with real events.