SIMFIT: a microcomputer software-toolkit for modelistic studies in biochemistry

A software package suitable for personal computers and designedto handle simulation and fitting problems related to the studyof biomolecules under pre-steady and steady state conditionsis presented, and its overall architecture as well as the implementedalgorithms illustrated. The peculiar features of the packageare: (i) integrated capability of simulating dynamic modelsand fitting to them experimental data; (ii) handling of stiffproblems; (iii) free use of algebraic as well as differentialequations; (iv) objective comparison of models of differentcomplexity. The above features are discussed through a numberof examples taken from the direct experience of the authorsin enzyme kinetics and ligand binding. Received on May 23, 1989; accepted on August 21, 1989     

Binding free energy calculations of adenosine deaminase inhibitors

The interactions between four inhibitors and adenosine deaminase (ADA) were examined by calculating their binding free energies after molecular dynamics simulations. A bonded model was used to represent the electrostatic potentials of the zinc coordination site. The charge distribution of the model was derived by using a two-stage electrostatic potential fitting calculations. The calculated binding free energies are in good agreement with the experimental data and the ranking of binding affinities is well reproduced. Notably, our findings suggest that non-polar contributions play an important role for ADA-inhibitor interactions.     

Fluorescence intensity and lifetime distribution analysis: toward higher accuracy in fluorescence fluctuation spectroscopy

Fluorescence fluctuation methods such as fluorescence correlation spectroscopy and fluorescence intensity distribution analysis (FIDA) have proven to be versatile tools for studying molecular interactions with single molecule sensitivity. Another well-known fluorescence technique is the measurement of the fluorescence lifetime. Here, we introduce a method that combines the benefits of both FIDA and fluorescence lifetime analysis. It is based on fitting the two-dimensional histogram of the number of photons detected in counting time intervals of given width and the sum of excitation to detection delay times of these photons. Referred to as fluorescence intensity and lifetime distribution analysis (FILDA), the technique distinguishes fluorescence species on the basis of both their specific molecular brightness and the lifetime of the excited state and is also able to determine absolute fluorophore concentrations. The combined information yielded by FILDA results in significantly increased accuracy compared to that of FIDA or fluorescence lifetime analysis alone. In this paper, the theory of FILDA is elaborated and applied to both simulated and experimental data. The outstanding power of this technique in resolving different species is shown by quantifying the binding of calmodulin to a peptide ligand, thus indicating the potential for application of FILDA to similar problems in the life sciences.     

Cooperative association of plasminogen with fibrinogen

We have examined the association of both Glu- and Lys-plasminogen to fibrinogen by ultracentrifugal sedimentation equilibrium in neutral isotonic buffer in the presence of aprotinin. The fibrinogen and plasminogens, which were homogenous, did not exhibit any self-association. In each association study, five different molar ratios of plasminogen to fibrinogen were examined. The data were analyzed by mathematical modeling using nonlinear least-squares curve fitting. Analyses of molecular species present showed that up to 4 mol of either Glu- or Lys-plasminogen was associated with each mol of fibrinogen. For the binding of Glu-plasminogen, the values of the successive changes of the standard free energy of association were found to be -5.48, -7.73, -8.88, and -11.41 kcal/mol (Ka = 2.16 X 10(4), 1.32 X 10(6), 1.06 X 10(7), and 1.08 X 10(9) M-1). For the experimental conditions used here, the association of Lys-plasminogen appears to be described by virtually the same fitting parameters. The very marked cooperativity of association found here would appear to imply that there are significant alterations of the structure of fibrinogen as a result of each successive molecule of plasminogen bound.     

Explicit formulation of titration models for isothermal titration calorimetry

