Recent Developments and Applications of the MMPBSA Method

The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach has been widely applied as an efficient and reliable free energy simulation method to model molecular recognition, such as for protein-ligand binding interactions. In this review, we focus on recent developments and applications of the MMPBSA method. The methodology review covers solvation terms, the entropy term, extensions to membrane proteins and high-speed screening, and new automation toolkits. Recent applications in various important biomedical and chemical fields are also reviewed. We conclude with a few future directions aimed at making MMPBSA a more robust and efficient method.


INTRODUCTION
It is widely accepted that high-level quantum mechanical (QM) methods provide the most detailed and accurate description of molecular structures, dynamics, and functions. However, for many biochemical systems that are often too complex, and/or biochemical processes that are too long, classical approaches are more commonly employed due to their efficiency and reasonable accuracy. To model biochemical systems classically, both long-range polar and short-range nonpolar interactions are important for accurate and transferrable models (Perutz, 1978;Davis and McCammon, 1990;Honig and Nicholls, 1995). Among the classical simulation methods, the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach (Srinivasan et al., 1998;Kollman et al., 2000;Gohlke and Case, 2004;Yang et al., 2011;Miller et al., 2012; has emerged as an efficient and reliable method to model molecular recognition, such as for protein-ligand binding interactions. The development of binding free energy calculation methods has been a central focus in molecular simulations. Indeed, theoretically rigorous, but computationally expensive, free energy perturbation, and thermodynamic integration methods were both proposed much earlier than the MMPBSA method (Zwanzig, 1954;Bennett, 1976;Straatsma and McCammon, 1991). For simple and small systems where they can be reliably applied, these "exact" methods have been shown to be more accurate than the MMPBSA method. However, their applications to typically large and complex biomolecular recognition problems are quite limited. This is due to low efficiency and slow convergence. There are several key approximations utilized in the MMPBSA method that allow it to be used as an efficient and reasonable approximation for free energy simulations. The PBSA model is used so that the solvation contribution to the free energy is approximated by using a continuum solvent model. In addition, the method separately approximates enthalpy and entropy contributions to the free energy, with varying degrees of success as discussed below. There are other approximate methods that are similar to MMPBSA, such as the Mining Minima (M2) method and the Linear Interaction Energy (LIE) method. These are all categorized as end-point methods that focus on the end states of processes to compute the free energy changes (Head et al., 1997;Luo et al., 1999Luo et al., , 2001Luo and Gilson, 2000;Aqvist and Marelius, 2001;Chen et al., 2004;Chang et al., 2007;Moghaddam et al., 2011;Mikulskis et al., 2012;Muddana and Gilson, 2012;Muddana et al., 2014).
In this review, we focus on recent developments and applications of the MMPBSA method that have been reported since 2014, which was roughly the time when the last major MMPBSA review was published (Genheden and Ryde, 2015). In the following, we first review the improvements made to the MMPBSA method over the last few years. This is followed by recent applications of the method in various fields, with a focus on biomedical applications.

Overview of MMPBSA
The MMPBSA method is most often applied to the calculation of binding free energies ( G bind ) of small molecule ligands bound to large biomolecule receptors, although large inter-biomolecular recognitions are also often reported as reviewed below. The binding free energy of the bound ligand-receptor complex in an aqueous solvent ( G bind,aq ) can be approximated as (Srinivasan et al., 1998;Kollman et al., 2000): G bind,solv = G polar + G non−polar , where E MM , G bind,solv , and − T S represent the gas-phase molecular mechanical energy change, the solvation free energy change, and the conformational entropy change upon binding, respectively. All of these changes are computed via ensemble averaging over a large set of sampled conformations. E MM includes three terms calculated using molecular mechanics (MM): the covalent energy change ( E covalent ), the electrostatic energy change ( E electrostatic ), and the van der Waals energy change ( E vdW ). E covalent consists of changes in the bond terms ( E bond ), the angle terms ( E angle ), and the torsion terms ( E torsion ). The solvation free energy change ( G bind,solv ) is usually separated into polar and non-polar contributions ( G polar and G non-polar ). The entropy term is the most difficult to compute, and it is often approximated with a normal mode method using a few selected snapshots.
To compute all the ensemble averages in Equations (1-4) for a binding affinity calculation, the MMPBSA method usually begins with a molecular dynamics (MD) simulation of the complex using the single-trajectory approach, or three separate MD simulations of the complex, receptor, and ligand, respectively, in the multi-trajectory approach. A snapshot of a structure is taken at various time points during the production portion of the MD simulations. These snapshots are then used to calculate average values and uncertainties of various quantities of interest. The MD simulations are almost always conducted in an explicit solvent model to obtain the most accurate snapshots possible before carrying out any calculations that make use of them. It is important to obtain many different conformations, or snapshots, over a suitable timeframe for use in later statistical analysis as it is often not trivial to observe converged averaging, even in the relatively easier single-trajectory approach .
The PBE is based on the more fundamental Poisson equation: where ε ( r) is a predefined dielectric distribution function for the solvated molecular system, φ ( r) is the potential distribution function, and ρ f ( r) is the fixed atomic charge density. To model the electrostatic interaction due to the additional presence of salt in the aqueous solution, the electrostatic potential φ( r) can be described by the PBE as: Here, λ( r) is a predefined ion-exclusion function with a value of 0 within the Stern layer and the molecular interior and a value of 1 outside the Stern layer. The salt-related term f (φ ( r)) is a function of the potential, φ ( r) , the valence, z i , of ion type i, and the bulk concentration, c i , at a given temperature T. When both the ionic strength and electric field are weak, the PBE can be linearized for easier numerical solutions: where κ 2 = 8π e 2 I ε v k B T . Here ε ν denotes the solvent dielectric constant, and I represents the ionic strength of the solution.
Over the past few years, a few new algorithm developments were reported for the numerical solution of the PBE (Xie, 2014;Fisicaro et al., 2016;Xie and Jiang, 2016). To deal with the singularity and nonlinearity of the PBE, Xie proposed a new decomposition and minimization scheme, together with a new proof on the existence and uniqueness of the PBE solution. A new PBE finite element solver was developed based on these solution decomposition and minimization techniques (Xie, 2014). Fisicaro et al. presented a preconditioned conjugate gradient technique to solve the generalized Poisson problem, and the linear regime of the PBE, in some 10 iterations. In combination with a self-consistent procedure, this technique was able to solve the non-linear Poisson-Boltzmann problem in a formulation including ionic steric effects (Fisicaro et al., 2016). Later Xie et al. incorporated nonlocal dielectric effects into the classic PBE for a protein in ionic solvent to derive a nonlocal modified Poisson-Boltzmann equation (NMPBE) and developed a finite element algorithm with a related package for solving the NMPBE (Xie and Jiang, 2016). Their results demonstrate the potential for the NMPBE to be a better predictor of electrostatic solvation and binding free energies compared to the standard PBE. It is worth noting that there has been a community wide push to explore alternative hardware for biomolecular simulations, such as the graphics processing units (GPU), which have a parallel architecture and are suited for high-performance computation with dense data parallelism (Colmenares et al., 2014a,b;Qi R. et al., 2017). A finite difference scheme with the successive over-relaxation method was implemented on the CUDA-based GPUs in the DelPhi package, which achieved a speedup of ∼10 times in the linear and non-linear cases (Colmenares et al., 2014b). More recently, Qi et al. implemented and analyzed commonly used linear PBE solvers on CUDA GPUs for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners (Qi R. et al., 2017). After extensive testing, the optimal GPU performance was observed using the Jacobi-preconditioned CG solver with a significant speedup that was up to 50 times faster than the standard CG solver on CPU. These progressive efforts on efficient numerical PBE solvers show great potential for accelerating MMPBSA computation.
Since the prior review (Genheden and Ryde, 2015), the numerical procedure and related factors for the widely used finite-difference method were also investigated for their impact on the MMPBSA method . This study showed that the impact of grid spacing on the quality of MMPBSA calculations is small in protein-ligand binding calculations; the agreement with experiment changed by a negligible amount when the grid spacing was changed from 0.50 to 0.25 Å. This indicated that the widely adopted default value of 0.50 Å used by the community was sufficient. The impact of different atomic radius sets and different molecular surface definitions was also analyzed, and weak influences were found on the agreement with experiment . This is probably due to the use of high protein dielectrics for the often-charged ligands and/or active sites as discussed below.
The effect of the solute dielectric constant was also investigated. A higher solute dielectric constant (using 2 or 4 instead of 1) was found to perform better in the virtual screening of ligands for tyrosine kinases (Sun et al., 2014a). Our own analysis of six groups of receptors reached a similar conclusion; the binding affinities using high dielectric constants (4 and 20) agreed better with experiment. The difference between calculations using dielectric constants of 4 and 20 was not very apparent except for the case of a highly charged binding pocket in one receptor . Aside from the study of higher solute dielectric constants, a residue-dependent dielectric model was also developed for use in an alanine scanning protocol with the MMPBSA method (Simoes et al., 2017). An attempt to modify the solute dielectric environment by incorporating structurally important, explicit water molecules in protein-ligand pockets for MMPBSA calculations was also reported, and it was found to improve the modeling of binding affinities for a series of JNK3 kinase inhibitors .
A hybrid QM/MM solute was also used in MMPBSA applications for predicting the binding affinities of FabI inhibitors (Su et al., 2015). The study suggested that the prediction results are sensitive to radii sets, GB methods, QM Hamiltonians, sampling protocols, and simulation length. The finding here appears to contradict our prior discussion. This is because the solute dielectric constant is set to 1 (vaccum) due to the explicit consideration of polarization in QM/MM studies. In the study of , a high solute dielectric constant was used to mimic polarization implicitly. In general, the use of high dielectric constants washes out the sensitivity to radii and surface definitions.

The Non-polar Solvation Term
The non-polar (non-electrostatic) solvation free energy contribution ( G non-polar ) to the solvation free energy change ( G bind,solv ) arises from the solute cavity formation within the solvent and van der Waals interactions between the solute and the solvent around the cavity (Weeks et al., 1971;Smith and Tanford, 1973;Chandler, 1977, 1980;Widom, 1982;Kang et al., 1987;Floris and Tomasi, 1989;Floris et al., 1991;Ashbaugh et al., 1999;Hummer, 1999;Lum et al., 1999;Gallicchio et al., 2000Gallicchio et al., , 2002Levy et al., 2003;Zacharias, 2003;Gallicchio and Levy, 2004;Su and Gallicchio, 2004;Dzubiella et al., 2006;Wagoner and Baker, 2006;Tan et al., 2007). Until recent times, the non-polar solvation free energy has been simply estimated to be proportional to the solvent accessible surface area (SASA) of the solute: This is referred to as the classical approach below. The surface tension γ and the correction term b are usually set to be constant for all solute molecules; for example, these are 0.00542 kcal/mol-Å 2 and 0.92 kcal/mol, respectively, in the AMBER package (Case et al., 2016). In more modern approaches (Tan et al., 2007), the cavity formation free energy and van der Waals (dispersion) free energy are modeled as separate terms because they scale differently vs. solute size. One way to correlate the cavity formation free energy is to use the volume (SAV) enclosed by the solvent accessible surface. A solvent accessible volume integration, or a solvent accessible surface integration, can be used to compute the dispersion term ( G dispersion ), so the total non-polar solvation free energy can be estimated as: These scaling factors apparently depend on the choices of atomic and solvent probe radii, for example, they are set to 0.0378 kcal/mol-Å 2 and −0.569 kcal/mol in the AMBER package (Case et al., 2016). It is interesting to note that the coefficient is much higher compared to the classical model since the dispersion term in solvation is always a negative term. The overall nonpolar solvation free energies are similar between the classical and modern models, at least for small molecules where all atoms are exposed. The performance of both the classical and modern nonpolar solvation models was analyzed, and it was found that the modern approach reduced the root-mean-square deviations of computed binding affinities (both relative and absolute) from experimental values while the correlations changed little from those computed using the classical approach .

The Entropy Term
The configurational entropy (S) in Equation (1) is often approximated by normal mode or quasi-harmonic analysis. Many proposed methods exist to calculate entropy (Kassem et al., 2015), but it is notoriously difficult to obtain a converged quantity. Thus, further approximations are often utilized; for example, usually only residues within a small sphere (radius of 8-12 Å) centered at the ligand and a limited number of snapshots (<100) are used for normal mode analysis. Furthermore, the single trajectory approach (i.e., that for the complex) is most often used in MMPBSA studies that do not consider any bindinginduced structural changes. In such an approach, configurational entropy computed by the normal mode analysis is often omitted completely in the ranking of relative binding affinities as its inclusion often does not improve the agreement with experiment . A few new ideas to approximate entropy were reported recently. For example, a new method termed BEERT (Binding Entropy Estimation of Rotation and Translation) was proposed to approximate configurational entropy changes in terms of the reduction in translational and rotational freedom of the ligand upon protein-ligand binding, starting from the flexible molecule approach (Ben-Shalornit et al., 2017). An interaction entropy (IE) method was also proposed to investigate the entropy change upon binding Duan et al., 2017). The interaction entropy (IE) contribution to binding free energy was defined as where E int pl is the fluctuation of the protein-ligand interaction energy for both electrostatic and van der Waals interactions.
The ensemble average of e β E int pl can be extracted from MD simulations, avoiding normal mode calculations. Both developers of BEERT and IE claim that these methods are highly efficient.

