The Role of Force Fields and Water Models in Protein Folding and Unfolding Dynamics

Protein folding is a fascinating, not fully understood phenomenon in biology. Molecular dynamics (MD) simulations are an invaluable tool to study conformational changes in atomistic detail, including folding and unfolding processes of proteins. However, the accuracy of the conformational ensembles derived from MD simulations inevitably relies on the quality of the underlying force field in combination with the respective water model. Here, we investigate protein folding, unfolding, and misfolding of fast-folding proteins by examining different force fields with their recommended water models, i.e., ff14SB with the TIP3P model and ff19SB with the OPC model. To this end, we generated long conventional MD simulations highlighting the perks and pitfalls of these setups. Using Markov state models, we defined kinetically independent conformational substates and emphasized their distinct characteristics, as well as their corresponding state probabilities. Surprisingly, we found substantial differences in thermodynamics and kinetics of protein folding, depending on the combination of the protein force field and water model, originating primarily from the different water models. These results emphasize the importance of carefully choosing the force field and the respective water model as they determine the accuracy of the observed dynamics of folding events. Thus, the findings support the hypothesis that the water model is at least equally important as the force field and hence needs to be considered in future studies investigating protein dynamics and folding in all areas of biophysics.

Addressing Challenges of Macrocyclic Conformational Sampling in Polar and Apolar Solvents: Lessons for Chameleonicity

We evaluated a workflow to reliably sample the conformational space of a set of 47 peptidic macrocycles. Starting from SMILES strings, we use accelerated molecular dynamics simulations to overcome high energy barriers, in particular, the cis-trans isomerization of peptide bonds. We find that our approach performs very well in polar solvents like water and dimethyl sulfoxide. Interestingly, the protonation state of a secondary amine in the ring only slightly influences the conformational ensembles of our test systems. For several of the macrocycles, determining the conformational distribution in chloroform turns out to be considerably more challenging. Especially, the choice of partial charges crucially influences the ensembles in chloroform. We address these challenges by modifying initial structures and the choice of partial charges. Our results suggest that special care has to be taken to understand the configurational distribution in apolar solvents, which is a key step toward a reliable prediction of membrane permeation of macrocycles and their chameleonic properties.

PEP-Patch: Electrostatics in Protein-Protein Recognition, Specificity, and Antibody Developability

The electrostatic properties of proteins arise from the number and distribution of polar and charged residues. Electrostatic interactions in proteins play a critical role in numerous processes such as molecular recognition, protein solubility, viscosity, and antibody developability. Thus, characterizing and quantifying electrostatic properties of a protein are prerequisites for understanding these processes. Here, we present PEP-Patch, a tool to visualize and quantify the electrostatic potential on the protein surface in terms of surface patches, denoting separated areas of the surface with a common physical property. We highlight its applicability to elucidate protease substrate specificity and antibody-antigen recognition and predict heparin column retention times of antibodies as an indicator of pharmacokinetics.

Cross-linking disulfide bonds govern solution structures of diabodies

In the last years, antibodies have emerged as a promising new class of therapeutics, due to their combination of high specificity with long serum half-life and low risk of side-effects. Diabodies are a popular novel antibody format, consisting of two Fv domains connected with short linkers. Like IgG antibodies, they simultaneously bind two target proteins. However, they offer altered properties, given their smaller size and higher rigidity. In this study, we conducted the-to our knowledge-first molecular dynamics (MD) simulations of diabodies and find a surprisingly high conformational flexibility in the relative orientation of the two Fv domains. We observe rigidifying effects through the introduction of disulfide bonds in the Fv -Fv interface and characterize the effect of different disulfide bond locations on the conformation. Additionally, we compare VH -VL orientations and paratope dynamics between diabodies and an antigen binding fragment (Fab) of the same sequence. We find mostly consistent structures and dynamics, indicating similar antigen binding properties. The most significant differences can be found within the CDR-H2 loop dynamics. Of all CDR loops, the CDR-H2 is located closest to the artificial Fv -Fv interface. All examined diabodies show similar VH -VL orientations, Fv -Fv packing and CDR loop conformations. However, the variant with a P14C-K64C disulfide bond differs most from the Fab in our measures, including the CDR-H3 loop conformational ensemble. This suggests altered antigen binding properties and underlines the need for careful validation of the disulfide bond locations in diabodies.

Lessons for Oral Bioavailability: How Conformationally Flexible Cyclic Peptides Enter and Cross Lipid Membranes

