Adequate prediction for inhibitor affinity of Aβ40 protofibril using the linear interaction energy method

The search for efficient inhibitors targeting Aβ oligomers and fibrils is an important issue in Alzheimer's disease treatment. As a consequence, an accurate and computationally cheap approach to estimate the binding affinity for many ligands interacting with Aβ peptides is very important. Here, the calculated binding free energies of 30 ligands interacting with 12Aβ11–40 peptides using the linear interaction energy (LIE) approach are found to be in good correlation with experimental data (R = 0.79). The binding affinities of these complexes are also calculated by using free energy perturbation (FEP) and molecular mechanic/Poisson–Boltzmann surface area (MM/PBSA) methods. The time-consuming FEP method provides results with similar correlation (R = 0.72), whereas MM/PBSA calculations show very low correlation with experimental data (R = 0.27). In all complexes, van der Waals interactions contribute much more than electrostatic interactions. The LIE model, which is much less time-consuming than both the FEP and MM/PBSA methods, opens the door to accurate and rapid affinity prediction of ligands with Aβ peptides and the design of new ligands.


Introduction
Alzheimer's disease (AD) is a common neurodegenerative disorder in the senior population and has strong negative affection on intellectual abilities. 1-3 More than 44 million people worldwide suffer from AD leading to a social welfare burden. 4 Despite extensive biophysical and clinical studies, all drugs targeting Ab oligomers and amyloid plaques made primarily of the 40 (Ab 40 ) and 42 (Ab 42 ) amino acids have failed and there is currently no efficient treatment against AD. [5][6][7][8][9][10] The mechanism of Ab peptide accumulation in the extracellular space with aging and the memory impairment by Ab aggregation are still poorly understood. Based on the cascade amyloid hypothesis that Ab oligomers are the most toxic species, 11 nding efficient compounds that inhibit or retard the formation of Ab oligomers and brils is a common strategy, for instance natural compounds, 12 metal chelators, 13 chaperones, 14 RNA aptamers, 15 and short peptides. 16,17 Computer simulations at different levels of protein representations are an efficient tool to provide deeper insights into the aggregation of Ab peptides. 18 Another goal of computational studies is to accurately estimate the binding affinity of protein-ligand complexes facilitating drug development. 19 In this context, several methodologies have been developed, including free energy perturbation (FEP), 20,21 thermodynamic integration (TI), 22,23 molecular mechanic/Poisson-Boltzmann surface area (MM/PBSA), [24][25][26] linear interaction energy (LIE), [27][28][29][30] and docking approaches. 31,32 Because of the large number of ligands to be tested, one have to consider both the accuracy of the calculated binding affinities and the computational costs. Docking methods use simple scoring functions to estimate binding affinities of a large ligand library with low CPU costs. However, because of the lack of exibility of ligands and proteins and the low accuracy of scoring functions, docking oen gives poor results compared with experimental data. 33,34 On the other hand, FEP/TI method, which calculates absolute binding free energies based on conformational sampling from extensive molecular dynamics (MD) simulations, can provide accurate results, but extremely expensive computational costs have hindered its application to large-scale screening. The end-point MM/ PBSA approach, which combines molecular mechanics and continuum solvents to calculate binding free energies, has been successfully applied for several systems. [35][36][37][38] However, the results of this method strongly depend on several factors, such as continuum-solvation method, dielectric constant, charge model, and entropic calculations. [39][40][41] The MM/PBSA method is also time consuming for calculating the vibrational entropy of large systems using normal mode analysis. Last but not least, the end-point free energy calculation method, LIE, has also been successfully applied to various systems. [42][43][44][45][46][47][48][49] In this method, the binding free energy is calculated based on the average van der Waals and electrostatic interaction differences of ligand with its surrounding environments upon association, i.e. the free ligand in solvent (free statedenoted as subscript f) and the ligand in complex with solvated protein (bound statedenoted as subscript b).
The coefficients a and b are the scaling factors for nonpolar and polar terms, whereas g is a constant term, which correlates to the change of the hydrophobic nature of the binding site according to different types of ligands. 50 In some modied versions, the effects of entropic contributions, 42 solventaccessible surface area (SASA), 43,51 intramolecular energies, 52 and structural descriptors of ligands have also been considered. 47 Moreover, continuum solvent models have been used rather than explicit solvent models in the original LIE method to reduce computational costs. [53][54][55][56] The performances of the LIE and MM/PBSA methods have been compared in several reports. Their accuracies strongly depend on the systems, 25,46,[57][58][59] and LIE is less time consuming than MM/PBSA because it does not require any entropy calculations. 41,60 Computer-aided drug design, which focuses on developing new inhibitors for Ab peptides and elaborates an efficient method for absolute binding free energy calculations for Ab peptide systems, has attracted a lot of attention. [61][62][63][64][65] While in some studies MM/PBSA method provided reasonable results in ranking ligand affinity, the absolute values somewhat did not correlate well with the experimental results. [66][67][68][69] This discrepancy may result from the choice of the dielectric constant and/or the calculation of entropy, which is usually considered for only the last MDgenerated structure. [39][40][41] To the best of our knowledge, LIE has not been used to calculate binding free energy for Ab peptide systems. This could be due to the lack of binding poses for Ab oligomers leading to the difficulty in nding a general parameter set for LIE model. In this work we have calculated the absolute binding free energies of 30 inhibitors interacting with a 12Ab 11-40 oligomer using molecular dynamics (MD) simulations and LIE, and compared the results with experimental data and the binding affinities calculated by using the FEP and MM/PBSA methods.