Isothermal titration calorimetry (ITC) produces a differential heat signal with respect to the total titrant concentration. This feature gives ITC excellent sensitivity for studying the thermodynamics of complex biomolecular interactions in solution. Currently, numerical methods for data fitting are based primarily on indirect approaches rooted in the usual practice of formulating biochemical models in terms of integrated variables. Here, a direct approach is presented wherein ITC models are formulated and solved as numerical initial value problems for data fitting and simulation purposes. To do so, the ITC signal is cast explicitly as a first-order ordinary differential equation (ODE) with total titrant concentration as independent variable and the concentration of a bound or free ligand species as dependent variable. This approach was applied to four ligand-receptor binding and homotropic dissociation models. Qualitative analysis of the explicit ODEs offers insights into the behavior of the models that would be inaccessible to indirect methods of analysis. Numerical ODEs are also highly compatible with regression analysis. Since solutions to numerical initial value problems are straightforward to implement on common computing platforms in the biochemical laboratory, this method is expected to facilitate the development of ITC models tailored to any experimental system of interest.     

Understanding microscopic binding of macrophage migration inhibitory factor with phenolic hydrazones by molecular docking, molecular dynamics simulations and free energy calculations

Macrophage migration inhibitory factor (MIF), an immunoregulatory protein, is a potential target for a number of inflammatory diseases. In the current work, the interactions between MIF and a series of phenolic hydrazones were studied by molecular docking, molecular dynamics (MD) simulations, binding free energy calculations, and binding energy decomposition analysis to determine the structural requirement for achieving favorable biological activity of phenolic hydrazones. First, molecular docking was used to predict the binding modes of inhibitors in the binding site of MIF. The good correlation between the predicted docking scores and the experimental activities shows that the binding conformations of the inhibitors in the active site of MIF are well predicted. Moreover, our results suggest that the flexibility of MIF is essential in ligand binding process. Then, MD simulations and MM/GBSA free energy calculations were employed to determine the dynamic binding process and compare the binding modes of the inhibitors with different activities. The predicted binding free energies given by MM/GBSA are not well correlated with the experimental activities for the two subsets of the inhibitors; however, for each subset, a good correlation between the predicted binding free energies and the experimental activities is achieved. The MM/GBSA free energy decomposition analysis highlights the importance of hydrophobic residues for the MIF binding of the studied inhibitors. Based on the essential factors for MIF-inhibitor interactions derived from the theoretical predictions, some derivatives were designed and the higher inhibitory activities of several candidates were confirmed by molecular docking studies. The structural insights obtained from our study are useful for designing potent inhibitors of MIF.     

Mean-field analysis of protein-protein interactions

Calculations were performed on the D1.3-E5.2 antibody-antibody complex estimating the binding affinities of the wild-type and 16 alanine substitutions. Analyzed were structural models of the interfacial region containing a zinc ion and crystallographic waters. A continuum approach was used to evaluate the electrostatic free energies and the hydrophobic effect was calculated by employing a buried molecular surface area relationship. Estimates of the absolute binding affinity reproduced the experimental value within the uncertainty of assessing entropic and strain energy contributions. The best correlation for mutants with experimental data was achieved when the hydrophilicity of created cavities were considered, and yielded a correlation coefficient of 0.7 and an average error of +/-1.4 kcal/mol. Empirically fitting the free energy function produced a smaller error of +/-1.0 kcal/mol. Depending on the electrical potential and electrostatic reorganization, scaling the 'protein dielectric constant' to approximately 10 may improve the accuracy of continuum models for evaluating amino acid substitutions.     

Molecular mass determination by sedimentation velocity experiments and direct fitting of the concentration profiles.

A new method for the direct molecular mass determination from sedimentation velocity experiments is presented. It is based on a nonlinear least squares fitting procedure of the concentration profiles and simultaneous estimation of the sedimentation and diffusion coefficients using approximate solutions of the Lamm equation. A computer program, LAMM, was written by using five different model functions derived by Fujita (1962, 1975) to describe the sedimentation of macromolecules during centrifugation. To compare the usefulness of these equations for the analysis of hydrodynamic results, the approach was tested on data sets of Claverie simulations as well as experimental curves of some proteins. A modification for one of the model functions is suggested, leading to more reliable sedimentation and diffusion coefficients estimated by the fitting procedure. The method seems useful for the rapid molecular mass determination of proteins larger than 10 kDa. One of the equations of the Archibald type is also suitable for compounds of low molecular mass, probably less than 10 kDa, because this model function requires neither the plateau region nor a meniscus free of solute.     

