Study of SQ109 analogs binding to mycobacterium MmpL3 transporter using MD simulations and alchemical relative binding free energy calculations

Conformational analysis of monoprotonated ethylenediamine unit in SQ109 (1a)

The ethylenediamine unit in SQ109 (1a) involves three dihedrals affecting the drug’s orientation inside the MmpL3 pore through rotation around (2-Ad)NHCH2–CH2NH2+Ger, (2-Ad)NHCH2CH2–NH2+Ger and C1,AdC2,Ad–NHCH2CH2NH2+Ger bonds (where C1,Ad or C2,Ad are adamantyl carbons at position-1, or -2 in 2-adamantyl group) described in  Tables 13. Using the B3LYP functional with dispersion interactions correction (B3LYP-D3) [36,37,38,39] and the 6-31G(d,p) basis set calculations we performed full geometry optimization and calculated the free energies of the conformational minima of SQ109 (1a), by manual rotation around these three dihedrals in a hydrophilic and a hydrophobic environment. Dispersion correction improves the calculation of the forces acting on the atoms in distances > 3 Å and the accuracy of relative conformational energies calculation which are shown in Table 1. The hydrophilic environment was simulated with an implicit water environment and a dielectric constant (ε) = 80 and the hydrophobic environment was simulated with an implicit chloroform environment and ε = 4.8 using the polarizable continuum model (PCM) [41] and taking advantage of a smooth switching function [42].

Table 1 Calculations of conformational free energies by rotation of (2-Ad)CH2–CH2NH2+Ger in SQ109 (1a) using the B3LYP-D3 /6-31G(d,p) and PCM for water (ε = 80) or chloroform (ε = 4.8)Table 2 Calculations of conformational energies by rotation of (2-Ad)CH2CH2–NH2+Ger in SQ109 (1a) using the B3LYP-D3 /6-31G(d,p) and PCM for water (ε = 80) or chloroform (ε = 4.8)Table 3 Calculations of conformational energies by rotation around C1,AdC2,Ad–NHCH2 in SQ109 (1a) using the B3LYP-D3 /6-31G(d,p) and PCM for water (ε = 80) and chloroform (ε = 4.8)

The rotation around the carbon–carbon bond, (2-Ad)NHCH2–CH2NH2+Ger, in SQ109 (1a), (Fig. 2A) generated the gauche(-), gauche( +), anti and eclipsed conformations and changed the relative orientation of the (2-Ad)NH and NH2+Ger groups and the hydrogen bonding profile of the ethylenediamine unit inside the MmpL3 receptor. The gauche conformations had the lowest free energy following the anti conformation, which was seriously destabilized, while the eclipsed conformer corresponds to the rotation barrier and was not stable. Both gauche conformations were stabilized since a hydrogen bond was formed between protonated and unprotonated nitrogen atoms in the ethylenediamine unit.

Fig. 2figure 2

Important conformational features of SQ109 (1a). In (A) are shown the conformations gauche(-) or gauche( +) or anti by rotation around the ethylenediamine’s unit carbon–carbon bond, (2-Ad)CH2–CH2NH2+CH2Ger (see Table 1)), which defined the relative orientation of 2-adamantyl and geranyl groups. The conformation around the (2-Ad)CH2CH2–NH2+CH2Ger dihedral was anti. In (B) are shown the conformations generated by rotation around the C1,AdC2,Ad–NHCH2 dihedral (see Table 3) which defined the position of axial NHCH2 as regards the cyclohexane subunit of 2-adamantyl group (geranyl group is shown with a light blue sphere); in the structure shown, in the right-hand part, the axial NHCH2 brings CH2 above the cyclohexane subunit of 2-adamantyl increasing steric repulsion (shown with two cross circular arc lines)

The rotation around (2-Ad)CH2CH2–NH2+CH2Ger dihedral favored the anti orientation (Table 2) since in the gauche conformation the steric energy increased due to repulsion between the geranyl and 2AdNH2+ groups. In these gauche conformations the geranyl group of the SQ109 (1a) could not adopt an extended structure that fits inside the narrow pore of the transporter.