Materials and methods 12Ab oligomer complexes preparations
The structure of 12Ab  oligomer was taken from the protein data base (PDB ID: 2LMN) based on numerous constraints from solid state NMR and electron microscopy. 70 Available inhibitors for Ab 40 peptide were searched referring the binding database. 71 The 3D structures of 30 inhibitors were downloaded from PubChem database (see Table S1 in the ESI † for the list of inhibitors). The geometrics of these compounds was optimized using chemical quantum calculation with B3LYP functional at 6-31G(d) level. The molecular docking protocol using Autodock Vina 72 was then applied to obtain the Ab bril complexes as starting structures for MD simulations. The binding poses of these inhibitors to the 12Ab 11-40 oligomer are depicted in Fig. 1.

Molecular dynamics (MD) simulations
The GROMACS version 5.1.3 (ref. 73) was employed to simulate the solvated complexes. The Amber99SB-ILDN force eld 74 was used to parameterize the protein and counter ions. The general Amber force eld 75 was used to represent the ligands with the help of Ambertools17, in which the atomic charges were obtained using RESP method 76 with quantum chemical calculations at B3LYP/6-311G(d,p) level. The ligand topologies were then converted to GROMAS format using ACPYPE package. 77 Finally, the complexes were inserted into the dodecahedron periodic boundary condition water boxes using TIP3P water model. 78 The smallest distance between a complex to the boundary of water box was chosen as 1. The soluble complexes were rst minimized using the steepest descent scheme. The minimized congurations were then relaxed in NVT and NPT ensembles with 100 ps MD length per simulations. The complexes were restrained by NVT simulations using a small harmonic force and free of restraints by NPT MD simulations. The relaxed system was then used as initial conformation of MD simulations during 20 ns. A total of four MD simulations for each complex were carried out with the same starting structure and different initial velocities. The MD parameters are described in previous studies, 79,80 in which the non-bonded van der Waals cut-off is 1.0 nm and the PME method is used for electrostatic interactions.

Data analysis and inhibitor affinity calculation
The root-mean-square-deviations (RMSDs) of Ca atoms were calculated with respect to the initial conguration. A hydrogen bond (HB) is dened based on geometric criterion: distance between acceptor A and donor D is less than 3.5Å and the A-H-D angle is >135 . A non-bonded contact (NBC) is dened if the distance between any hydrophobic atom, i.e. C or S atom, of the ligand and any other atom of Ab peptides is in range of 2.9-3.9Å. Diagrams for HB and NBC networks are carried out using LigPlot+ program. 81 The van der Waals and electrostatic interactions were calculated to build prediction models for LIE method (eqn (1)) with linear regression analysis. Free binding energy of inhibitors with 12Ab 11-40 oligomer using the LIE model (eqn (1)) were also compared to MM/PBSA and FEP approaches, successfully applied to several Ab peptides complexes. [64][65][66] Computational details for MM/PBSA and FEP calculations are extensively described in our previous papers. 82,83 Results and discussion