Ligand binding by antibody IgE Lb4: assessment of binding site preferences using microcalorimetry, docking, and free energy simulations.

Antibody IgE Lb4 interacts favorably with a large number of different compounds. To improve the current understanding of the structural basis of this vast cross-reactivity, the binding of three dinitrophenyl (DNP) amino acids (DNP-alanine, DNP-glycine, and DNP-serine) is investigated in detail by means of docking and molecular dynamics free energy simulations. Experimental binding energies obtained by isothermal titration microcalorimetry are used to judge the results of the computational studies. For all three ligands, the docking procedure proposes two plausible subsites within the binding region formed by the antibody CDR loops. By subsequent molecular dynamics simulations and calculations of relative free energies of binding, one of these subsites, a tyrosine-surrounded pocket, is revealed as the preferred point of complexation. For this subsite, results consistent with experimental observations are obtained; DNP-glycine is found to bind better than DNP-serine, and this, in turn, is found to bind better than DNP-alanine. The suggested binding mode makes it possible to explain both the moderate binding affinity and the differences in binding energy among the three ligands.     

Maximum likelihood estimation of molecular motor kinetics from staircase dwell-time sequences

Molecular motors, such as kinesin, myosin, or dynein, convert chemical energy into mechanical energy by hydrolyzing ATP. The mechanical energy is used for moving in discrete steps along the cytoskeleton and carrying a molecular load. High resolution single molecule recordings of motor steps appear as a stochastic sequence of dwells, resembling a staircase. Staircase data can also be obtained from other molecular machines such as F1 -ATPase, RNA polymerase, or topoisomerase. We developed a maximum likelihood algorithm that estimates the rate constants between different conformational states of the protein, including motor steps. We model the motor with a periodic Markov model that reflects the repetitive chemistry of the motor step. We estimated the kinetics from the idealized dwell-sequence by numerical maximization of the likelihood function for discrete-time Markov models. This approach eliminates the need for missed event correction. The algorithm can fit kinetic models of arbitrary complexity, such as uniform or alternating step chemistry, reversible or irreversible kinetics, ATP concentration and mechanical force-dependent rates, etc. The method allows global fitting across stationary and nonstationary experimental conditions, and user-defined a priori constraints on rate constants. The algorithm was tested with simulated data, and implemented in the free QuB software.     

The nature of cation-pi binding: interactions between tetramethylammonium ion and benzene in aqueous solution.

A combined quantum mechanical and molecular mechanical Monte Carlo simulation method was used to determine the free energy of binding between tetramethylammonium ion (TMA+) and benzene in water. The computed free energy as a function of distance (the potential of mean force) has two minima that represent contact and solvent-separated complexes. These species are separated by a broad barrier of about 3 kJ/mol. The results are in good accord with experimental data and suggest that TMA+ binds to benzene more favorably than to chloride ion, with an association constant of about 0.8 M-1.     

Computational analysis of plasmepsin IV bound to an allophenylnorstatine inhibitor

The plasmepsin proteases from the malaria parasite Plasmodium falciparum are attracting attention as putative drug targets. A recently published crystal structure of Plasmodium malariae plasmepsin IV bound to an allophenylnorstatine inhibitor [Clemente, J.C. et al. (2006) Acta Crystallogr. D 62, 246-252] provides the first structural insights regarding interactions of this family of inhibitors with plasmepsins. The compounds in this class are potent inhibitors of HIV-1 protease, but also show nM binding affinities towards plasmepsin IV. Here, we utilize automated docking, molecular dynamics and binding free energy calculations with the linear interaction energy LIE method to investigate the binding of allophenylnorstatine inhibitors to plasmepsin IV from two different species. The calculations yield excellent agreement with experimental binding data and provide new information regarding protonation states of active site residues as well as conformational properties of the inhibitor complexes.     

