Two‐dimensional replica‐exchange method for predicting protein–ligand binding structures

We have developed a two‐dimensional replica‐exchange method for the prediction of protein–ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein‐ligand complex systems (MDM2 and HSP 90‐alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica‐exchange umbrella sampling method that we have previously developed. © 2013 Wiley Periodicals, Inc.     

Absolute Protein Binding Free Energy Simulations for Ligands with Multiple Poses,a Thermodynamic Path That Avoids Exhaustive Enumeration of the Poses

We propose a free energy calculation method for receptor–ligand binding, which have multiple binding poses that avoids exhaustive enumeration of the poses. For systems with multiple binding poses, the standard procedure is to enumerate orientations of the binding poses, restrain the ligand to each orientation, and then, calculate the binding free energies for each binding pose. In this study, we modify a part of the thermodynamic cycle in order to sample a broader conformational space of the ligand in the binding site. This modification leads to more accurate free energy calculation without performing separate free energy simulations for each binding pose. We applied our modification to simple model host–guest systems as a test, which have only two binding poses, by using a single decoupling method (SDM) in implicit solvent. The results showed that the binding free energies obtained from our method without knowing the two binding poses were in good agreement with the benchmark results obtained by explicit enumeration of the binding poses. Our method is applicable to other alchemical binding free energy calculation methods such as the double decoupling method (DDM) in explicit solvent. We performed a calculation for a protein–ligand system with explicit solvent using our modified thermodynamic path. The results of the free energy simulation along our modified path were in good agreement with the results of conventional DDM, which requires a separate binding free energy calculation for each of the binding poses of the example of phenol binding to T4 lysozyme in explicit solvent. © 2019 Wiley Periodicals, Inc.     

Efficient calculation of relative binding free energies by umbrella sampling perturbation

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide‐protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1–2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation. © 2014 Wiley Periodicals, Inc.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

FDS: flexible ligand and receptor docking with a continuum solvent model and soft-core energy function

The docking of flexible small molecule ligands to large flexible protein targets is addressed in this article using a two-stage simulation-based method. The methodology presented is a hybrid approach where the first component is a dock of the ligand to the protein binding site, based on deriving sets of simultaneously satisfied intermolecular hydrogen bonds using graph theory and a recursive distance geometry algorithm. The output structures are reduced in number by cluster analysis based on distance similarities. These structures are submitted to a modified Monte Carlo algorithm using the AMBER-AA molecular mechanics force field with the Generalized Born/Surface Area (GB/SA) continuum model. This solvent model is not only less expensive than an explicit representation, but also yields increased sampling. Sampling is also increased using a rotamer library to direct some of the protein side-chain movements along with large dihedral moves. Finally, a softening function for the nonbonded force field terms is used, enabling the potential energy function to be slowly turned on throughout the course of the simulation. The docking procedure is optimized, and the results are presented for a single complex of the arabinose binding protein. It was found that for a rigid receptor model, the X-ray binding geometry was reproduced and uniquely identified based on the associated potential energy. However, when side-chain flexibility was included, although the X-ray structure was identified, it was one of three possible binding geometries that were energetically indistinguishable. These results suggest that on relaxing the constraint on receptor flexibility, the docking energy hypersurface changes from being funnel-like to rugged. A further 14 complexes were then examined using the optimized protocol. For each complex the docking methodology was tested for a fully flexible ligand, both with and without protein side-chain flexibility. For the rigid protein docking, 13 out of the 15 test cases were able to find the experimental binding mode; this number was reduced to 11 for the flexible protein docking. However, of these 11, in the majority of cases the experimental binding mode was not uniquely identified, but was present in a cluster of low energy structures that were energetically indistinguishable. These results not only support the presence of a rugged docking energy hypersurface, but also suggest that it may be necessary to consider the possibility of more than one binding conformation during ligand optimization.     

Assessment of programs for ligand binding affinity prediction

The prediction of the binding free energy between a ligand and a protein is an important component in the virtual screening and lead optimization of ligands for drug discovery. To determine the quality of current binding free energy estimation programs, we examined FlexX, X-Score, AutoDock, and BLEEP for their performance in binding free energy prediction in various situations including cocrystallized complex structures, cross docking of ligands to their non-cocrystallized receptors, docking of thermally unfolded receptor decoys to their ligands, and complex structures with "randomized" ligand decoys. In no case was there a satisfactory correlation between the experimental and estimated binding free energies over all the datasets tested. Meanwhile, a strong correlation between ligand molecular weight-binding affinity correlation and experimental predicted binding affinity correlation was found. Sometimes the programs also correctly ranked ligands' binding affinities even though native interactions between the ligands and their receptors were essentially lost because of receptor deformation or ligand randomization, and the programs could not decisively discriminate randomized ligand decoys from their native ligands; this suggested that the tested programs miss important components for the accurate capture of specific ligand binding interactions.     

BLEEP—potential of mean force describing protein–ligand interactions: I. Generating potential

We have developed BLEEP (biomolecular ligand energy evaluation protocol), an atomic level potential of mean force (PMF) describing protein–ligand interactions. The pair potentials for BLEEP have been derived from high-resolution X-ray structures of protein–ligand complexes in the Brookhaven Protein Data Bank (PDB), with a careful treatment of homology. The use of a broad variety of protein–ligand structures in the derivation phase gives BLEEP more general applicability than previous potentials, which have been based on limited classes of complexes, and thus represents a significant step forward. We calculate the distance distributions in protein–ligand interactions for all 820 possible pairs that can be chosen from our set of 40 different atom types, including polar hydrogen. We then use a reverse Boltzmann methodology to convert these into energy-like pair potential functions. Two versions of BLEEP are calculated, one including and one excluding interactions between protein and water. The pair potentials are found to have the expected forms; polar and hydrogen bonding interactions show minima at short range, around 3.0 Å, whereas a typical hydrophobic interaction is repulsive at this distance, with values above 4.0 Å being preferred. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1165–1176, 1999     

Structural mining: self-consistent design on flexible protein-peptide docking and transferable binding affinity potential

A flexible protein-peptide docking method has been designed to consider not only ligand flexibility but also the flexibility of the protein. The method is based on a Monte Carlo annealing process. Simulations with a distance root-mean-square (dRMS) virtual energy function revealed that the flexibility of protein side chains was as important as ligand flexibility for successful protein-peptide docking. On the basis of mean field theory, a transferable potential was designed to evaluate distance-dependent protein-ligand interactions and atomic solvation energies. The potential parameters were developed using a self-consistent process based on only 10 known complex structures. The effectiveness of each intermediate potential was judged on the basis of a Z score, approximating the gap between the energy of the native complex and the average energy of a decoy set. The Z score was determined using experimentally determined native structures and decoys generated by docking with the intermediate potentials. Using 6600 generated decoys and the Z score optimization criterion proposed in this work, the developed potential yielded an acceptable correlation of R(2) = 0.77, with binding free energies determined for known MHC I complexes (Class I Major Histocompatibility protein HLA-A(*)0201) which were not present in the training set. Test docking on 25 complexes further revealed a significant correlation between energy and dRMS, important for identifying native-like conformations. The near-native structures always belonged to one of the conformational classes with lower predicted binding energy. The lowest energy docked conformations are generally associated with near-native conformations, less than 3.0 Angstrom dRMS (and in many cases less than 1.0 Angstrom) from the experimentally determined structures.     

A discriminative Ramachandran potential of mean force aimed at minimizing secondary structure bias

We introduce PMF*, a novel potential of mean force (PMF) for the Ramachandran ?/Ψ dihedral plot of the 20 standard amino acids and assess its relevance to the conformation of polypeptides by scoring structures in the protein data bank and decoy datasets. The new energy function is a linear combination of the conventional, unreferenced PMF and the ΔPMF relative to the free energy of all amino acids in the parameterization set of structures, effectively removing their respective biases toward α‐helix and β‐strand. It is shown that low‐resolution crystal structures, NMR structures, and theoretical models have on average significantly higher energies than high‐resolution crystal structures; also PMF* is more discriminative for structure quality than the individual PMF and ΔPMF energy functions. PMF* may be well suited for use as a restraint energy term in the refinement of experimental structures and theoretical models. © 2012 Wiley Periodicals, Inc.     

Free energy calculation using molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. II. Thermodynamic integration along a spatial reaction coordinate

We propose the thermodynamic integration along a spatial reaction coordinate using the molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. This method provides a free energy calculation in solution along the reaction coordinate defined by the Cartesian coordinates of the solute atoms. The proposed method is based on the blue moon algorithm which can, in principle, handle any reaction coordinate as far as it is defined by the solute atom positions. In this article, we apply the present method to the complex formation process of the crown ether 18-Crown-6 (18C6) with the potassium ion in an aqueous solution. The separation between the geometric centers of these two molecules is taken to be the reaction coordinate for this system. The potential of mean force (PMF) becomes the maximum at the separation between the molecular centers being ～4 A?, which can be identified as the free energy barrier in the process of the molecular recognition. In a separation further than the free energy barrier, the PMF is slightly reduced to exhibit a plateau. In the region closer than the free energy barrier, approach of the potassium ion to the center of 18C6 also decreases the PMF. When the potassium ion is accommodated at the center of 18C6, the free energy is lower by -5.7 ± 0.7 kcal/mol than that at the above mentioned plateau or converged state. By comparing the results with those from the free energy calculation along the coupling parameters obtained in our previous paper [T. Miyata, Y. Ikuta, and F. Hirata, J. Chem. Phys. 133, 044114 (2010)], it is found that the effective interaction in water between 18C6 and the potassium ion vanishes beyond the molecular-center-separation of 10 A?. Furthermore, the conformation of 18C6 is found to be significantly changed depending upon the 18C6-K(+) distance. A proper conformational sampling and an accurate solvent treatment are crucial for realizing the accurate PMF, and we believe that the proposed method is useful to evaluate the PMF in a solution. A discussion upon the PMF in terms of the three-dimensional distribution function for the solvent is also presented.     

Molecular dynamics study of the inclusion of cholesterol into cyclodextrins

The interaction of three cyclodextrins (CDs), viz. beta-CD, heptakis (2,6-di-O-methyl)-beta-CD (DM-beta-CD), and 2-hydroxypropyl-beta-CD (HP-beta-CD), with cholesterol was investigated using molecular dynamics (MD) simulations. The free energy along the reaction pathway delineating the inclusion of cholesterol into each CD was computed using the adaptive biasing force method. The association constant and the corresponding association free energy were derived by integrating the potential of mean force (PMF) over a representative ordering parameter. The results show that the free energy profiles possess two local minima corresponding to roughly equally probable binding modes. Among the three CDs, DM-beta-CD exhibits the highest propensity to associate with cholesterol. Ranking for binding cholesterol, viz. DM-beta-CD > HP-beta-CD > beta-CD, agrees nicely with experiment. Partitioning of the PMF into free energy components illuminates that entering of cholesterol into the CD cavity is driven mainly by electrostatic interactions, whereas deeper inclusion results from van der Waals forces and solvation effects. Additional MD simulations were performed to investigate the structural stability of the host-guest complexes near the free energy minima. The present results demonstrate that association of cholesterol and CDs follows two possible binding modes. Although the latter are thermodynamically favorable for all CDs, one of the two inclusion complexes appears to be preferred kinetically in the case of DM-beta-CD.     

A water-swap reaction coordinate for the calculation of absolute protein-ligand binding free energies

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.     

Protein:Ligand binding free energies: A stringent test for computational protein design

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, we consider tyrosyl‐tRNA synthetase (TyrRS), which has been extensively redesigned experimentally. We consider its specificity for its substrate l ‐tyrosine (l ‐Tyr), compared to the analogs d ‐Tyr, p‐acetyl‐, and p‐azido‐phenylalanine (ac‐Phe, az‐Phe). We simulate l ‐ and d ‐Tyr binding to TyrRS and six mutants, and compare the structures and binding free energies to a more rigorous “MD/GBSA” procedure: molecular dynamics with explicit solvent for structures and a Generalized Born + Surface Area model for binding free energies. Next, we consider l ‐Tyr, ac‐ and az‐Phe binding to six other TyrRS variants. The titration results are sensitive to the precise rotamer definition, which involves a short energy minimization for each sidechain pair to help relax bad contacts induced by the discrete rotamer set. However, when designed mutant structures are rescored with a standard GBSA energy model, results agree well with the more rigorous MD/GBSA. As a third test, we redesign three amino acid positions in the substrate coordination sphere, with either l ‐Tyr or d ‐Tyr as the ligand. For two, we obtain good agreement with experiment, recovering the wildtype residue when l ‐Tyr is the ligand and a d ‐Tyr specific mutant when d ‐Tyr is the ligand. For the third, we recover His with either ligand, instead of wildtype Gln. © 2015 Wiley Periodicals, Inc.     

Delineating Protein–Protein Curvilinear Dissociation Pathways and Energetics with Naïve Multiple-Walker Umbrella Sampling Simulations

The protein–protein interaction energetics can be obtained by calculating the potential of mean force (PMF) from umbrella sampling (US) simulations, in which samplings are often enhanced along a predefined vector as the reaction coordinate. However, any slight change in the vector may significantly vary the calculated PMF, and therefore the energetics using a random choice of vector may mislead. A non-predefined curve path-based sampling enhancement approach is a natural alternative, but was relatively less explored for protein–protein systems. In this work, dissociation of the barnase–barstar complex is simulated by implementing non-predefined curvilinear pathways in US simulations. A simple variational principle is applied to determine the lower bound PMF, which could be used to derive the standard free energy of binding. Two major dissociation pathways, which include interactions with the RNA-binding loop and the Val 36 to Gly 40 loop, are observed. Further, the proposed approach was used to discriminate the decoys from protein–protein docking studies. © 2019 Wiley Periodicals, Inc.     

An analysis of the interactions between the Sem-5 SH3 domain and its ligands using molecular dynamics, free energy calculations, and sequence analysis.

The Src-homology-3 (SH3) domain of the Caenorhabditis elegans protein Sem-5 binds proline-rich sequences. It is reported that the SH3 domains broadly accept amide N-substituted residues instead of only recognizing prolines on the basis of side chain shape or rigidity. We have studied the interactions between Sem-5 and its ligands using molecular dynamics (MD), free energy calculations, and sequence analysis. Relative binding free energies, estimated by a method called MM/PBSA, between different substitutions at sites -1, 0, and +2 of the peptide are consistent with the experimental data. A new method to calculate atomic partial charges, AM1-BCC method, is also used in the binding free energy calculations for different N-substitutions at site -1. The results are very similar to those obtained from widely used RESP charges in the AMBER force field. AM1-BCC charges can be calculated more rapidly for any organic molecule than can the RESP charges. Therefore, their use can enable a broader and more efficient application of the MM/PBSA method in drug design. Examination of each component of the free energy leads to the construction of van der Waals interaction energy profiles for each ligand as well as for wild-type and mutant Sem-5 proteins. The profiles and free energy calculations indicate that the van der Waals interactions between the ligands and the receptor determine whether an N- or a Calpha-substituted residue is favored at each site. A VC value (defined as a product of the conservation percentage of each residue and its van der Waals interaction energy with the ligand) is used to identify several residues on the receptor that are critical for specificity and binding affinity. This VC value may have a potential use in identifying crucial residues for any ligand-protein or protein-protein system. Mutations at two of those crucial residues, N190 and N206, are examined. One mutation, N190I, is predicted to reduce the selectivity of the N-substituted residue at site -1 of the ligand and is shown to bind similarly with N- and Calpha-substituted residues at that site.     

Hamiltonian and distance replica exchange method studies of Met-enkephalin

The conformational states of the zwitterionic form of the pentapeptide Met-enkephalin were explored with the use of explicit solvent molecular dynamics (MD). The N and C termini are ionized, as appropriate to polar solvent conditions, and consequently, there is a competition between open forms driven by polar solvation of the ammonium and carboxylate groups and closed forms driven by their salt-bridge formation. Normal MD started from an open state does not sample closed conformations. Sampling was enhanced with a distance replica exchange method (DREM) and with a Hamiltonian replica exchange method (HREM). The potential of mean force (PMF) along an end-to-end distance reaction coordinate was obtained with the DREM. The PMF shows a stable salt-bridge state and the presence of a large region of open states, as hypothesized for conformationally promiscuous small opiate peptides. The HREM systems differ by scaling the peptide-peptide and peptide-solvent electrostatic and Lennard-Jones potentials, with the goal of improving the sampling efficiency with a limited number of systems. A small number of systems were found to be sufficient to sample closed and open states. A principal component analysis (PCA) shows that the HREM-generated fluctuations are dominated by the first two principal modes. The first corresponds to the end-to-end reaction coordinate found in the DREM, and the first mode PMF is similar to the DREM PMF. The second mode describes the presence of two conformations, both of which correspond to the salt-bridge state distance. The conformers differ in the values of neighboring psi and phi dihedral angles, since such psi/phi compensation can still produce the same end-to-end distance. The two-dimensional PMF constructed from the first two PCA modes captures most of the significant backbone conformational space of Met-enkephalin.     

苯磺酰胺从碳酸酐酶II中脱离过程的分子动力学模拟

综合运用分子动力学模拟和自由能计算方法研究了苯磺酰胺分子从碳酸酐酶II (CA II)的活性位点脱离过程中底物与酶之间的动态相互作用. 脱离过程的平均力势(PMF)显示, 底物脱离时存在一个特殊的结合状态. 其中, 静电相互作用占据了主导地位. 轨迹分析显示, 除了金属离子的配位作用之外, 底物脱离路径上的关键残基Leu198、Thr199和Thr200通过与底物磺胺基的氢键作用阻碍了底物从酶中的脱离. 当前的研究对于深入认识磺胺类药物与CA II的详细结合过程和相关的药物改良与设计具有重要的指导意义.     

Identifying ligand binding sites and poses using GPU-accelerated Hamiltonian replica exchange molecular dynamics

We present a method to identify small molecule ligand binding sites and poses within a given protein crystal structure using GPU-accelerated Hamiltonian replica exchange molecular dynamics simulations. The Hamiltonians used vary from the physical end state of protein interacting with the ligand to an unphysical end state where the ligand does not interact with the protein. As replicas explore the space of Hamiltonians interpolating between these states, the ligand can rapidly escape local minima and explore potential binding sites. Geometric restraints keep the ligands from leaving the vicinity of the protein and an alchemical pathway designed to increase phase space overlap between intermediates ensures good mixing. Because of the rigorous statistical mechanical nature of the Hamiltonian exchange framework, we can also extract binding free energy estimates for all putative binding sites. We present results of this methodology applied to the T4 lysozyme L99A model system for three known ligands and one non-binder as a control, using an implicit solvent. We find that our methodology identifies known crystallographic binding sites consistently and accurately for the small number of ligands considered here and gives free energies consistent with experiment. We are also able to analyze the contribution of individual binding sites to the overall binding affinity. Our methodology points to near term potential applications in early-stage structure-guided drug discovery.     

Inhibitor docking screened by the modified SAFE_p scoring function: application to cyclic urea HIV-1 PR inhibitors

Our laboratory has in the past developed a method for the prediction of ligand binding free energies to proteins, referred to as SAFE_p (Solvent free energy predictor). Previously, we have applied this protocol for the prediction of the binding free energy of peptidic and cyclic urea HIV-1 PR inhibitors, whose X-ray structures bound to enzyme are known. In this work, we present the first account of a docking simulation, where the ligand conformations were screened and inhibitor ranking was predicted on the basis of a modified SAFE_p approach, for a set of cyclic urea-HIV-1 PR complexes whose structures are not known. We show that the optimal dielectric constant for docking is rather high, in line with the values needed to reproduce some protein residue properties, like pKa's. Our protocol is able to reproduce most of the observed binding ranking, even in the case that the components of the equation are not fitted to experimental data. Partition of the binding free energy into pocket and residue contributions sheds light into the importance of the inhibitor's fragments and on the prediction of "hot spots" for resistance mutations.     

WATsite: Hydration site prediction program with PyMOL interface

Water molecules that mediate protein–ligand interactions or are released from the binding site on ligand binding can contribute both enthalpically and entropically to the free energy of ligand binding. To elucidate the thermodynamic profile of individual water molecules and their potential contribution to ligand binding, a hydration site analysis program WATsite was developed together with an easy‐to‐use graphical user interface based on PyMOL. WATsite identifies hydration sites from a molecular dynamics simulation trajectory with explicit water molecules. The free energy profile of each hydration site is estimated by computing the enthalpy and entropy of the water molecule occupying a hydration site throughout the simulation. The results of the hydration site analysis can be displayed in PyMOL. A key feature of WATsite is that it is able to estimate the protein desolvation free energy for any user specified ligand. The WATsite program and its PyMOL plugin are available free of charge from http://people.pnhs.purdue.edu/~mlill/software . © 2014 Wiley Periodicals, Inc.