Structures and energies of complexes during MD simulations
The time-evolution of Ca RMSD values of the peptides in the 30 complexes is shown in Fig. S1 in ESI. † The RMSD values increased rapidly at the beginning of simulations and all complexes reached their equilibrium states aer 5 ns. In most of cases, the RMSD values are smaller than 0.8 nm, which clearly indicates that these systems are quite stable during the last 15 ns.
Averaged van der Waals and electrostatic interactions between each ligand and Ab 11-40 oligomer (bound state) or water (free state) as well as the experimental binding affinity of each ligand are listed in Table S2. † Full references and experimental K i values for these inhibitors are provided in ESI. † Except for Thioavin T (Pubchem ID 16954), the magnitude of electrostatic energies of all other ligands in the bound state are smaller than in the free state, indicating that the 12Ab  oligomer mainly interacts with the ligands by van der Waals interactions. To obtain deeper insights into the binding modes of each inhibitor, 2-dimension protein-ligand interaction diagrams using the last MD-generated structures are drawn in Fig. S2. † All diagrams show that strong van der Waals interactions in all complexes, and a small number of average hydrogen bonds in the 4 MD simulations (Fig. S3 †).

Optimization of LIE equations
We aimed at building a LIE model that can exactly predict binding free energies. Based on previous theoretical studies on the complex of P450 CAM , 50 potassium channel, 84 and aspartic proteases; 85 the parameters a was set to be 0.18, and the parameter b was set in the range of 0.33-0.50. 52,86 Using standard values (a ¼ 0.18 and b ¼ 0.5), no correlation (R z À0.57) was found between the calculated and experimental binding free energies. This is probably due to the lack of a specic binding site for Ab oligomers, and/or, as mentioned above, the loss of electrostatic interaction of ligand from the free to bound states, which is not observed in other protein complexes. [42][43][44][45][46][47][48][49] Because of the failure of these standard a and b parameters, a new set of parameters for Ab peptides is needed. Based on the average interaction energies of training set comprising of 20 inhibitors taken randomly listed in Table S4, † the a, b, and g parameters are calculated to have the values of 0.288, À0.049, and À5.880 kcal mol À1 , respectively (eqn (2)) giving a Pearson correlation R and a standard error of 0.79 and of 0.95 kcal mol À1 , respectively (Fig. 2).
In addition,Åqvist and co-workers proposed that for a very small or negative b value, the electrostatic interactions in water and in complex with protein can be scaled differently, giving thus different b b (for bound state) and b f (for free state) values. 52 Although very small or even negative values of b have been found in some studies, 46,87,88 a small b value of À0.049 found in our calculations is still somewhat unexpected. This may result from either the loss of electrostatic interactions of ligands upon association (see Table S2 †) or the effect of functional groups of ligands. 89 Observed results are in good consistent with previous studies that the Ab inhibitors favorably form domination of vdW interaction energy. 90,91 It also is in good consistent with the obtained 2-D protein-ligand interaction diagrams (Fig. S2 †) where the vdW interaction was found to dominate over the electrostatic interaction.
In addition, the obtained g value was calculated of À5.880 kcal mol À1 . The large negative value of g indicates small hydrophobic desolvation effects of ligands and very strong hydrophobic interactions between the ligands and 12Ab  oligomer, leading to the deformation of Ab peptides structures. 50,56 This is consistent with the fact the van der Waals interactions are dominant in the interaction of ligands and 12Ab  oligomer (see Table S2 and additional data of ESI for details †). The difference is probably caused by the erroneousness of imitating the interaction among constituent molecules, including protein, ligand, and solvated molecules. 92,93 Validate the approach We emphasize that our a, b, and g parameters and our Pearson correlations do not change if we apply on the test set comprising of 10 inhibitors taken randomly (Table S5 †), verifying therefore the correlative power of our model and indicating that our parameter set for the LIE model is rather robust. In particular, the Pearson correlation and standard error were obtained as the values of 0.72 and 1.09 kcal mol À1 , respectively (Fig. 3). Absolutely, the model accuracy is appropriate to estimate the binding affinity of trial inhibitor to Ab 40 system. Whereas, the small error means that is able to categorize ligands revealing similar binding affinities.

