Secondary neutrons are the predominant source of stray and leakage radiation outside the proton fields and therefore pose a risk to patients for the development of second cancers. The accuracy of the nuclear physics model used to predict stray neutron radiation fields in proton radiotherapy is not well understood. The general-purpose Monte Carlo N-Particle eXtended (MCNPX) code was used to calculate the therapeutic absorbed dose and neutron spectral fluence from a proton treatment unit using three nuclear physics models: the Bertini model, the cascade-exciton model (CEM), and the Li�ge cascade model (INCL4). The purpose of this study was to compare and quantify differences in predictions from these models for an un-modulated and a range modulated 160 MeV proton therapy beam in 1) the dosimetric characteristics of the therapeutic absorbed dose and 2) stray neutron fields produced by a proton radiotherapy unit using the default Bertini model as the baseline of comparison. The therapeutic absorbed dose was calculated in a water phantom 20 cm downstream of the nozzle exit, whereas the neutron spectral fluence was calculated in air. The ambient dose equivalent per therapeutic absorbed dose (H*(10)/D) from secondary neutrons produced by the nozzle components was calculated for each model. Based on these calculations, we determined H*(10)/D at the isocenter, 100 cm downstream from the isocenter, and at lateral distances of 100 cm from the isocenter. Our results indicated that calculations of the therapeutic absorbed dose ratios were in close agreement (~1% difference) for all three nuclear models. The H*(10)/D values differed only by 1 � 2 mSv Gy-1 at the isocenter with or without range modulation. However, the neutron spectral fluence calculations typically varied by about a factor of two or more using the two alternative models for intranuclear cascade processes.