Selective ion binding by human lysozyme and its mutants is probed with the three-dimensional interaction site model theory which is the statistical mechanical integral equation theory. Preliminary and partial results of the study have been already published (Yoshida, N. et al. J. Am. Chem. Soc. 2006, 128, 12042-12043). The calculation was carried out for aqueous solutions of three different electrolytes, CaCl2, NaCl, and KCl, and for four different mutants of the human lysozyme: wild type, Q86D, A92D, and Q86D/A92D, which have been studied experimentally. The discussion of this article focuses on the cleft that consists of amino acid residues from Q86 to A92. For the wild type of protein in the aqueous solutions of all the electrolytes studied, there are no distributions observed for the ions inside the cleft. The Q86D mutant shows essentially the same behavior with that of the wild type. The A92D mutant shows strong binding ability to Na+ in the recognition site, which is in accord with the experimental results. There are two isomers of the Q86D/A92D mutant, e.g., apo-Q86D/A92D and holo-Q86D/A92D. Although both isomers exhibit the binding ability to the Na+ and Ca2+ ions, the holo isomer shows much greater affinity compared with the apo isomer. Regarding the selective ion binding of the holo-Q86D/A92D mutant, it shows greater affinity to Ca2+ than to Na+, which is also consistent with the experimental observation.