The rotation around the C1,AdC2,Ad–NHCH2 dihedral defined the position of the alkyl and the 2-adamantyl ammonium groups and the DFT calculations results are shown in Table 3. The C1,AdC2,Ad–NHCH2 moiety can adopt only the equivalent two gauche conformations (which place the nitrogen of the ammonium group at C2-adamantyl position and Cn-1 at symmetrical positions) due to the severe steric repulsions of the axial NHCH2 in the anti conformation with the cyclohexane subunit group of the 2-adamantyl group (Fig. 2B).

In summary, the DFT study showed that rotation of SQ109 around (2-Ad)NHCH2–CH2NH2+Ger dihedral favored two gauche conformations as minima while other conformers were unpopulated in both dielectric media as well as inside the transporter’s pore demonstrated by our MD simulations (see the next paragraph). Rotation around (2-Ad)NHCH2CH2–NH2+Ger bond favored an anti orientation and rotation around (2-Ad)–NHCH2CH2NH2+Ger bond showed that the position of NHCH2 group above the cyclohexane subunit was prohibited.

MD simulations of the SQ109–MmpL3 complex

MmpL3 protein consists by 12 transmembrane (ΤΜ) α-helices and the TM region also contains two extra α-helices attached to the cytoplasmic membrane surface. We used the structure with PDB ID 6AJG [17] of the protein after excluding the C-terminus that has residues M1-E749. The TM domain consists by the following α-helices and their residues: TM1 (14–33), TM2 (174–199), TM3 (208–224), TM4 (238–264), TM5 (271–301), TM6 (306–338), TM7 (396–415), TM8 (552–576), TM9 (583–601), TM10 (625–648), TM11 (660–690), TM12 (697–728).

In the PDB ID 6AJG, [17] SQ109 (1a) binds to the transmembrane domain of the MmpL3 transporter, from G641 to F649. It disrupts the hydrogen bonding interactions between the two D-Y pairs, where proton translocation takes place, blocking activation of the transporter. In particular, the ethylenediamine group of SQ109 (1a) intervenes between the D256-Y646 and D645-Y257 pairs [14] and forms hydrogen bonds with Asp645. The X-ray structure (PDB ID 6AJG) [17] showed that the complex is stabilized through numerous van der Waals interactions between the geranyl-ethylenediamine moiety of SQ109 (1a) and surrounding amino acids of MmpL3, including I249, I253, I297 in the upper part and L642, Y257, Y646, L686 in the bottom part of the pocket while the 2-adamantyl group lies close to F260 and F649. Ethylenediamine molecule has pKa,1 = 10.71 and pKa,2 = 7.56 [43] and at pH 7.4 the mono and diprotonated species will be equally populated. However, since the basic amino group close to adamantyl group is more hindered to protonation, we assumed that the monoprotonated species must be predominated inside the hydrated MmpL3 pore, although actual protonation states are not known and would probably require neutron diffraction structures. Strikingly, in the X-ray structure (PDB ID 6AJG) [17] the ethylenediamine unit adopts a high energy conformation by rotation around its central carbon‒carbon bond, described in (2-Ad)NHCH2‒CH2NH2+Ger, having eclipsed the C-N bonds. Protein − ligand coordinates obtained from crystallography or cryo-em experiments provide a static model corresponding to a snapshot. While is not uncommon to find ligand conformers in complex with a protein that differ significantly from the lowest energy conformation in solution, due to stabilizing interactions inside the receptor, it is appropriate to judge the stability of the ligand’s conformation in the experimental structure using MD simulations, as previously suggested [44,45,46,47,48].

We performed 100 ns-MD simulations (two repeats, see Fig. S3) using the SQ109 (1a)-MmpL3 complex (PDB ID 6AJG) [17] embedded in a hydrated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer using the amber99sb force field (ff99sb) [49] and a 500 ns-MD simulation (Fig. S4) for testing the stability of the protein complex in longer time scale. For each simulation, initial atom velocities were assigned randomly and independently. Structure of SQ109 (1a)-MmpL3 complex with PDB ID 6AJG [17] has ~ 2.6 Å resolution which is higher than other structures resolved in refs [21, 22] while no other structure of MmpL3 in complex with an SQ109 (1a) analog has been resolved.

