The response functions of four widely used Bonner sphere spectrometers (BSS) with an Lil scintillator and different 3 He detectors were calculated by means of the three-dimensional Monte Carlo neutron transport code, MCNP, taking into consideration a detailed description of the detector set-up; they were then compared with experimental calibration data. The calculated response functions reproduce the calibration data in shape. A common adjustment factor is only needed to take into account the neutron detection efficiency of the central detector. An interpolation method proposed to correct the BSS response functions to another polyethylene density was verified.