Validation of an automated procedure for the prediction of relative free energies of binding on a set of aldose reductase inhibitors

Among the available methods for predicting free energies of binding of ligands to a protein, the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) and molecular mechanics generalized Born surface area (MM-GBSA) approaches have been validated for a relatively limited number of targets and compounds in the training set. Here, we report the results of an extensive study on a series of 28 inhibitors of aldose reductase with experimentally determined crystal structures and inhibitory activities, in which we evaluate the ability of MM-PBSA and MM-GBSA methods in predicting binding free energies using a number of different simulation conditions. While none of the methods proved able to predict absolute free energies of binding in quantitative agreement with the experimental values, calculated and experimental free energies of binding were significantly correlated. Comparing the predicted and experimental ΔG of binding, MM-PBSA proved to perform better than MM-GBSA, and within the MM-PBSA methods, the PBSA of Amber performed similarly to Delphi. In particular, significant relationships between experimental and computed free energies of binding were obtained using Amber PBSA and structures minimized with a distance-dependent dielectric function. Importantly, while free energy predictions are usually made on large collections of equilibrated structures sampled during molecular dynamics in water, we have found that a single minimized structure is a reasonable approximation if relative free energies of binding are to be calculated. This finding is particularly relevant, considering that the generation of equilibrated MD ensembles and the subsequent free energy analysis on multiple snapshots is computationally intensive, while the generation and analysis of a single minimized structure of a protein–ligand complex is relatively fast, and therefore suited for high-throughput virtual screening studies. At this aim, we have developed an automated workflow that integrates all the necessary steps required to generate structures and calculate free energies of binding. The procedure is relatively fast and able to screen automatically and iteratively molecules contained in databases and libraries of compounds. Taken altogether, our results suggest that the workflow can be a valuable tool for ligand identification and optimization, being able to automatically and efficiently refine docking poses, which sometimes may not be accurate, and rank the compounds based on more accurate scoring functions.     

Molecular dynamics study of the binding of phenylalanine stereoisomers to thermolysin

The stereospecificity in binding of phenylalanine as inhibitor in the active site of the thermolysin, has been investigated by means of molecular dynamics simulations using free energy integration techniques. The difference in the free energy of binding was found to be 2.0 ± 1.8 kJ/mol in favour of the D-form. This agrees with the experimental value, 2.8 kJ/mol. The result was obtained using a standard empirical force field (that of GROMOS). A different force field with 30% bigger charges (more like ab initio charges) was also tried. This resulted in less fluctuations and a more precise binding, but in a free energy. difference that was clearly larger than the experimental one. The phenylalanine backbone is located close to the zinc atom and the ring stays in the hydrophobic pocket in both the cases. The two stereoisomers differ mainly in the orientation of the backbone plane with respect to the active site and the rotational state of the dihedral around the C---CO^β bond.     

A comparison of weight average and direct boundary fitting of sedimentation velocity data for indefinite polymerizing systems

Analysis of sedimentation velocity data for indefinite self-associating systems is often achieved by fitting of weight average sedimentation coefficients (s(20,w)) However, this method discriminates poorly between alternative models of association and is biased by the presence of inactive monomers and irreversible aggregates. Therefore, a more robust method for extracting the binding constants for indefinite self-associating systems has been developed. This approach utilizes a set of fitting routines (SedAnal) that perform global non-linear least squares fits of up to 10 sedimentation velocity experiments, corresponding to different loading concentrations, by a combination of finite element simulations and a fitting algorithm that uses a simplex convergence routine to search parameter space. Indefinite self-association is analyzed with the software program isodesfitter, which incorporates user provided functions for sedimentation coefficients as a function of the degree of polymerization for spherical, linear and helical polymer models. The computer program hydro was used to generate the sedimentation coefficient values for the linear and helical polymer assembly mechanisms. Since this curve fitting method directly fits the shape of the sedimenting boundary, it is in principle very sensitive to alternative models and the presence of species not participating in the reaction. This approach is compared with traditional fitting of weight average data and applied to the initial stages of Mg(2+)-induced tubulin self-associating into small curved polymers, and vinblastine-induced tubulin spiral formation. The appropriate use and limitations of the methods are discussed.     