Comparing LIE model with MM/PBSA and FEP methods
To establish the performance of our LIE calculation with respect to the FEP and MM-PBSA methods, the binding free energies of all ligands interacting with Ab 11-40 using the three methods are listed in Table S3. † The correlations between the calculated and experimental data using the three approaches are shown in Fig. 4. It should be emphasized that our groups have been successfully applied FEP and MM/PBSA to calculated binding free energy for various ligand-Ab peptides complexes. 64,66,91,94,95 When all of ligands are considered, the time-consuming FEP method, expected to give accurate binding affinity of proteinligand complex, generates results with a correlation to experimental data (R ¼ 0.72) similar to the LIE method with our parameter set (R ¼ 0.79). In contrast, the binding affinity calculated using MM/PBSA method gives a poor correlation with experiment data (R ¼ 0.27). In Ab peptide systems, as there is no unique binding pose, the choice of the most appropriate dielectric constant value for accurate electrostatic calculations is problematic. The dependence of the MM/PBSA results on the dielectric constant values has already been discussed out on other systems; however, no systematic approach to improve the results was reported. 41,59 A poor correlation can also come from the entropy estimation. In order to reduce the computational time, we only chose the last structures of MD simulation to calculate the entropy. Further studies are required to gain deeper insights and improve the results of MM/PBSA approach.

CPU time consumption
The required computing resource for ligand affinity prediction to 12Ab 11-40 peptide using LIE approach is much less than that of MM/PBSA and/or FEP protocols. The ratio is up to ca. 50 times. Indeed, the MM/PBSA approach requires a huge of CPU time for estimating DG PB and ÀTDS values due to the 12Ab 40ligand systems consist of more than 5700 atoms. It takes more than 5700 node-hours (>8 months) to estimate the ligandaffinity for a single complexed system using MM/PBSA approach. In particular, it costs ca. 1200 computing nodehours for estimating DG PB value over 150 snapshots of each complex. In another hand, the entropic approximation is required more computational costs, although the normal model calculation is very approximate. The entropic evaluation is serial calculation, whereas the simulated time is up to several weeks for every normal model estimation. However, it is noted that the computational time was not included the MD production time, which is of ca. 3.6 node-days for three independent trajectories of a complex. That also is total time to predict the binding affinity of a ligand to 12Ab 11-40 peptide using LIE approach, because non-covalent bond interaction energy between two molecules was recorded over simulation time intervals. Furthermore, the FEP results were obtained throughout the alchemical process, the CPU time consumption is extremely huge that it is slightly larger than MM/PBSA approach in this case. Fig. 3 Correlation between experimental binding free energies and that calculated using LIE model (eqn (2)) of testing set comprising of 10 complexes taken randomly listed in Table S5. † Fig. 4 Correlations between calculated and experimental binding free energies calculated using different approaches. Three methods were applied on 30 complexes.

Conclusions
In this study, a LIE model has been successfully optimized to predict the binding affinity of ligands to 12Ab  protobril showing good correlation between calculated and experimental values (R ¼ 0.79) with a small error (d ¼ 0.95). We found that the standard LIE parameters cannot provide a good correlation with experimental data, which is likely due to the lack of a unique binding pose for the 12Ab 11-40 oligomer.
The van der Waals interactions are dominant and much more important than electrostatic interaction, which is consistent with the small number of H-bonds found in the complexes of the ligands and 12Ab  peptide during MD simulations. The magnitude of electrostatic interactions of the ligand in the free state is larger than that in the bound state in most of the cases, which indicates small desolvation effects of ligands.
Binding free energies of all complexes are calculated by using the FEP and the MM/PBSA methods. Compared to our LIE model, the results from the FEP approach do not provide better correlation (R ¼ 0.72) with experimental data, while the MM/ PBSA method provides very low correlation (R ¼ 0.27).
An accurate and fast method for computing binding affinities of molecules with Ab peptides is of great interest. Our LIE model, much less time-consuming than both FEP and MM/ PBSA methods, opens the door to accurate and rapid affinity prediction of ligands with Ab peptides and may help design new ligands.

Conflicts of interest
There are no conicts to declare.