Cyclic peptides extend the druggable target space due to their size, flexibility, and hydrogen-bonding capacity. However, these properties impact also their passive membrane permeability. As the "journey" through membranes cannot be monitored experimentally, little is known about the underlying process, which hinders rational design. Here, we use molecular simulations to uncover how cyclic peptides permeate a membrane. We show that side chains can act as "molecular anchors", establishing the first contact with the membrane and enabling insertion. Once inside, the peptides are positioned between headgroups and lipid tails─a unique polar/apolar interface. Only one of two distinct orientations at this interface allows for the formation of the permeable "closed" conformation. In the closed conformation, the peptide crosses to the lower leaflet via another "anchoring" and flipping mechanism. Our findings provide atomistic insights into the permeation process of flexible cyclic peptides and reveal design considerations for each step of the process.

Challenges in antibody structure prediction

Advances in structural biology and the exponential increase in the amount of high-quality experimental structural data available in the Protein Data Bank has motivated numerous studies to tackle the grand challenge of predicting protein structures. In 2020 AlphaFold2 revolutionized the field using a combination of artificial intelligence and the evolutionary information contained in multiple sequence alignments. Antibodies are one of the most important classes of biotherapeutic proteins. Accurate structure models are a prerequisite to advance biophysical property predictions and consequently antibody design. Specialized tools used to predict antibody structures based on different principles have profited from current advances in protein structure prediction based on artificial intelligence. Here, we emphasize the importance of reliable protein structure models and highlight the enormous advances in the field, but we also aim to increase awareness that protein structure models, and in particular antibody models, may suffer from structural inaccuracies, namely incorrect cis-amide bonds, wrong stereochemistry or clashes. We show that these inaccuracies affect biophysical property predictions such as surface hydrophobicity. Thus, we stress the importance of carefully reviewing protein structure models before investing further computing power and setting up experiments. To facilitate the assessment of model quality, we provide a tool "TopModel" to validate structure models.

Assessing developability early in the discovery process for novel biologics

Beyond potency, a good developability profile is a key attribute of a biological drug. Selecting and screening for such attributes early in the drug development process can save resources and avoid costly late-stage failures. Here, we review some of the most important developability properties that can be assessed early on for biologics. These include the influence of the source of the biologic, its biophysical and pharmacokinetic properties, and how well it can be expressed recombinantly. We furthermore present in silico, in vitro, and in vivo methods and techniques that can be exploited at different stages of the discovery process to identify molecules with liabilities and thereby facilitate the selection of the most optimal drug leads. Finally, we reflect on the most relevant developability parameters for injectable versus orally delivered biologics and provide an outlook toward what general trends are expected to rise in the development of biologics.

Bispecific antibodies-effects of point mutations on CH3-CH3 interface stability

A new format of therapeutic proteins is bispecific antibodies, in which two different heavy chains heterodimerize to obtain two different binding sites. Therefore, it is crucial to understand and optimize the third constant domain (CH3-CH3) interface to favor heterodimerization over homodimerization, and to preserve the physicochemical properties, as thermal stability. Here, we use molecular dynamics simulations to investigate the dissociation process of 19 CH3-CH3 crystal structures that differ from each other in few point mutations. We describe the dissociation of the dimeric interface as a two-steps mechanism. As confirmed by a Markov state model, apart from the bound and the dissociated state, we observe an additional intermediate state, which corresponds to an encounter complex. The analysis of the interdomain contacts reveals key residues that stabilize the interface. We expect that our results will improve the understanding of the CH3-CH3 interface interactions and thus advance the developability and design of new antibodies formats.

pH-dependent structural diversity of profilin allergens determines thermal stability

The family of profilin allergens is a common class of proteins found in plants, viruses and various eukaryotes including mammals. Profilins are characterized by an evolutionary conserved structural fold, which is responsible for their cross-reactive nature of Immunoglobulin E (IgE) antibodies. Despite their high overall structural similarity, they exhibit substantial differences in their biophysical properties, such as thermal and pH stability. To understand the origin of these functional differences of Amb a 8, Art v 4 and Bet v 2, we performed constant pH molecular dynamics simulation in combination with Gaussian accelerated MD simulations. Depending on the respective protonation at different pH levels, we find distinct differences in conformational flexibility, which are consistent with experimentally determined melting temperatures. These variations in flexibility are accompanied by ensemble shifts in the conformational landscape and quantified and localized by residue-wise B-factors and dihedral entropies. These findings strengthen the link between flexibility of profilin allergens and their thermal stability. Thus, our results clearly show the importance of considering protonation dependent conformational ensembles in solution to elucidate biophysical differences between these structurally similar allergens.

