Position- and direction-dependent diffusion coefficients of spherical, nonpolar solutes as represented by noble gases in a model lipid bilayer are investigated by nonequilibrium and equilibrium molecular dynamics (MD) simulations. The nonequilibrium MD method, which is based on linear response theory, is developed by determining an optimal range of external forces in which the response of solute mean velocities is linear and the perturbation to bilayer properties is minimal and by comparing with equilibrium MD simulation and experiment. Solute diffusivity is found to be smaller in a more ordered acyl-chain region and in the plane normal to the bilayer interface and exhibits markedly different dependences on solute size, depending on location in the bilayer interior and diffusion direction. The diffusion coefficients calculated closely follow those obtained from free-flight-length distributions of solute, a novel measure of free-volume anisotropy, based on the Einstein relation, suggesting that free volume is largely responsible for molecular diffusion in lipid bilayers. Theoretical formulas for solute diffusivity in lipid bilayers where the free-volume redistribution is governed primarily by overall rotation and local isomerization of Lipid molecules is proposed in the framework of dynamic free-volume theory. This theory predicts variable degrees of changes in molecular diffusivity with solute size and membrane structure depending on relative contributions of overall rotation and local isomerization of lipid molecules to free-volume redistribution and can explain the observed size dependences provided that lipid overall rotation (local isomerization) plays a major role in determining transverse (lateral) diffusivity of solute molecules.