Extension to Membrane Proteins
The development of implicit membrane models based on existing continuum solvent approaches has advanced, making it possible to extend the MMPBSA method to biological membrane systems. The presence of an implicit membrane adds a complication to the numerical solution of the underlying PBE that is brought about by dielectric inhomogeneities that appear on the boundary surfaces of the computation grid. This issue can be alleviated by employing the periodic boundary condition, which is a common practice for electrostatic computations in MD simulations. The conjugate gradient and successive over-relaxation methods can be adjusted to take into account periodic calculations, but the convergence rate using either method is quite low. This limits their application to MMPBSA calculations which require that a large number of conformations be processed. To improve the convergence rate for use in biomolecular applications, the Incomplete Cholesky preconditioning method and the geometric multigrid method were both extended to incorporate periodicity (Botello- Smith and Luo, 2015). Applications to protein-ligand binding utilizing the newly developed membrane MMPBSA method were also reported (Greene et al., 2016;Xiao et al., 2017).

Extension to the High-Speed Screening of Ligands
The MMPBSA method has also been utilized as a rescoring method for docking applications. In protein-protein docking, it was found that MMPBSA rescoring was more capable at distinguishing correct complex structures from decoys than ZDOCK scoring in a test set of 46 protein-protein complexes for certain combinations of force fields, solvation models, and dielectric constants . A new binding affinity estimator, PBSA_E, that is based on MMPBSA inputs, was also optimized for protein-ligand binding. The optimization made use of a training set consisting of high-quality experimental data that was gathered from 145 complex structures. When the predicted binding affinities using PBSA_E were compared with the predicted binding affinities using other popular scoring functions such as GlideXP, GlideSP, and SYBYL_F, the PBSA_E method exhibited improved accuracy in terms of both achieving higher correlations with measured binding affinities and lower root-mean-square deviations . A study on MMPBSA in hierarchical virtual screening (HVS) was reported. This study demonstrated the predictability and validity of using the MMPBSA method for lead discovery as it identified novel inhibitors of the p38 MAP kinase by employing a physics-based Frontiers in Molecular Biosciences | www.frontiersin.org scoring function combined with a knowledge-based structural filter (Cao et al., 2014).
A study on the performance of MMPBSA in docking applications found that it depends on the choice of the solvent models among many factors that were analyzed. Specifically, the authors found that the choice of solvent models plays a minor role for one-protein-family/one-ligand cases which represent the unbiased protein-ligand complex sampling. However, for the total dataset with biased sampling, where some proteins and their homologs have an overabundant presence in the dataset, they found that numerical PB solvent methods do not perform as well as GB solvent methods. In addition, they showed that numerical PB methods were more sensitive to whether MD simulations were used for averaging. Such methods may be currently more suitable for individual protein binding free energy rankings where MD simulations can be easily conducted. This study also demonstrated that the numerical noise from screening applications that utilize only one or a few structures should be addressed in future developments of the dielectric model for numerical PB methods (Sun et al., 2014b).