Comparison of hydrophobicity scales for predicting biophysical properties of antibodies

While antibody-based therapeutics have grown to be one of the major classes of novel medicines, some antibody development candidates face significant challenges regarding expression levels, solubility, as well as stability and aggregation, under physiological and storage conditions. A major determinant of those properties is surface hydrophobicity, which promotes unspecific interactions and has repeatedly proven problematic in the development of novel antibody-based drugs. Multiple computational methods have been devised for in-silico prediction of antibody hydrophobicity, often using hydrophobicity scales to assign values to each amino acid. Those approaches are usually validated by their ability to rank potential therapeutic antibodies in terms of their experimental hydrophobicity. However, there is significant diversity both in the hydrophobicity scales and in the experimental methods, and consequently in the performance of in-silico methods to predict experimental results. In this work, we investigate hydrophobicity of monoclonal antibodies using hydrophobicity scales. We implement several scoring schemes based on the solvent-accessibility and the assigned hydrophobicity values, and compare the different scores and scales based on their ability to predict retention times from hydrophobic interaction chromatography. We provide an overview of the strengths and weaknesses of several commonly employed hydrophobicity scales, thereby improving the understanding of hydrophobicity in antibody development. Furthermore, we test several datasets, both publicly available and proprietary, and find that the diversity of the dataset affects the performance of hydrophobicity scores. We expect that this work will provide valuable guidelines for the optimization of biophysical properties in future drug discovery campaigns.

Grid inhomogeneous solvation theory for cross-solvation in rigid solvents

Grid Inhomogeneous Solvation Theory (GIST) has proven useful to calculate localized thermodynamic properties of water around a solute. Numerous studies have leveraged this information to enhance structure-based binding predictions. We have recently extended GIST toward chloroform as a solvent to allow the prediction of passive membrane permeability. Here, we further generalize the GIST algorithm toward all solvents that can be modeled as rigid molecules. This restriction is inherent to the method and is already present in the inhomogeneous solvation theory. Here, we show that our approach can be applied to various solvent molecules by comparing the results of GIST simulations with thermodynamic integration (TI) calculations and experimental results. Additionally, we analyze and compare a matrix consisting of 100 entries of ten different solvent molecules solvated within each other. We find that the GIST results are highly correlated with TI calculations as well as experiments. For some solvents, we find Pearson correlations of up to 0.99 to the true entropy, while others are affected by the first-order approximation more strongly. The enthalpy-entropy splitting provided by GIST allows us to extend a recently published approach, which estimates higher order entropies by a linear scaling of the first-order entropy, to solvents other than water. Furthermore, we investigate the convergence of GIST in different solvents. We conclude that our extension to GIST reliably calculates localized thermodynamic properties for different solvents and thereby significantly extends the applicability of this widely used method.

Explicit solvation thermodynamics in ionic solution: extending grid inhomogeneous solvation theory to solvation free energy of salt-water mixtures

Hydration thermodynamics play a fundamental role in fields ranging from the pharmaceutical industry to environmental research. Numerous methods exist to predict solvation thermodynamics of compounds ranging from small molecules to large biomolecules. Arguably the most precise methods are those based on molecular dynamics (MD) simulations in explicit solvent. One theory that has seen increased use is inhomogeneous solvation theory (IST). However, while many applications require accurate description of salt-water mixtures, no implementation of IST is currently able to estimate solvation properties involving more than one solvent species. Here, we present an extension to grid inhomogeneous solvation theory (GIST) that can take salt contributions into account. At the example of carbazole in 1 M NaCl solution, we compute the solvation energy as well as first and second order entropies. While the effect of the first order ion entropy is small, both the water-water and water-ion entropies contribute strongly. We show that the water-ion entropies are efficiently approximated using the Kirkwood superposition approximation. However, this approach cannot be applied to the water-water entropy. Furthermore, we test the quantitative validity of our method by computing salting-out coefficients and comparing them to experimental data. We find a good correlation to experimental salting-out constants, while the absolute values are overpredicted due to the approximate second order entropy. Since ions are frequently used in MD, either to neutralize the system or as a part of the investigated process, our method greatly extends the applicability of GIST. The use-cases range from biopharmaceuticals, where many assays require high salt concentrations, to environmental research, where solubility in sea water is important to model the fate of organic substances.

