The Kirkwood-Buff solution theory to give the relations between macroscopic thermodynamic properties and the fluctuation integrals (G_ij) was utilized to predict solubility of solutes in supercritical fluids. The solvent-solute fluctuation integral (G₂₁) in the derivation for solubility of solute is expressed in terms of the solvent-solvent fluctuation integral (G₁₁) using the hard sphere expansion (HSE) conformal solution method with the modified mean density approximation (MMDA) where the scaling factor (R₁₂) represents the ratio of the first peak heights of the radial distribution functions for the mixture and the reference fluid having the mean density determined from the mean density approximation (MDA). The values of R₁₂ were evaluated by considering it as an adjustable parameter and solving the Ornstein-Zernike equation with the hypernetted chain (HNC) closure, and were compared. It is shown that solubility can be evaluated with an equation of state for pure supercritical fluid, three molecular parameters, and the scaling factor (R₁₂) without knowledge of critical properties of solutes, which can not be measured precisely for some organic solids. This model based on the molecular theory leads to better results in solubility calculations than both the Peng-Robinson equation of state with the classical mixing rule and the previous method with the original MDA instead of the MMDA. It might be due to the superiority of the MMDA over the original MDA.