Evaluating the energetics of empty cavities and internal mutations in proteins.

The energetics of cavity formation in proteins is evaluated with two different approaches and results are analyzed and compared to experimental data. In the first approach, free energy of cavity formation is extracted by RMS fitting from the distribution of numbers of cavities, N, with different volumes, Vcav, in 80 high-resolution protein structures. It is assumed that the distribution of number of cavities according to their volume follows the Boltzmann law, N(Vcav) = exp [(-a.Vcav-b)/kT], or its simplified form. Specific energy cost of cavity formation, a, extracted by RMS fitting from these distributions is compared to a values extracted from experimental free energies of cavity formation in T4 lysozyme fitted to similar expressions. It is found that fitting of both sets of data leads to similar magnitudes and uncertainties in the calculated free energy values. It is shown that Boltzmann-like distribution of cavities can be derived for a simple model of an equilibrium interconversion between mutants in an extracellular system. We, however, suggest that a partitioning into cavity-dependent and cavity-independent terms may lose meaning when one attempts to describe mutation effects on protein stability in terms of specific free energy contributions. As an alternative approach, a direct molecular mechanics evaluation is attempted of T4 lysozyme destabilization by five single cavity-creating mutations. The calculations are based on the approach used in calculations of the energetics of packing defects in crystals. For all mutations calculated destabilizations agree with the corresponding experimental values within +/-0.6 kcal/mol. A computational relaxation of the mutant was most difficult to achieve for the mutation producing the smallest cavity. However, calculations do not always reproduce crystallographically observed contraction/expansion of cavities. It is suggested that this may be related to usually observed large RMS differences (> 1 A) between crystallographic and energy-minimized protein structures, and thus correct energetics might be easier to calculate than the correct geometry.     

Probing ligand protein binding equilibria with fluorescence fluctuation spectroscopy

We examine the binding of fluorescent ligands to proteins by analyzing the fluctuation amplitude g(0) of fluorescence fluctuation experiments. The normalized variance g(0) depends on the molecular brightness and the concentration of each species in the sample. Thus a single g(0) measurement is not sufficient to resolve individual species. Titration of the ligand with protein establishes the link between molecular brightness and concentration by fitting g(0) to a binding model and allows the separation of species. We first apply g(0) analysis to binary dye mixtures with brightness ratios of 2 and 4 to demonstrate the feasibility of this technique. Next we consider the influence of binding on the fluctuation amplitude g(0). The dissociation coefficient, the molecular brightness ratio, and the stochiometry of binding strongly influence the fluctuation amplitude. We show that proteins with a single binding site can be clearly differentiated from proteins with two independent binding sites. The binding of fluorescein-labeled digoxigenin to a high-affinity anti-digoxin antibody was studied experimentally. A global analysis of the fluctuation amplitude and the fluorescence intensity not only recovered the dissociation coefficient and the number of binding sites, but also revealed the molecular heterogeneity of the hapten-antibody complex. Two species were used to model the molecular heterogeneity. We confirmed the molecular heterogeneity independently by fluorescence lifetime experiments, which gave fractional populations and molecular brightness values that were virtually identical to those of the g(0) analysis. The identification and characterization of molecular heterogeneity have far-reaching consequences for many biomolecular systems. We point out the important role fluctuation experiments may have in this area of research.     

Curcumin binds in silico to anti-cancer drug target enzyme MMP-3 (human stromelysin-1) with affinity comparable to two known inhibitors of the enzyme

