A significant improvement in Flare™ FEP is the ability to perform perturbations with ligands which have different charge states within the same FEP graph. In this article, we present the results from an FEP calculation demonstrating a charge change between ligands in the TYK2 dataset¹ and highlight the latest analysis features which are new in Flare V9.

After careful preparation of the system, our dataset consisted of 18 ligands (Figure 1), including ligand 32 which had a negatively charged carboxylate (Figure 2). To perform an FEP calculation with a charge perturbation, a water molecule is converted into an ion to maintain the neutral net charge of the whole system. The water molecule is transformed into a sodium ion (if the perturbated ligand has –1 charge change) or chloride ion (if the perturbated ligand has +1 charge change). For the non-bonded method, we used the Particle-Mesh-Ewald (PME) method, as it is more advanced than reaction field, and handles both short-range and long-range electrostatics more accurately.

Figure 1. The Flare FEP map generated for the 18-ligand, TYK2 dataset¹, highlighting the cycle which included “ligand 32” and thus the charge change.

Figure 2. The charged ligand included in the Flare FEP experiment from the TYK2 dataset¹, ligand 32, with a carboxylate group.

Flare FEP ΔQ Results and New Analysis Features

The eighteen-ligand benchmark validation study had 54 perturbations in total. We used the following parameters to optimize the results; Custom Parameters (ligand torsions), a longer equilibration time of 1500 ps, applied Grand Canonical Nonequilibrium Candidate Monte Carlo (GCNCMC) water sampling during the equilibration stage², set the solvent ionic strength to 0.150 M, and utilized the new Flare FEP feature whereby the simulation length is variable between charged and neutral ligands, the default settings of which are shown in Figure 3. Longer simulations are needed for charged ligands, as studies have shown that charge perturbations require more extensive sampling for the free energy to converge. This feature also allows you to independently extend the simulation length for any of the two legs of your simulation (free, bound). For example, if you are studying a highly flexible ligand, with several rotatable bonds, you may need to run longer in the free state, to allow the ligand to adopt all the different conformations. On the other hand, at the bound state, this ligand may be relatively rigid, because it is constrained by interactions within the binding site, so the default options for the bound state simulation length may suffice.

Figure 3. Options available for Flare FEP calculations, specifically highlighting the separate simulation length options for charged and neutral ligands.

Using these settings, the results shown in the Activity Plot provided statistics which included an R² = 0.80 and MUE = 0.45 kcal/mol; high accuracy for a dataset with mixed charges. New to Flare V9, we can also investigate the contribution associated with each link to the overall error of the benchmark map (Figure 4), coloring each link as either red, gray or blue, that signify high, average or low error contribution respectively. This allows the user to view which links are mostly contributing to the overall error, providing the opportunity to investigate and add instances which can fine tune these perturbations. In such examples, we can change the number of lambda windows, increase simulation length, and/or remove redundant links with high hysteresis. By adding an instance to individual perturbation links, we can retain the initial result for comparison, rather than overriding and losing previous calculations.

Figure 4. The benchmark map as colored by the new Error Analysis tool (highlighted in the red box) in Flare FEP. Each link is either red, blue or gray, where red is a higher error, gray is average error and blue is a low error.

Another tool we will consider is the Torsion Plot, which allows the exploration of the torsion angles of rotatable bonds throughout the FEP simulation. We use it to identify conformational changes that occur between the star and the end molecules. As shown in Figure 5, when the link between ligand 32 and 30 is selected, we can compare various torsional rotations in the chosen ligands. For example, we observe the torsional angle C12-N2-C10-N1 (highlighted in Figure 5) as having similar torsion plots in both the free and bound states for ligand 32 (right-hand graph). However, the same torsional angle of ligand 30 has more rotational freedom in the free state (blue plot, left-hand graph). Such a feature provides clear analysis around rotatable bonds within a dataset and how varying functional groups may affect torsional rotation. It can be very insightful to compare the torsions of the lambda 0 and 1 portions of each perturbation, free and bound, to help identify if there may be different conformational space available, or to investigate if custom torsions are required.

Figure 5. New Torsion Plot analysis panel. Provides graphs of the torsional angles that are modeled during the calculation and how these compare between the selected perturbed ligands. Here, we highlight the perturbation between the nitrile group (ligand 30 – left-hand graph) and a charged carboxylate (ligand 32 – right-hand graph).

The last analysis tool that we are going to show is contact analysis. When a link is selected, we can now view the contacts that are made between the protein and each of the ligand atoms which are involved in the perturbation. For example, in Figure 6 we can see when λ=0.0, the contacts that occur most of the time throughout the simulations of all lambda windows are hydrogen bonds between VAL94 of the protein and N1 or H9 of the ligands. We can analyze each lambda window individually, looking for differences in the contacts made between the ligands, or as a group, to determine which contacts are consistent throughout the simulation.

Figure 6. The new contact analysis in Flare FEP, highlighting contacts made between the protein and two ligands of the selected FEP map link across the lambda windows.

You can learn more about the new features which are now available in Flare V9 in our recent article and webinar.

References

1. J. Liang et al. Lead Optimization of a 4Aminopyridine Benzamide Scaffold To Identify Potent, Selective, and Orally Bioavailable TYK2 Inhibitors. J. Med. Chem. 2013, 56, 4521−4536
2. O. J. Melling et al. Enhanced Grand Canonical Sampling of Occluded Water Sites Using Nonequilibrium Candidate Monte Carlo, J. Chem. Theory Comput., 2023, 19, 1050-1062