Molecular Dynamics in Flare™: The Importance of Ionic Strength when Modeling Charged Ligands

Including a salt concentration to match the experimental conditions of your binding assay to your Molecular Dynamics (MD) or Free Energy Perturbation (FEP) model can significantly improve accuracy¹. Binding affinity predictions can be improved: particularly for charged ligands. Recent work at Cresset has focused on the inclusion of ligand charge change transformations in FEP calculations, now possible in Flare V8, specifically considering the TYK2 dataset². This dataset incorporates two negatively charged ligands: one with a carboxylate (ligand '32') and one with a tetrazole group (ligand '33') into an otherwise neutral set of ligands (Figure 1). Here we will focus on the Molecular Dynamics (MD) investigations of ligand 33.

Figure 1. The two charged ligands from the TYK2 dataset², ligand 32 (left) which has a carboxylate group, and ligand 33 (right) with a charged tetrazole.

Molecular dynamics of charged ligands

Here we investigate the effect on the dynamics of ligand 33 of setting the ionic strength (i.e. NaCl salt concentration) to 150mM compared to simply neutralizing the system (with no added concentration of salts). Within the 'More Options' of the Dynamics calculation panel, there are options for advanced settings, as shown in Figure 2. One choice is the default option to "Just neutralize" the system, while another allows manual adjustment. In this instance, it was set to 150mM to better mimic assay conditions³. The Dynamics experiment for this system was run for 10ns, using three different ways: in solution, in complex set to “Just neutralize”, and in complex with an ionic strength of 150mM.

Figure 2. Options available for molecular dynamics within Flare, specifically highlighting the Solvent ionic strength.

MD analysis in Flare

As the dynamics run completes, an MD trajectory is shown, providing an easy way to visualize the RMSD of various components within the simulation. Checking the ligand RMSD serves as a reliable way to ensure the ligand/s in question are stable in complex and performing as expected, especially in preparation for FEP calculations. For example, large or sudden shifts in RMSD value, demonstrating high deviation from the initial ligand frame, may indicate substantial movements in the ligand or even potential dissociation from the binding site. In this case, comparing the RMSD of ligand 33 with and without the ionic strength set to 150mM, as shown in Figure 3, we can identify frames in which the tetrazole ring is likely flipping within the simulation.

Figure 3. Graphical depiction of the ligand heavy atoms RMSD from two 10ns Dynamics calculations for ligand 33, with ionic strength set to both "Just neutralize" (blue) and 150mM (orange).

In addition to RMSD, Dynamics analysis within Flare provides a breakdown of the number of contacts which occur between the ligand and nearby residues, including water molecules. Focusing on the two MD simulations which include the protein-ligand complex, key binding interactions between ligand 33 and an Arginine residue (ARG13) are identified. By comparing contacts made between the tetrazole ring and ARG13, it was shown that although fewer interactions occurred when an ionic strength of 150mM was applied, these interactions were more stable, as they were present in a higher percentage of MD frames, compared to those made when using the "Just neutralize" function. This behaviour is shown in Table 1, parts A and B.

Table 1. Percentage of frames in which hydrogen bonds and salt bridges were present during the 10ns dynamics simulation between the four tetrazole nitrogens of ligand 33 and the ARG13 protein residue, when (A) "Just neutralize" or (B) 150 mM were implemented for the ionic strength setting.

(A) "Just neutralize"

(B) 150mM Ionic Strength

To explore the ligand RMSD differences and reduction in the number of contacts further, we analyzed the trajectory output in the 3D view. We could use the ligand RMSD graph to determine frames at which tetrazole ring flipping may have occurred. The interactions between both N5 and N6 of the tetrazole to the ARG13 residue were most prominent from Table 1, and thus a measurement was taken between these points at the start of both in-complex MD simulations, as seen in Figure 4.

Figure 4. The measurements taken at 0, 2, 8 and 10ns between N5 and N6 of the tetrazole and the ARG13 residue in both the (A) “Just neutralize” and the (B) 150mM ionic strength simulations.

By investigating the protein-ligand measurements throughout both simulations, we were able to demonstrate the tetrazole ring was less stable during the MD simulation where the system was just neutralized, compared to when the ionic strength was set to 150mM. This was evident as increased flipping of the tetrazole ring was observed within the 3D view window, thus producing more varied contacts, with N4 and N3 on the other side of the ring (Table 1, A).

Finally, to provide further evidence that the tetrazole ring was less stable at the binding site when "Just neutralize" was used, the torsional angle (C14-C15-C16-N6) highlighted in Figure 5 was considered. For this comparison, another 10ns Dynamics simulation was run with the ligand in solution. When plotted, as shown in Figure 6, it is evident once more that the charged tetrazole ring is significantly more stable when the ionic strength was set to 150mM as the torsional angle between 144-180o is present in significantly more frames than any other angle. Comparatively, when "Just neutralize" was used, the torsional angles observed were more variable, not dissimilar to when ligand 33 was in solution.

Figure 5. The torsional angle (C14-C15-C16-N6) was measured throughout the three MD simulations for comparison.

Figure 6. Bar graph which depicts the torsional angle (C14-C15-C16-N6) of the tetrazole ring of ligand 33 throughout the three MD simulations, "Just neutralize" (dark blue), 150mM ionic strength (orange) and ligand in solution (light blue).

Overall, in this study we have been able to demonstrate the importance of mimicking the assay conditions of ionic strength, especially when considering charged ligands, utilizing the many molecular Dynamics analysis tools provided by Flare.

For a visual demonstration of Flare's Molecular Dynamics environment, workflows and analysis tools, watch our webinar.

To try out the full functionality of Flare on your projects, request an evaluation today.

References

1. W. Chen et al. Accurate Calculation of Relative Binding Free Energies between Ligands with Different Net Charges. J. Chem. Theory Comput. 2018, 14, 6346-6358
2. J. Liang et al. Lead Optimization of a 4‑Aminopyridine Benzamide Scaffold To Identify Potent, Selective, and Orally Bioavailable TYK2 Inhibitors. J. Med. Chem. 2013, 56, 4521−4536
3. R. E. Farrell Jr. RNA and the Cellular Biochemistry Revisited. RNA Methodologies, 3rd edition. 2005, 1-16. Academic Press