In silico interaction of curcumin with the enzyme MMP-3 (human stromelysin-1) was studied by molecular docking using AutoDock 4.2 as the docking software application. AutoDock 4.2 software serves as a valid and acceptable docking application to study the interactions of small compounds with proteins. Interactions of curcumin with MMP-3 were compared to those of two known inhibitors of the enzyme, PBSA and MPPT. The calculated free energy of binding (ΔG binding) shows that curcumin binds with affinity comparable to or better than the two known inhibitors. Binding interactions of curcumin with active site residues of the enzyme are also predicted. Curcumin appears to bind in an extendended conformation making extensive VDW contacts in the active site of the enzyme. Hydrogen bonding and pi-pi interactions with key active site residues is also observed. Thus, curcumin can be considered as a good lead compound in the development of new inhibitors of MMP-3 which is a potential target of anticancer drugs. The results of these studies can serve as a starting point for further computational and experimental studies.     

An efficient computational method for calculating ligand binding affinities

Virtual compound screening using molecular docking is widely used in the discovery of new lead compounds for drug design. However, the docking scores are not sufficiently precise to represent the protein-ligand binding affinity. Here, we developed an efficient computational method for calculating protein-ligand binding affinity, which is based on molecular mechanics generalized Born/surface area (MM-GBSA) calculations and Jarzynski identity. Jarzynski identity is an exact relation between free energy differences and the work done through non-equilibrium process, and MM-GBSA is a semimacroscopic approach to calculate the potential energy. To calculate the work distribution when a ligand is pulled out of its binding site, multiple protein-ligand conformations are randomly generated as an alternative to performing an explicit single-molecule pulling simulation. We assessed the new method, multiple random conformation/MM-GBSA (MRC-MMGBSA), by evaluating ligand-binding affinities (scores) for four target proteins, and comparing these scores with experimental data. The calculated scores were qualitatively in good agreement with the experimental binding affinities, and the optimal docking structure could be determined by ranking the scores of the multiple docking poses obtained by the molecular docking process. Furthermore, the scores showed a strong linear response to experimental binding free energies, so that the free energy difference of the ligand binding (ΔΔG) could be calculated by linear scaling of the scores. The error of calculated ΔΔG was within ≈±1.5 kcal•mol⁻¹ of the experimental values. Particularly, in the case of flexible target proteins, the MRC-MMGBSA scores were more effective in ranking ligands than those generated by the MM-GBSA method using a single protein-ligand conformation. The results suggest that, owing to its lower computational costs and greater accuracy, the MRC-MMGBSA offers efficient means to rank the ligands, in the post-docking process, according to their binding affinities, and to compare these directly with the experimental values.     

Evaluation and Minimization of Uncertainty in ITC Binding Measurements: Heat Error,Concentration Error,Saturation, and Stoichiometry

BackgroundIsothermal titration calorimetry (ITC) is uniquely useful for characterizing binding thermodynamics, because it straightforwardly provides both the binding enthalpy and free energy. However, the precision of the results depends on the experimental setup and how thermodynamic results are obtained from the raw data.MethodsExperiments and Monte Carlo analysis are used to study how uncertainties in injection heat and concentration propagate to binding enthalpies in various scenarios. We identify regimes in which it is preferable to fix the stoichiometry parameter, N, and evaluate the reliability of uncertainties provided by the least squares method.ResultsThe noise in the injection heat is mainly proportional in character, with ~ 1% and ~ 3% uncertainty at 27C and 65C, respectively; concentration errors are ~ 1%. Simulations of experiments based on these uncertainties delineate how experimental design and curve fitting methods influence the uncertainty in the final results.ConclusionsIn most cases, experimental uncertainty is minimized by using more injections and by fixing N at its known value. With appropriate technique, the uncertainty in measured binding enthalpies can be kept below ~ 2% under many conditions, including low C values.General SignificanceWe quantify uncertainties in ITC data due to heat and concentration error, and identify practices to minimize these uncertainties. The resulting guidelines are important when ITC data are used quantitatively, such as to test computer simulations of binding. Reproducibility and further study are supported by free distribution of the new software developed here.