Extension to Other High-Performance Analyses
MMPBSA was also adapted for high-performance mutational analysis. Single Amino Acid Mutation based change in Binding free Energy (SAAMBE) was developed to predict changes in binding free energy that are brought about by point mutations. SAAMBE utilized 3D structures of protein-protein complexes in a sequence-and structure-based approach. The method was centered around two components: a MMPBSAbased component, and a set of statistical terms obtained from the physico-chemical properties of protein complexes. A better agreement with experimentally determined binding free energy changes over a set of 1,300 mutations in 43 proteins indicated a significant improvement for predictions made using SAAMBE Petukh et al., 2015).
Additionally, MMPBSA was adapted into a novel structurebased multiscale approach to identify the key specificity determining residues (SDRs) of PDZ domains that appeared in explicit solvent MD simulations on PDZ-peptide complexes. SDRs were then used together with knowledge-based scoring functions in a proteome-wide search to locate their interaction partners (Tiwari and Mohanty, 2014).

New Toolkits
Over the past few years, several new toolkits were released to facilitate the use of MMPBSA calculations. A free energy workflow tool, FEW, was developed for AMBER to assist in the setup of molecular dynamics simulations in explicit membrane environments. FEW also assists in the setup and execution of effective binding free energy calculations for a single-trajectory implicit solvent/implicit membrane MMPBSA approach that involves multiple ligands binding to the same membrane protein (Homeyer and Gohlke, 2015). g_mmpbsa was developed for GROMACS, and it implemented the MMPBSA approach using subroutines written in-house or sourced from the GROMACS and APBS packages (Kumari et al., 2014).
GMXPBSA is another user-friendly suite of Bash/Perl scripts for streamlining MMPBSA calculations for GROMACS users (Paissoni et al., 2014). An easy-to-use pipeline tool named Calculation of Free Energy (CaFE) was published to facilitate both MMPBSA and LIE calculations. CaFe is capable of handling numerous static structure and molecular dynamics trajectory file formats generated by different molecular simulation packages, and it also supports various force field parameters (Liu and Hou, 2016).
A third common group of applications are for studies of targets in neural disorders. Targets such as dopamine D2 receptor , monoamine oxidase enzymes (Marsavelski and Vianello, 2017), and opioid G protein-coupled receptors  were analyzed as antipsychotic compounds. Here, Leonis et al. employed the recently reported crystal structure of the human κ-opioid receptor (κ-OR) to explain the binding mechanism with its antagonist JDTic and agonist SalA. Both JDTic and SalA show that they are capable of forming a favorable complex in the MMPBSA analysis, which was later confirmed by experiment . For Parkinson's disease, ligand-binding to the SUR1 receptor (Santos et al., 2016) and the adenosine receptor  was analyzed. Many targets were studied for Alzheimer's disease including β-secretase (Koukoulitsa et al., 2016), angiotensin-converting enzyme (Bhavaraju et al., 2016), acetylcholinesterase and butyrylcholinesterase (Kurt et al., 2017), and inhibitors directly targeting amyloid aggregation (Berhanu and Masunov, 2015). Ataxin-2 protein was studied for the treatment of spinocerebellar ataxia , superoxide dismutase for amyotrophic lateral sclerosis (Zhuang et al., 2016), and Niemann-Pick type Cl and C2 proteins (Poongavanam et al., 2016) that occur in rare neurodegenerative diseases.
There were a wide range of applications to targets outside the scope of pharmaceutical research. Although these studies are not related to drug discovery, accurate modeling of protein-ligand binding is also important in a wide variety of other contexts. For example, Starovoytov et al. calculated MMPBSA binding affinities of BPA-A, BPA-C, and BPA-D bound to human estrogenrelated receptor γ in their toxicology research. Their analysis showed that BPA-A was the strongest binder to the receptor (Starovoytov et al., 2014). Other similar studies include globin, for its catalytic mechanism in the hydrolysis of substituted phenyl hexanoates (Ercan et al., 2014), the role of conserved residues in substrate binding to Brassica rapa auxin amidohydrolase (Smolko et al., 2016), the effect of hydrophobic interactions in substrate binding to recombinant enzyme carboxylesterase (Shao et al., 2014), the substrate-enzyme interactions of endo-1,4-β-xylanase , azoreductase protein for the biodegradation of azo dyes Haghshenas et al., 2016), cytochrome P450 2A6 for nicotine addiction , Cel48F for producing bioethanol via fiber degradation (Qian M. D. et al., 2016), rubisco for biofuel production (Siqueira et al., 2016), streptavidin-biotin complex in biochemical sensing , sorotidine 5-monophosphate decarboxylase for its impressive rate enhancement (Jamshidi et al., 2014), tyrosyl-tRNA synthetases for genetic encoding of unnatural amino acids (Ren et al., 2015), T7 RNA polymerase for generating RNA labels (Borkotoky et al., 2016), AF9 in the YEATS family for the recognition of H3K9ac , cysteine protease 1 precursor from Zea for the hydrolysate of corn gluten meal , folate receptor alpha for producing milk with high folate concentration (Sahoo et al., 2014c), and both acid amido synthetase  and brassinosteroid (Lei et al., 2015) for plant growth.

