We recently reported the equilibrium maximum solubility of cholesterol in a lipid bilayer, χ(chol)/* to be 0.66 in four different phosphatidylcholines, and 0.51 in a phosphatidylethanolamine (Huang, J., J. T. Buboltz, and G. W. Feigenson, 1999. Biochim. Biophys. Acta. in press). Here we present a model of cholesterol-phospholipid mixing that explains these observed values of χ(chol)/*. Monte Carlo simulations show that pairwise-additivity of nearest-neighbor interactions is inadequate to describe all the χ(chol)/* values. Instead, if cholesterol multibody interactions are assigned highly unfavorable energy, then jumps occur in cholesterol chemical potential that lead to its precipitation from the bilayer. Cholesterol precipitation is most likely to occur near three discrete values of cholesterol mole fraction, 0.50, 0.57, and 0.67, which correspond to cholesterol/phospholipid mole ratios of 1/1, 4/3, and 2/1, respectively. At these solubility limits, where cholesterol chemical potential jumps, the cholesterol-phospholipid bilayer mixture forms highly regular lipid distributions in order to minimize cholesterol-cholesterol contacts. This treatment shows that dramatic structural and thermodynamic changes can occur at particular cholesterol mole fractions without any stoichiometric complex formation. The physical origin of the unfavorable cholesterol multibody interaction is explained by an 'umbrella model': in a bilayer, nonpolar cholesterol relies on polar phospholipid headgroup coverage to avoid the unfavorable free energy of cholesterol contact with water. Thus, at high cholesterol mole fraction, this unfavorable free energy, not any favorable cholesterol-phospholipid interaction, dominates the mixing behavior. This physical origin also explains the 'cholesterol condensing effect' and the increase in acyl chain order parameter in cholesterol- phospholipid mixtures.