We observed a rotation of the carbon–carbon bond in the ethylenediamine unit that converts the eclipsed conformation of SQ109 (1a) in the MmpL3 X-ray structure (PDB ID 6AJG), [17] into the two possible gauche conformations, with no population of the eclipsed conformer over the whole trajectory. In each of the two gauche conformations the monoprotonated ethylenediamine unit formed stabilizing hydrogen bond interactions with both D256 and D645, as compared to the eclipsed conformation which can form hydrogen bonds only with D645. Thus, in the gauche conformation the NH2+ group of the monoprotonated ethylenediamine unit formed direct ionic hydrogen bonds with both D256 and D645 while the unprotonated NH acted as a donor in one direct and one water-mediated hydrogen bond with D256 (Fig. 3). We observed also that water molecules entered the pore forming hydrogen bonds between the monoprotonated ethylenediamine unit and side chains of Y257/Y646, D256/D645, S293 or carbonyl groups backbone peptide bonds. The geranyl-ethylenediamine moiety was tightly enclosed in a narrow area of the transporter’s pore with the geranyl chain surrounded by  amino acids Y257, I297, L642, Y646 while at the bottom wider part of the binding area the 2-adamantyl group had van der Waals contacts with F260 and F649.

Fig. 3figure 3

A Structure of MmpL3 in complex with SQ109 (1a) after 100 ns MD simulations with ff99sb [49]. B Close-up view of the binding site. SQ109 (1a) formed hydrogen bond interactions with the two D-Y pairs, D256-Y646 and D645-Y257. Color scheme for frames: Ligand = petrol sticks, receptor = purple ribbons, residues as light purple sticks. Hydrogen bonding interactions = green dashes. For the protein, the experimental structure of SQ109 (1a) in complex with MmpL3 (PDB ID 6AJG [17]) from Ms, was used as the starting structure for the MD simulations after excluding C-terminus, consisting by M1-E749 residues. The TM region included the following helices and their residues: TM1 (14–33), TM2 (174–199), TM3 (208–224), TM4 (238–264), TM5 (271–301), TM6 (306–338), TM7 (396–415), TM8 (552–576), TM9 (583–601), TM10 (625–648), TM11 (660–690), TM12 (697–728)

Simulations of the complexes of MmpL3 with SQ109 analogs 1b-i, 2Docking calculations results

The X-ray structure of the MmpL3 − SQ109 (1a) complex (PDB ID 6AJG [17]) was used as the template structure for the docking calculations of ligands SQ109(1a), 1b-i, 2 with MmpL3 after excluding C-terminus, consisting of M1-E749 residues as previously described. The highest score docking pose of SQ109 (1a) inside the MmpL3 pore (ChemScore [50] scoring function) had a root-mean-square-deviations of heavy atoms (RMSDligand) 1.7 Å compared to the structure with PDB ID 6AJG [17]. This suggested that the calculated docking poses could describe accurately the binding orientation of SQ109 analogues. The docking poses of the monoprotonated ethylenediamine SQ109 analogues 1b-i, 2 showed that the new adamantyl moieties can fit in the empty region at the bottom of the binding area where for example the bigger alkyl adducts linked at C-2 adamantyl position can be accommodated. The docking calculations showed that addition of alkyl, aryl or heteroaryl adducts, e.g. Me, Et, nPr, nBu, Ph, Bn, Hex, 5-phenyl thiazol-2-yl (Ph-Thz) at the adamantyl C-2 of SQ109 (1a), in compounds 1b-i, respectively, or replacement of the 2-adamantyl group by a 1-adamantyl-dimethylmethylene (C-1 dimethylmethylene) group in compound 2, resulted in highest score docking poses having a gauche(−) or gauche( +) by rotation of (2-Ad)CH2–CH2NH2+CH2Ger that bind MmpL3 pore.

The docking algorithm produced 30 docking solutions. We visually inspected and realized that the first 4–5 highest docking solutions corresponded to a similar conformation of the ligand inside the receptor.

For the diprotonated SQ109 (1a) analogs, mainly the gauche(−) or gauche( +) by rotation of (2-Ad)CH2–CH2NH2+CH2Ger were obtained as highest score docking poses, but also in two cases the anti and eclipsed conformations were also observed (Table S2), since these later conformers stabilized inside the receptor with strong hydrogen bonding interactions despite their much lower stability in solution compared to the gauche conformation (Table 1). It was not unusual that the highest docking pose for each ligand was different since this reflected the high flexibility of the ligands and the random nature of the docking algorithm (genetic algorithm runs). For the MD simulations we used the highest scoring docking pose as starting conformation. We also tested as starting structure a conformation for the ligands 1b-i resulted from superposition with the last snapshot from 100 ns-MD simulation of MmpL3-SQ109 (1a) complex.

