Gearboxes are one of critical components of railroad vehicles, aero-engines and other rotating machinery, and their health condition is crucial for the safe operation of these rotating machines. Fault signature extraction is a vital step for gearbox health condition monitoring. To this end, knowing the characteristics of fault signature is a prerequisite. Dynamic model is an effective tool for analysing the fault features of gearbox. However, in reality, the vibrations measured from the gearbox housing are usually coupled with the vibrations of the rolling bearings, the gear mesh and the gearbox housing, which are influenced by the friction, lubrication, surface roughness, etc., of each contact pair. Consequently, any dynamic model of individual components cannot accurately reflect the actual vibration behavior of the gearbox and also the resulted vibration response cannot match the vibrations measured from the gearbox housing. To resolve this problem, this paper develops a comprehensive dynamic model to analyze the vibration response of a gearbox with multiple localized defects by coupling the effects of spur gear mesh, deep groove ball bearings, time-varying displacement excitation of the bearing caused by localized defects, time-varying mesh stiffness (TVMS) of the gear, and mixed elastohydrodynamic (EHD) lubrication. The coupled vibrations of the gearbox housing can then be obtained to analyze the vibration behavior of a gearbox with multiple localized defects. In addition, the influences of the relative motion of the bearing rollers and flexible cage on the gearbox vibration are analyzed. At last, the proposed dynamic model is validated by experiments.