Methods using both machine learning (ML) and statistics have been successfully applied to large data sets obtained from MD simulations of biomolecular systems, with these approaches typically involving the generation of large features sets through determining all of the Cα-Cα (or Cβ-Cβ) distances or backbone dihedrals on a per frame basis. One such application of these features sets is to use them to describe slow timescale processes occurring in biological systems, enabling the development of collective variables (i.e., reaction coordinates) that can describe this process in a coarse-grained manner. With these collective variables, complex processes can be well described with 1 or 2 dimensions and/or applied to methods such as metadynamics, in order to enhance the sampling of molecular dynamics simulations. Further, prior work has demonstrated how these large data sets can be exploited to identify (conformationally) important residues or groups of residues distributed throughout an entire biomolecular system. Here, we describe a Python package called "KIF" (short for Key Interactions Finder), which we have developed to provide a largely automated workflow that allows researchers to identify the key non-covalent interactions and residues that describe a change in any user defined conformational state. As we will demonstrate in this manuscript, these changes can be anything from relatively large-scale loop conformational changes all the way through to much more subtle changes in a single side chain conformation. To do this, we built upon the popular MD simulation analysis tool PyContact, which is an MD simulation engine-agnostic method to identify and score all the non-covalent interactions present in an MD simulation. By combining the per frame interaction data with a user defined target variable, one can apply KIF to identify the key non-covalent interactions distributed throughout the protein that are the most sensitive to the target variable's value. Further, KIF allows for the conformational change of interest to be described as a classification problem (i.e., State A vs State B) or as a regression problem (i.e., a continuous target variable that describes the conformational change), and the end user has complete flexibility in terms of deciding how to generate this target variable to describe their conformational change of interest. The software outputs a relative score for each interaction based on its relationship to the target variable, and alongside this, per residue scores can also be determined (by summing and normalizing the per interactions scores obtained for each residue). The per feature and per residue scores generated by our program can be outputted in tabular form or as PyMOL compatible scripts in order to visualize the results on a 3D structure.