Preferential interaction coefficients of proteins in mixed solvents are bulk thermodynamic parameters that relate molecular characteristics of protein solvation with solvent effects on protein thermodynamics. Because of their bulk nature, they give no insight in the molecular level nature of protein solvation. In this study, we develop a methodology which provides insight into the molecular anatomy of preferential interaction coefficients by elucidating protein solvation in mixed solvents. Our methodology makes use of molecular simulations and reveals intricacies of solvent-protein interactions which are not accounted for by less detailed models for solvent effects on protein thermodynamics. This is demonstrated for lysozyme in 30 vol% aqueous glycerol. We find that solvent regions near protein O- and N-atoms that favor the formation of multiple hydrogen-bonds with glycerol positively contribute to the preferential interaction coefficient (15 +/- 4) due to the preferential solvation by glycerol molecules with long residence times (>2 ns). Yet, the overall value of the preferential interaction coefficient is negative as solvent regions near protein surface loci with similar affinities for glycerol and water have a stronger negative contribution (-22 +/- 4). On the basis of these results, we discuss the current scope and future prospects of our methodology to understand solvent effects on protein thermodynamics.