Correcting cis-trans-transgressions in macromolecular structure models

Many macromolecular X-ray and cryo-EM structure models deposited in the PDB contain biologically relevant small molecule ligands with unsaturated fatty acid acyl chains, whose cis-trans stereochemistry is incorrect. The molecules are either not properly defined in their stereochemical restraint files, or the proper stereochemistry is neglected during model building. Often, the same molecules appear in deposited models in both isomeric configurations, one of which is almost always incorrect, and the use of the same moiety (HET) identifier and restraint files in model refinement is wrong. We present case studies of frequently occurring molecules and a compilation of identified cases of C-C=C-C cis-trans geometry in the deposited structure models. Full listings of cis/trans torsion angles are provided for models with commonly occurring molecules to assist identification and correction of cis-trans errors and prevent inadvertent use of incorrect models. Caveats for users, advice for modellers and suggestions for remediation efforts with a simple but effective restraint file modification are provided.

Conformational Shifts of Stacked Heteroaromatics: Vacuum vs. Water Studied by Machine Learning

Stacking interactions play a crucial role in drug design, as we can find aromatic cores or scaffolds in almost any available small molecule drug. To predict optimal binding geometries and enhance stacking interactions, usually high-level quantum mechanical calculations are performed. These calculations have two major drawbacks: they are very time consuming, and solvation can only be considered using implicit solvation. Therefore, most calculations are performed in vacuum. However, recent studies have revealed a direct correlation between the desolvation penalty, vacuum stacking interactions and binding affinity, making predictions even more difficult. To overcome the drawbacks of quantum mechanical calculations, in this study we use neural networks to perform fast geometry optimizations and molecular dynamics simulations of heteroaromatics stacked with toluene in vacuum and in explicit solvation. We show that the resulting energies in vacuum are in good agreement with high-level quantum mechanical calculations. Furthermore, we show that using explicit solvation substantially influences the favored orientations of heteroaromatic rings thereby emphasizing the necessity to include solvation properties starting from the earliest phases of drug design.

Surprisingly Fast Interface and Elbow Angle Dynamics of Antigen-Binding Fragments

Fab consist of a heavy and light chain and can be subdivided into a variable (V H and V L ) and a constant region (C H 1 and C L ). The variable region contains the complementarity-determining region (CDR), which is formed by six hypervariable loops, shaping the antigen binding site, the paratope. Apart from the CDR loops, both the elbow angle and the relative interdomain orientations of the V H -V L and the C H 1-C L domains influence the shape of the paratope. Thus, characterization of the interface and elbow angle dynamics is essential to antigen specificity. We studied nine antigen-binding fragments (Fab) to investigate the influence of affinity maturation, antibody humanization, and different light-chain types on the interface and elbow angle dynamics. While the CDR loops reveal conformational transitions in the micro-to-millisecond timescale, both the interface and elbow angle dynamics occur on the low nanosecond timescale. Upon affinity maturation, we observe a substantial rigidification of the V H and V L interdomain and elbow-angle flexibility, reflected in a narrower and more distinct distribution. Antibody humanization describes the process of grafting non-human CDR loops onto a representative human framework. As the antibody framework changes upon humanization, we investigated if both the interface and the elbow angle distributions are changed or shifted. The results clearly showed a substantial shift in the relative V H -V L distributions upon antibody humanization, indicating that different frameworks favor distinct interface orientations. Additionally, the interface and elbow angle dynamics of five antibody fragments with different light-chain types are included, because of their strong differences in elbow angles. For these five examples, we clearly see a high variability and flexibility in both interface and elbow angle dynamics, highlighting the fact that Fab interface orientations and elbow angles interconvert between each other in the low nanosecond timescale. Understanding how the relative interdomain orientations and the elbow angle influence antigen specificity, affinity, and stability has broad implications in the field of antibody modeling and engineering.

Conformational Ensembles of Antibodies Determine Their Hydrophobicity