Protein-Protein Binding Interactions
There are many applications of MMPBSA that involve the calculation of protein-protein or protein-peptide binding affinities. For anti-cancer applications, we saw studies on complexes of an ankyrin repeat with integrin-linked kinase binding to PINCH1 (Gautam et al., 2014(Gautam et al., , 2015, cyclindependent kinase 8 with cyclin C , Hsp90 with Cdc37-derived peptides , and p53 with MDMX . As an illustration, Wang et al. designed an eleven-residue peptide that was able to bind to Hsp90 with a predicted affinity of 6.9 mM, which was comparable to the experimental value of 3.0 mM . For anti-viral applications, research efforts were focused on analyzing interactions for ankyrin and domain III of the envelope protein of dengue virus II Dubey et al., 2017), a modular capsomere of a murine polyomavirus (MPV) VLP designed to protect against influenza , antibody recognition of immunoglobulin 2D1 for influenza virus H1N1 (Leong et al., 2015), differential structural dynamics and antigenicity of two HA-specific CTL epitopes binding to HLA-A * 0201 for influenza virus H5N1 (Sun and Liu, 2015), a complex of histone deacetylase 6 and ubiquitinspecific protease 5 for influenza virus A (Passos et al., 2016), and for a complex of an antibody and epitope variant of HIV-1 p24 capsid protein (Karim et al., 2015). For antibacterial applications, there were studies of salamander PGRP1 with its splice variant adPGRP1a for the innate immune system (Qi Z. T. et al., 2017), and an adsorption mechanism of human betadefensin-3 on bacterial membranes (Lee et al., 2016). Based on their MMPBSA analysis, Lee et al. found that the binding affinity of human beta-defensin-3 bound to a gram-positive membrane is over 3 times higher than when it is bound to a gram-negative membrane. Interestingly, this difference is mostly derived from electrostatic interactions, consistent with a net charge that is three times larger for gram-positive membranes compared to gram-negative membranes (Lee et al., 2016).
Studies of protein-protein interactions were also found for other diseases and disorders including serum paraoxonase 1 with high-density lipoprotein for antiatherosclerotic activity (Patra et al., 2014), vascular endothelial growth factor A with binding domains of anti-angiogenic agents for retinal neovascular degenerative diseases (Platania et al., 2015), angiotensinconverting enzyme with Angiotensin-II for hypertension , titin with T-cap/telethonin for dilated cardiomyopathy (Kumar D. T. et al., 2017), CC chemokine ligand 5 (CCL5) with human neutrophil peptide-1 (HNP1) for chronic inflammatory diseases (Wichapong et al., 2016), and ANKS6 with ANKS3 for polycystic kidney disease (Kan et al., 2016). In the study of CCL5 binding with human HNP1, Wichapong et al. reported a correlation coefficient of r = 0.66 between experiment and MMPBSA results for different species of CCL5. In addition, they confirmed that the entropy change upon binding was negligible in this study, so the entropy term could be ignored when a relative binding free energy was considered (Wichapong et al., 2016).
In addition, there were research efforts that focused on basic mechanisms such as protein stability and conformational dynamics (Bhavaraju and Hansmann, 2015;Getov et al., 2016). Many other studies were also reported, including an analysis of full length amylin oligomer aggregation (Berhanu and Masunov, 2014), effects of single point mutations on amyloid formation (Bhavaraju and Hansmann, 2015;Getov et al., 2016;Petukh et al., 2016), interactions between methylated histone H3 and effector domains of the PHD family in pursuit of a molecular mechanism of epigenetics (Grauffel et al., 2015), the binding mechanism of actin-depolymerizing factor 1 and G-actin (Du et al., 2016), the interaction between phosphotyrosine binding domains and peptides for neuronal development, immune responses, tissue homeostasis, and cell growth (Sain et al., 2016), the cognate transducer complex srII-htrII for the downstream signaling mechanism of sensory rhodopsin (Sahoo and Fujiwara, 2017), the binding of CBP to c-Myb for understanding the exact function of CBP and its interaction with c-Myb (Odoux et al., 2016), the catalytic stability of the tetrameric complex of cystathionine gamma-lyase (El-Sayed et al., 2015), the study of the mechanism of a molecular chaperone using acid-stress chaperone HdeA and its substrate protein (Zhou et al., 2017), and the binding of thiopeptide to a ribosomal subunit in order to understand the structure-activity relationship of thiostrepton derivatives (Wolf et al., 2014).
There are also several non-biomedical applications such as a study on the interactions between two adjacent gamma tubulins within the gamma-tubulin ring complex for growing yeast (Suri et al., 2014), between beta-sheet regions in corneous beta proteins of sauropsids to explain its stability and polymerization into filaments (Calvaresi et al., 2016), and the recognition of Avr protein by eukaryotic transcription factor xa5 of rice to understand the gene-for-gene mechanism that governs the direct interaction of R-Avr protein (Dehury et al., 2015). In the study by Suri et al. computational alanine scanning was employed to determine hotspot contributions in the interaction between gamma tubulins. Their analysis showed that most hotspot mutations reduce affinity by 1.16 kcal/mol, while for very crucial amino acids, such as Asp252 and Arg341, the affinity was decreased by more than 10 kcal/mol, which correlated well with experiment (Suri et al., 2014).

Complexes Involving Nucleic Acids
For DNA-protein interactions, there were mechanistic studies involving proteins such as HU, one of the major nucleoidassociated proteins for stabilizing DNA bending (Kim et al., 2014), highly conserved chromatin protein Cren7 for cellular processes such as transcription, replication, and repair , MutS, which recognizes mismatched DNA in DNA repair using ATP (Ishida and Matsumoto, 2016), copper nucleases for predicting electrostatic interactions with B-DNA (Liu C. M. et al., 2016), and metallopeptides binding to the Drew-Dickerson dodecamer (Galindo-Murillo and Cheatham, 2014). For RNA-protein interactions, the binding of toll-like receptor 3 and 22 was studied in relation to fish viral diseases , and the RNA-recognition motif of RNAbinding protein was analyzed for anti-cancer actitivies (Chang et al., 2016). Two studies on DNA-DNA interactions were also reported, including the calculation of binding affinities of short double stranded oligonucleotides (Yesudas et al., 2015), of three-quartet intramolecular human telomeric DNA G-quadruplexes (Islam et al., 2016), and of strand-strand interactions in a human telomeric tetrameric quadruplex (Chaubey et al., 2015). Yesudas et al. studied oligonucleotides (9-20 mers DNA) with the "3-trajectory" approach, without counter ions, by using GBSA, PBSA, and 3D-RISM-KH methods to calculate binding free energies. They showed that the 3D-RISM-KH method was in better agreement with the experimental data for larger oligonucleotides while the GBSA method performed better for smaller oligonucleotides (Yesudas et al., 2015).
There were also several studies involving small ligands binding to nucleic acids. Interactions of anticancer compounds with DNA, such as cisplatin and oxaliplatin (Jalili et al., 2016), distamycin with the DNA minor groove (Jalili and Maddah, 2017), and plant alkaloid chelerythrine with the human telomere sequence (Ghosh et al., 2015) were reported. Aminobenzimidazole binding to an internal ribosome entry site was studied for its anti-hepatitis C virus effect (Henriksen et al., 2014). The interactions of histone-derived antimicrobial peptides buforin II and DesHDAP1 with DNA were also investigated (Sim et al., 2017). Basic mechanistic studies were carried out, including cationic porphyrin-anthraquinone hybrids binding to DNA duplexes (Arba and Tjahjono, 2015) and G-quadruplexes (Arba et al., 2016), as well as ligands binding to riboswitches upon mutation . The interactions of reactive metabolites of anticancer compounds with DNA were also analyzed (Tumbi et al., 2014). In the binding of ligands to riboswitches, Hu et al. analyzed the relative binding free energies for a guanine riboswitch (GR) and a GUA complex relative to three complexes: 6GU (3.4 kcal/mol), 2BP (5.48 kcal/mol), and XAN (6.19 kcal/mol). These values were in good agreement with experimental observations of 3.21, 4.12, and 5.47 kcal/mol, respectively .

Guest-Host and Nano Systems
For guest-host systems, several studies were reported. Two octa acid hosts complexed with six guest molecules were analyzed by Bhakat and Soderhjelm for resolving the problem of trapped water in binding cavities. They performed a well-tempered funnel metadynamics (WT-FM) and MMPBSA analyses for the two octa-acid hosts, OAH (without methyl groups) and OAMe (with methyl groups), using both GAFF and OPLS force fields. Their analyses showed that the binding affinities of WT-FM are basically similar to experimental values while MMPBSA results have errors in the range of 5-10 kcal/mol due to its approximate nature (Bhakat and Soderhjelm, 2017). A binding interaction was analyzed to show that electrostatic interactions have the largest contribution to the stability of the cucurbituril-pseudorotaxane complex (Malhis et al., 2015). Complex stability was analyzed for multiple systems, such as the 1:1 and 1:2 inclusion complexes formed by nor-Secocucurbit[10]uril and 1-adamantanmethylammonium in water , cyclodextrin-Ibuprofen complexes , E-selectin-oligosaccharide complexes (Barra et al., 2017), and the naringenin-2,6-dimethyl βcyclodextrin inclusion complex (Sangpheak et al., 2014). Finally, enantiomeric discrimination of chiral organic salts by chiral aza-15-crown-5 ether with C1 symmetry was reported (Kocakaya et al., 2015).
For nano systems, a mechanism explaining how C-60 can block potassium ion channels was proposed (Calvaresi et al., 2015). The authors showed that a new binding site for C-60 exists in the channel cavity at the intracellular entrance of the selectivity filter. The escape barrier from the binding site is ∼21 kcal/mol as calculated via the umbrella sampling method, in good agreement with the MMPBSA result. MMPBSA was also used in the theoretical design of the cyclic lipopeptide nanotube as a molecular channel in the lipid bilayer Khavani et al., 2015Khavani et al., , 2017, in the study of the enzyme immobilization mechanism of alpha-chymotrypsin onto carbon nanotubes in organic media (Zhang L. Y. et al., 2015), and the mechanism of carbon nanotube activation of subtilisin Carlsberg in polar and non-polar organic media .

Monomer Stability
For single biomolecule stability research, we saw reports on the structures and energies of the alternate frame folding calbindin-D 9k protein (Tong et al., 2015), the effect of copper ions in the stability and structural change of human growth protein (Tazikeh-Lemeski, 2016), and the role of potassium in stabilizing the human telomeric intra-molecular G-quadruplex structure Wang and Liu, 2017). In their study of calbindin-D 9k , Tong et al. found two transition states and an intermediate state with a first rate-controlling barrier of 4.7 kcal/mol and a second barrier of 1.7 kcal/mol using MMPBSA, both of which are in good agreement with experiment (Tong et al., 2015).

CURRENT LIMITATIONS AND FUTURE DIRECTIONS
MMPBSA methods are widely applied to calculate binding affinities at a reasonable computational cost. These computational analyses have provided a large volume of valuable predictive results in a wide variety of studies. Even though MMPBSA is known to be less accurate than some of the more computationally expensive methods, like the free energy perturbation and thermodynamic integration methods, the qualitative agreement is often good enough to aid collaborative efforts involving both computational and experimental researchers. Developers are also actively working to improve MMPBSA methods for higher accuracy and efficiency by introducing better solute and solvent models, by porting the expensive energy computation (mostly involving the solvation terms) to faster GPU platforms, and by improving entropy estimations. Efforts are also needed to extend the MMPBSA method for various screening purposes that involve a large number of ligands and/or mutations to achieve a higher overall level of accuracy and efficiency. It is apparent that applications of the MMPBSA method have grown considerably in many different areas of biomolecular study. Most of these applications involve protein-ligand binding affinity calculations due to their utility in drug discovery research efforts. There are also many applications in the study of biomacromolecular complexes. It is also noteworthy that a few guest-host and nano systems utilized MMPBSA calculations, indicating a wider development space for this method in the future.