The Reynolds lubrication equation (RLE) is widely used to design sliding contacts in mechanical machinery. While providing an excellent description of hydrodynamic lubrication, friction in boundary lubrication regions is usually considered by empirical laws, because continuum theories are expected to fail for lubricant film heights h 0 ≪ 10 nm, especially in highly loaded tribosystems with normal pressures p n ≫ 0.1 GPa. Here, the performance of RLEs is validated by molecular dynamics simulations of pressurized (with p n = 0.2 to 1 GPa) hexadecane in a gold converging-diverging channel with minimum gap heights h 0 = 1.4 to 9.7 nm. For p n ≤ 0.4 GPa and h 0 ≥ 5 nm, agreement with the RLE requires accurate constitutive laws for pressure-dependent density and viscosity. An additional nonlinear wall slip law relating wall slip velocities to local shear stresses extends the RLE’s validity to even the most severe loading condition p n = 1 GPa and h 0 = 1.4 nm. Our results demonstrate an innovative route for continuum modeling of highly loaded tribological contacts under boundary lubrication.