MD Simulations of SQ109 analogs in complex with MmpL3

Το explore the dynamic behavior between the ligands 1a-i, 2 and MmpL3 we performed 100 ns-MD simulations (two repeats, see Fig. S3) with starting structure the docking pose with the highest score. Thus, these docking poses were embedded in hydrated POPC bilayers of ~ 140 lipids and subjected to MD simulations using ff99sb. For each simulation, initial atom velocities were assigned randomly and independently.We performed the 100 ns-MD simulations using as starting structure a conformation for the ligands 1b-i resulted from superposition with SQ109 (1a) but the results were similar.

For the sizeable alkyl substituents at C-2 adamantyl position,, e.g. in in analogs 1f, g, i, 500 ns-MD simulations were also performed for testing the stability of the protein complexes in longer time scale (Fig. S4). For each simulation, initial atom velocities were assigned randomly and independently. We observed from the 100 ns-MD simulations or the 500 ns-MD simulations that the RMSD values of the protein TM Cα carbons converged after 10–40 ns with RMSDprotein (Ca TM) ≤ 2.1 Å, suggesting small changes compared to the X-ray structure (Fig. 4, S1, Table 4). A different stabilization between the full protein and the TM region was observed reflected by the different RMSDprotein (Ca full protein) and RMSDprotein (Ca TM) with values for the latter much smaller than the former. It is not unusual that the loops connecting the TM-helices are very flexible and increase the RMSD (Ca full protein). The ligands, which contained very flexible moieties such as the geranyl and ethylenediamine groups, shifted considerably from the starting docking pose, as revealed from the high values of RMSD ligand and in most cases the ligand binding conformation equilibrateds in a stable position inside the transporter’s pore after 70 ns (Fig. 4, S1, Table 4). The last frame described well the ligand–protein interaction frequency plots. Figure 4 shows last frames and ligand–protein interaction frequency plots for SQ109 (1a) and selected alkyl and aryl groups attached at the adamantyl C-2 of 1a (R = H) in Fig. 1 including analogs 1b (R = nBu), 1 h (R = Ph), 1i (R = Ph-Thz) while the last frames, ligand–protein interaction frequencies for the other analogs, as well as the RMSD plots from MD simulations can be found in Fig. S1.

The MD simulations showed that the monoprotonated ethylenediamine unit of the ligands 1b-i, 2 adopted a gauche conformation, which favored hydrogen bonding interactions with 256 and/or D645 of MmpL3 (Fig. 4, S1). According to the MD simulations the complexes between 1b-i, 2 and MmpL3 formed common van der Waals interactions with the protein’s residues along the narrow area of the transporter pore. Compared to SQ109 (1a) the geranyl-ethylenediamine moiety was surrounding similarly by the amino acid residues L642, Y646, Y257 while the 2-adamantyl group lied close to F260 and F649 at the bottom part of the binding area. However, the alkyl adducts increased the hydrophobic interactions at the bottom of the binding area with F260 and F649 and increased also the hydrophobic interactions with residues in the walls of narrow area of the pore, e.g. I253, see for example the analogs 1e, h, i in Fig. 4. Water molecules formed hydrogen bonds between the monoprotonated ethylenediamine unit and receptor’s residues, as described previously for SQ109 (1a).

We explored the disposition of the important residues for binding of the ligands (D256, D645, Y257, Y646, F260, F649) along the MD simulation and checked for other frames that could accommodate better the ligands. The RMSD (Ca) of these important residues for binding were plotted (Fig. S5). The RMSD (Ca) converged to values that range between 0.5 and 1.5 confirming not important disposition of these residues along the MD simulation.

While the monoprotonated form of SQ109 (1a) and analogs most likely predominated inside the hydrated MmpL3 pore, we also performed MD simulations for the diprotonated species (Fig. S2). The MD simulations showed that the ligands in the diprotonated ethylenediamine form had similar coordinates inside the MmpL3 (Fig. S2) compared to the monoprotonated form.

MM-GBSA binding free energy calculations