A major challenge in the development of antibody biotherapeutics is their tendency to aggregate. One root cause for aggregation is exposure of hydrophobic surface regions to the solvent. Many current techniques predict the relative aggregation propensity of antibodies via precalculated scales for the hydrophobicity or aggregation propensity of single amino acids. However, those scales cannot describe the nonadditive effects of a residue’s surrounding on its hydrophobicity. Therefore, they are inherently limited in their ability to describe the impact of subtle differences in molecular structure on the overall hydrophobicity. Here, we introduce a physics-based approach to describe hydrophobicity in terms of the hydration free energy using grid inhomogeneous solvation theory (GIST). We apply this method to assess the effects of starting structures, conformational sampling, and protonation states on the hydrophobicity of antibodies. Our results reveal that high-quality starting structures, i.e., crystal structures, are crucial for the prediction of hydrophobicity and that conformational sampling can compensate errors introduced by the starting structure. On the other hand, sampling of protonation states only leads to good results when combined with high-quality structures, whereas it can even be detrimental otherwise. We conclude by pointing out that a single static homology model may not be adequate for predicting hydrophobicity.

Antibodies exhibit multiple paratope states influencing VH-VL domain orientations

In the last decades, antibodies have emerged as one of the most important and successful classes of biopharmaceuticals. The highest variability and diversity of an antibody is concentrated on six hypervariable loops, also known as complementarity determining regions (CDRs) shaping the antigen-binding site, the paratope. Whereas it was assumed that certain sequences can only adopt a limited set of backbone conformations, in this study we present a kinetic classification of several paratope states in solution. Using molecular dynamics simulations in combination with experimental structural information we capture the involved conformational transitions between different canonical clusters and additional dominant solution structures occurring in the micro-to-millisecond timescale. Furthermore, we observe a strong correlation of CDR loop movements. Another important aspect when characterizing different paratope states is the relative VH/VL orientation and the influence of the distinct CDR loop states on the VH/VL interface. Conformational rearrangements of the CDR loops do not only have an effect on the relative VH/VL orientations, but also influence in some cases the elbow-angle dynamics and shift the respective distributions. Thus, our results show that antibodies exist as several interconverting paratope states, each contributing to the antibody's properties.

Solvation Thermodynamics in Different Solvents: Water-Chloroform Partition Coefficients from Grid Inhomogeneous Solvation Theory

Reliable information on partition coefficients plays a key role in drug development, as solubility decisively affects bioavailability. In a physicochemical context, the partition coefficient of a solute between two different solvents can be described as a function of solvation free energies. Hence, substantial scientific efforts have been made toward accurate predictions of solvation free energies in various solvents. The grid inhomogeneous solvation theory (GIST) facilitates the calculation of solvation free energies. In this study, we introduce an extended version of the GIST algorithm, which enables the calculation for chloroform in addition to water. Furthermore, GIST allows localization of enthalpic and entropic contributions. We test our approach by calculating partition coefficients between water and chloroform for a set of eight small molecules. We report a Pearson correlation coefficient of 0.96 between experimentally determined and calculated partition coefficients. The capability to reliably predict partition coefficients between water and chloroform and the possibility to localize their contributions allow the optimization of a compound's partition coefficient. Therefore, we presume that this methodology will be of great benefit for the efficient development of pharmaceuticals.

Local and Global Rigidification Upon Antibody Affinity Maturation

During the affinity maturation process the immune system produces antibodies with higher specificity and activity through various rounds of somatic hypermutations in response to an antigen. Elucidating the affinity maturation process is fundamental in understanding immunity and in the development of biotherapeutics. Therefore, we analyzed 10 pairs of antibody fragments differing in their specificity and in distinct stages of affinity maturation using metadynamics in combination with molecular dynamics (MD) simulations. We investigated differences in flexibility of the CDR-H3 loop and global changes in plasticity upon affinity maturation. Among all antibody pairs we observed a substantial rigidification in flexibility and plasticity reflected in a substantial decrease of conformational diversity. To visualize and characterize these findings we used Markov-states models to reconstruct the kinetics of CDR-H3 loop dynamics and for the first time provide a method to define and localize surface plasticity upon affinity maturation.

Protein-Protein Binding as a Two-Step Mechanism: Preselection of Encounter Poses during the Binding of BPTI and Trypsin

