PyAutoFEP: an automated Free Energy Perturbation workflow for GROMACS integrating enhanced sampling methods

CNPq - Conselho Nacional de Desenvolvimento Científico e Tecnológico FAPEMIG - Fundação de Amparo à Pesquisa do Estado de Minas Gerais CAPES - Coordenação de Aperfeiçoamento de Pessoal de Nível Superior Cálculos de Perturbações de Energia Livre (Free Energy Perturbation, FEP) são, na atualidade, rotineiramente aplicados na descoberta de fármacos para estimar a energia de Gibbs relativa de ligação de pequenas moléculas a alvos biomoleculares de interesse. A aplicação de técnicas de amostragem estendida pode melhorar a correlação entre predições e dados experimentais, especialmente em sistemas com mudanças conformacionais. Devido ao grande número de perturbações necessárias em uma campanha de descoberta de fármacos, a preparação manual de cálculos de FEP não é viável. Neste trabalho, apresentamos o PyAutoFEP, uma ferramenta flexível e de código aberto para automatizar a preparação de cálculos de FEP para energia de Gibbs relativa de ligação. PyAutoFEP foi escrito em Python3 e automatiza a geração de mapas de perturbação, topologias duplas, construção de sistemas para simulações de Dinâmica molecular e análise. O PyAutoFEP suporta diversos campos de força, incorpora diferentes versões de métodos de amostragem REST2 e esquemas flexíveis de janelas de λ para representar a perturbação. Para validar a ferramenta PyAutoFEP, esta foi aplicada para uma série de 14 ligantes do receptor Farnesóide X. Este sistema fez parte do desafio para a comunidade Grand Challenge 2, organizado pelo grupo Drug Design Data Resources. Uma média de 88% de sinais corretos das predições foi observado e 75% das predições alcançaram erros menores do que 1,5 kcal/mol. Resultados utilizando os campos de força Amber03/GAFF, CHARMM36m/CGenFF e OLSA-AA/M/LigParGen obtiveram r de Pearson de 0,71 ± 0,13, 0,30 ± 0,27 e 0,66 ± 0,20, respectivamente. Os resultados utilizando o Amber03/GAFF e OPSA-AA/M/LigParGen foram comparáveis às melhores submissões originais ao Grand Challenge 2. A aplicação da REST2 melhorou os resultados utilizando o CHARMM36m/CGenFF (r de Pearson = 0,43 ± 0,21), porém pouco impacto foi observado para os outros campos de força. Finalmente, utilizando um modelo estatístico, foi estimada a probabilidade de descoberta de moléculas com potência ao menos uma unidade de 1 pKi melhor do que um protótipo hipotético, quando uma ferramenta como o PyAutoFEP é aplicada para triar 10 ou 100 análogos. As probabilidades são até 7 vezes maiores se comparadas àquela obtida usando seleção aleatória, quando 100 moléculas são triadas, indicando que o PyAutoFEP pode ser útil em problemas de otimização de protótipos. O PyAutoFEP está disponível no repositório do grupo, https://github.com/lmmpf/PyAutoFEP. Free Energy Perturbation (FEP) calculations are now routinely used in drug discovery to estimate the relative free energy of binding (RFEB) of small molecules to a biomolecular target of interest. Using enhanced sampling can improve the correlation between predictions and experimental data, especially in systems with conformational changes. Due to the large number of perturbations required in drug discovery campaigns, manual setup of FEP calculations is no longer viable. Here, we introduce PyAutoFEP, a flexible and open-source tool to aid the setup of RFEB FEP. PyAutoFEP is written in Python3, and automates the generation of perturbation maps, dual-topologies, system building and molecular dynamics (MD), and analysis. PyAutoFEP supports multiple force fields, incorporates two flavors of REST2 enhanced sampling method, and allows flexible λ values along perturbation windows. To validate PyAutoFEP, it was applied to a set of 14 Farnesoid X receptor ligands, a system included in the Drug Design Data Resource Grand Challenge 2. A mean 88% correct sign prediction was achieved, and 75% of the predictions had an error below 1.5 kcal/mol. Results using Amber03/GAFF, CHARMM36m/CGenFF, and OLSA-AA/M/LigParGen had Pearson’s r values of 0.71 ± 0.13, 0.30 ± 0.27 and 0.66 ± 0.20, respectively. The Amber03/GAFF and OLSA-AA/M/LigParGen results were on par with the top Grand Challenge 2 submissions. Applying REST2 improved the results using CHARMM36m/CGenFF (Pearson’s r = 0.43 ± 0.21), but had little impact on the other force fields. Finally, we estimated the probability of finding a molecule 1 pKi better than a lead when using PyAutoFEP to screen 10 or 100 analogs. The probabilities, when comparing to random sampling, increased up 7-fold when 100 molecules were to be screened, suggesting that PyAutoFEP would likely be useful for lead optimization. PyAutoFEP is available on GitHub at https://github.com/lmmpf/PyAutoFEP.