We applied the MM-GBSA method [29,30,31] using the OPLS2005 force field [51] for the calculation of binding free energies (ΔGeff) of ligands 1a-i, 2 inside the MmpL3 pore, using ensembles from 20 ns-MD simulations and calculated binding free energies without or with considering the membrane environment of the protein complex. For each simulation, initial atom velocities were assigned randomly and independently. Thus, we tested the membrane protein – ligand systems using an implicit membrane, i.e. a hydrophobic slab, [52,53,54,55] and considering explicitly water molecules inside the binding area [56]. The range of molecular weight of the ligands is 400–500 Da which corresponds to tripeptides [30] and their carbon skeleton was long enough to interact with many residues inside the receptor area. If a correct ranking of the binding affinities of 1b-i, 2 for the MmpL3 pore could be accomplished with the MM-GBSA method this would be a significant reduction in computational resources, compared to the accurate binding free energy perturbation methods for 1b-i, 2—MmpL3 complexes in the POPC lipid bilayers containing ~ 105 atoms.

However, the calculated mean values of three repeats for the monoprotonated forms (Table 4) and diprotonated forms (Table S2) showed that the MM-GBSA method (applied with or without modifications to consider the membrane environment of the protein complex using a hydrophobic slab [52,53,54,55] and the explicit waters inside the binding area) [56] did not afford valuable results. The binding free energy values were dispersed and did not follow any trend, see Table 4, S2. Indeed, as a calculation method, MM-GBSA can show large deviations (e.g. 4 kcal mol−1) in standard binding free energies compared to the experimental binding free energies. The method normally can provide useful results [57] for ranking of substituents [29] in the same series of ligands when the experimental binding affinities range is ~ 1000 or higher [30] which is not the case for 1b-i, 2 since their KD values differ only by 4-fold.

Compared to the monoprotonated species, we observed that the diprotonated ethylenediamine moieties in compounds 1a-i,2 form stronger hydrogen bonding interactions with the receptor area. This agreed with the MM-GBSA binding free energy values for the monoprotonated forms compared to the diprotonated forms since in the latter case the values were significantly lower consistent with enhanced hydrogen bonding interactions between the ligands and the receptor’s binding site (Table 4, S2). Because of the stronger hydrogen bonding interactions, the ligands had smaller flexibility as was shown by the lower RMSDprotein (Ca TM) and lower RMSDligand (Fig. S1, S2, Table 4, S2).

Table 4 Ligand-MmpL3 binding free energies (ΔGeff) calculated using the MM-GBSA [29,30,31] method with the OPLS2005 [51] force field for the calculations of the intermolecular interactions without or with using a hydrophobic slab [52,53,54,55] to model the membrane environment of the protein, RMSDligand and RMSDprotein (Ca TMD) mean values from 100 ns-MD simulations for 1a-i, 2Fig. 4figure 4