Biomolecular recognition between proteins follows complex mechanisms, the understanding of which can substantially advance drug discovery efforts. Here, we track each step of the binding process in atomistic detail with molecular dynamics simulations using trypsin and its inhibitor bovine pancreatic trypsin inhibitor (BPTI) as a model system. We use umbrella sampling to cover a range of unbinding pathways. Starting from these simulations, we subsequently seed classical simulations at different stages of the process and combine them to a Markov state model. We clearly identify three kinetically separated states (an unbound state, an encounter state, and the final complex) and describe the mechanisms that dominate the binding process. From our model, we propose the following sequence of events. The initial formation of the encounter complex is driven by long-range interactions because opposite charges in trypsin and BPTI draw them together. The encounter complex features the prealigned binding partners with binding sites still partially surrounded by solvation shells. Further approaching leads to desolvation and increases the importance of van der Waals interactions. The native binding pose is adopted by maximizing short-range interactions. Thereby side-chain rearrangements ensure optimal shape complementarity. In particular, BPTI’s P1 residue adapts to the S1 pocket and prime site residues reorient to optimize interactions. After the paradigm of conformation selection, binding-competent conformations of BPTI and trypsin are already present in the apo ensembles and their probabilities increase during this proposed two-step association process. This detailed characterization of the molecular forces driving the binding process includes numerous aspects that have been discussed as central to the binding of trypsin and BPTI and protein complex formation in general. In this study, we combine all these aspects into one comprehensive model of protein recognition. We thereby contribute to enhance our general understanding of this fundamental mechanism, which is particularly critical as the development of biopharmaceuticals continuously gains significance.

Macrocycle Cell Permeability Measured by Solvation Free Energies in Polar and Apolar Environments

The relation of surface polarity and conformational preferences is decisive for cell permeability and thus bioavailability of macrocyclic drugs. Here, we employ grid inhomogeneous solvation theory (GIST) to calculate solvation free energies for a series of six macrocycles in water and chloroform as a measure of passive membrane permeability. We perform accelerated molecular dynamics simulations to capture a diverse structural ensemble in water and chloroform, allowing for a direct profiling of solvent-dependent conformational preferences. Subsequent GIST calculations facilitate a quantitative measure of solvent preference in the form of a transfer free energy, calculated from the ensemble-averaged solvation free energies in water and chloroform. Hence, the proposed method considers how the conformational diversity of macrocycles in polar and apolar solvents translates into transfer free energies. Following this strategy, we find a striking correlation of 0.92 between experimentally determined cell permeabilities and calculated transfer free energies. For the studied model systems, we find that the transfer free energy exceeds the purely water-based solvation free energies as a reliable estimate of cell permeability and that conformational sampling is imperative for a physically meaningful model. We thus recommend this purely physics-based approach as a computational tool to assess cell permeabilities of macrocyclic drug candidates.

Conformational selection of allergen-antibody complexes-surface plasticity of paratopes and epitopes

Antibodies have the ability to bind various types of antigens and to recognize different antibody-binding sites (epitopes) of the same antigen with different binding affinities. Due to the conserved structural framework of antibodies, their specificity to antigens is mainly determined by their antigen-binding site (paratope). Therefore, characterization of epitopes in combination with describing the involved conformational changes of the paratope upon binding is crucial in understanding and predicting antibody-antigen binding. Using molecular dynamics simulations complemented with strong experimental structural information, we investigated the underlying binding mechanism and the resulting local and global surface plasticity in the binding interfaces of distinct antibody-antigen complexes. In all studied allergen-antibody complexes, we clearly observe that experimentally suggested epitopes reveal less plasticity, while non-epitope regions show high surface plasticity. Surprisingly, the paratope shows higher conformational diversity reflected in substantially higher surface plasticity, compared to the epitope. This work allows a visualization and characterization of antibody-antigen interfaces and might have strong implications for antibody-antigen docking and in the area of epitope prediction.

Solvation Free Energy as a Measure of Hydrophobicity: Application to Serine Protease Binding Interfaces

Solvation and hydrophobicity play a key role in a variety of biological mechanisms. In substrate binding, but also in structure-based drug design, the thermodynamic properties of water molecules surrounding a given protein are of high interest. One of the main algorithms devised in recent years to quantify thermodynamic properties of water is the grid inhomogeneous solvation theory (GIST), which calculates these features on a grid surrounding the protein. Despite the inherent advantages of GIST, the computational demand is a major drawback, as calculations for larger systems can take days or even weeks. Here, we present a GPU accelerated version of the GIST algorithm, which facilitates efficient estimates of solvation free energy even of large biomolecular interfaces. Furthermore, we show that GIST can be used as a reliable tool to evaluate protein surface hydrophobicity. We apply the approach on a set of nine different proteases calculating localized solvation free energies on the surface of the binding interfaces as a measure of their hydrophobicity. We find a compelling agreement with the hydrophobicity of their substrates, i.e., peptides, binding into the binding cleft, and thus our approach provides a reliable description of hydrophobicity characteristics of these biological interfaces.