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 h0 ≪ 10 nm, especially in highly loaded tribosystems with normal pressures pn ≫ 0.1 GPa. Here, the performance of RLEs is validated by molecular dynamics simulations of pressurized (with pn = 0.2 to 1 GPa) hexadecane in a gold converging-diverging channel with minimum gap heights h0 = 1.4 to 9.7 nm. For pn ≤ 0.4 GPa and h0 ≥ 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 pn = 1 GPa and h0 = 1.4 nm. Our results demonstrate an innovative route for continuum modeling of highly loaded tribological contacts under boundary lubrication.