Last frames (left part) of monoprotonated ethylenediamine form for SQ109 (1a), and analogs 1e, h, i inside the MmpL3 transporter in a POPC lipid bilayer from 100 ns-MD simulations with ff99sb. The receptor-ligand interaction frequency histograms plotted for the whole trajectory are shown on the right (the alkyl group connected at the adamantyl C-2 is shown inside the parenthesis following the compound number. Color scheme for frames: Ligand = petrol or purple or brown or orange sticks, receptor = white ribbons, residues in light purple sticks, hydrogen bonding interactions = dark grey dashes. For the protein, the experimental structure of SQ109 (1a) in complex with MmpL3 (PDB ID 6AJG 17) was used as the reference structure for the MD simulations after excluding C-terminus which included M1-E749 residues

Structure–activity relationships of SQ109 analogs against MmpL3 using alchemical binding free energy calculations with TI/MD

We previously measured [13] using SPR the binding affinities of SQ109 (1a) and the 9 active analogs 1b-i, 2 against MtMmpL3 [23] with Kd values R = H (SQ109 (1a), 2060 μM); R = Me (1b, 248 μM); R = Et (1b, 190 μM); R = nPr (1d, 106 μM); R = nBu (1e, 108 μM); R = nHex (1d, 81 μM). It was observed that the Kd decreased showing tighter binding to MmpL3 as the R substituent at C-2 adamantyl (which is an H in SQ109) became larger and more hydrophobic. A similar effect was observed with phenyl (1 h, 136 μM); benzyl (1 g, 74 μM) (Table S1). The C-1 dimethylmethylene analog 2 had a Kd = 106 μΜ which was close to the isomeric 1d (n-Pr) which had a Kd = 120 μM and 1i (Ph-Thz) with a Kd = 91 μM.

The FEP/MD [58,59,60] or TI/MD [61, 33, 34] methods which are based on statistical mechanics can provide accurate results for relative binding free energies with an error 1 kcal mol−1 and have been applied in membrane protein–ligand complexes, e.g. GPCRs. [32, 35, 62,63,64,65]We applied the TI/MD method combined with a thermodynamic cycle method to examine if the binding profile of the ethylenediamine analogs 1b-i, 2 was the same with SQ109 (1a) in its complex with MmpL3 (PDB ID 6AJG [17]). This might be likely if there is an agreement between calculated and experimental relative binding free energies for alchemical transformations between pairs of compounds 1a-i, 2. We performed TI/MD calculations for alchemical transformations in selected pairs of diamine SQ109 analogs that were not accompanied with large changes in ligand’s alkyl. Thus, we calculated perturbations by one or two methylene groups in the C-2 alkyl adduct (Table 5).

Table 5 Free energy perturbation values computed with the MBAR method [66] for alchemical simulations performed with TI/MD for pairs of ligands bound to Mmpl3

The end states in the alchemical calculations tested were similar to the structure in the corresponding complexes resulted from the converged 100 ns-MD simulations. This was checked to certify that the 2 ns-MD simulation in each λ-state was enough for the complexes to converge in an equilibrium structure. Two repeats of TI/MD calculations were performed for each alchemical transformation.

The effect in binding free energy by increasing the length of the alkyl chain by one methylene, which was examined with the alchemical perturbations 1a (H) → 1b (Me) or 1b (Me) → 1c (Et) or 1c (Et) → 1d (n-Pr) or 1d (n-Pr) → 1e (n-Bu), was to increase binding affinity (Table 5). As noted previously, in the 100 ns-MD simulations of MmpL3—1a-e complexes the RMSDprotein (Ca TMD) was ≤ 2.1 Å (Table 4). Thus, the last snapshots of the MD simulations were suitable starting structures for the TI/MD simulations of the studied perturbations (Fig. 4, S1).

The biggest change in experimental binding free energy was noted when  H (SQ109) changed to Me (1b), ΔΔGb,exp =− 1.30 kcal mol−1 ± 0.79 kcal mol−1, and when Ph (1h) chaged to Bn (1 g), ΔΔGb,exp = − 0.38 ± 0.02 kcal mol−1 or Et (1c) changed to Pr (1d), ΔΔGb,exp =− 0.36 kcal mol−1 ± 0.29 kcal mol−1.

The binding free energy changes for 1a (H) → 1b (Me) were ΔΔGb,exp =− 1.30 kcal mol−1 ± 0.79 kcal mol−1, ΔΔGb,TI/MD = − 0.49 ± 0.06 kcal mol−1, for 1b (Me) → 1c (Et) were ΔΔGb,exp = − 0.16 ± 0.14 kcal mol−1, ΔΔGb,TI/MD = 0.06 ± 0.08 kcal mol−1, for 1c (Et) → 1d (n-Pr) were ΔΔGb,exp = − 0.36 ± 0.29 kcal mol−1, ΔΔGb,TI/MD = − 0.87 ± 0.09 kcal mol−1 and for 1d (n-Pr) → 1e (n-Bu) were ΔΔGb,exp = 0.01 ± 0.30 kcal mol−1, ΔΔGb,TI/MD = 0.20 ± 0.11 kcal mol−1.

We considered next the perturbation in the C-2 adamantyl alkyl by two methylenes in 1e (n-Bu) → 1f (n-Hex) and we studied the change by one methylene from phenyl to benzyl in C-2 substituent in 1h (Ph) → 1g (Bn) where the perturbation in conformational space should be relatively important. The binding free energy changes were ΔΔGb,exp =—0.18 kcal mol−1 ± 0.01 kcal mol−1, ΔΔGb,TI/MD = -1.42 ± 0.15 kcal mol−1 and ΔΔGb,exp = -0.38 ± 0.02 kcal mol−1, ΔΔGb,TI/MD = -1.83 ± 0.15 kcal mol−1, respectively. We did not test larger perturbations that were not consistent with the method’s principles. [61]

留言 (0)

沒有登入
gif