Molecular dynamics with explicit solvent is favored for its ability to more cor- rectly simulate aqueous biological processes, and has become routine thanks to in- creasingly powerful computational resources. However, analysis techniques including Markov state models (MSMs) ignore solvent atoms and focus solely on solute coordi- nates despite solvent being implicated in myriad biological phenomena. We present a unified framework called “Solvent-Shells Featurization” for including solvent degrees of freedom in analysis, and show that this method produces better models. We apply this method to simulations of dewetting in the two-domain protein BphC to generate a predictive MSM and identify functional water molecules. Furthermore, the proposed methodology could be easily extended for building MSMs of any systems with indis- tinguishable components.