Quantitative Structure-Activity Relationship, Activity Prediction, and Molecular Dynamics of Non-nucleotide Reverse Transcriptase Inhibitors

This study used in-silico strategies to identify Enumerated Etravirine as a promising therapeutic agent for HIV. Our findings on molecular interactions and dynamics support the rational design of novel NNRTIs as possible HIV treatment alternatives.

The increasing incidence of HIV-1 drug resistance presents a challenge to the effectiveness of combination antiretroviral therapy, particularly in Southern Africa. The development of resistance to non-nucleotide reverse transcriptase inhibitors (NNRTIs) threatens the long-term success of antiretroviral therapy. In 2019, antimicrobial resistance directly accounted for an estimated 1.27 million deaths globally. This study employed an in-silico approach to investigate NNTRI drugs and their derivatives. Techniques used included density functional theory calculations, molecular docking, enumeration, quantitative structure-activity relationship (QSAR) analysis, molecular dynamics simulation (MDS), and molecular mechanics with generalized Born and surface area methods. The analysis focused on various pyrimidine derivatives and six NNRTI drugs, examining their interactions with the HIV-1 protein (PDB code 1HQU).

A QSAR model was developed to predict the biological activity of the six NNRTIs under investigation. Using 94 pyrimidine derivatives, the QSAR model achieved an R² OF 0.822 and a Q² of 0.815, indicating a high level of predictive accuracy.

MDS were conducted to assess the stability of various ligands and their newly developed alternatives, ensuring they remained bound to the protein's active site over a 200-nanosecond simulation timeframe. Etravirine exhibited root mean square deviation (RMSD) fluctuations of approximately 4.5 Å, while its enumerated derivatives showed RMSD fluctuations of 3.5 Å. Through molecular docking, MDS, and free energy computations, enumerated Etravirine demonstrated the best performance, with an activity value of 7.373 and a docking score of -10.517 kcal/mol. Furthermore, the calculated free energy of binding for enumerated Etravirine was -89.684 kcal/mol, outperforming other ligands under investigation. The significant improvement suggests that the modified Etravirine holds promising potential as a novel agent in antiretroviral therapy.

The obtained lower RMSD value, the enhanced amino acid interactions, and the highest binding free energy indicate that enumerated Etravirine could serve as a viable alternative for HIV/AIDS treatment.

Despite remarkable strides in treatment, the grave global health menace of acquired immunodeficiency syndrome (AIDS), caused by human immunodeficiency virus-1 (HIV-1), remains a persistent threat¹. Antiretroviral medications known as reverse transcriptase inhibitors, RTIs, are used to treat HIV infection. Reverse transcriptase, a viral deoxyribonucleic acid (DNA) polymerase enzyme essential for retrovirus replication, is inhibited by reverse transcriptase inhibitors (RTIs). Nucleoside reverse transcriptase inhibitors (NRTIs) and non-nucleotide reverse transcriptase inhibitors (NNRTIs) are the main RTIs².

Highly active antiretroviral therapy (HAART), a combination of various antiviral medications, has emerged as the standard therapy for HIV, effectively controlling AIDS spread and transforming this once-fatal illness into a manageable chronic condition³. NNRTIs for HIV-1 are currently a significant part of the HAART regimen⁴. Antimicrobial resistance (AMR) is among the most critical global health threats, with an estimated 1,27 million fatalities⁵. As a result, there is an urgent need to troubleshoot AMR and identify novel antimicrobial agents. According to the World Health Organization (WHO), AMR has reached alarming levels in many parts of the world.

This poses a severe risk to achieving sustainable development goals, as it undermines food security, economic growth, and health security while also contributing to social and economic inequalities⁶. The prevalence of drug-resistant viruses is on the rise, even in antiretroviral medications designed to combat HIV. According to recent statistics, nearly 30 million people across the globe were taking antiretroviral therapy as of the end of 2022⁷.

WHO conducted 30 surveys and discovered that in 21 of those surveys, more than 10% of individuals who began their first-line antiretroviral therapy had resistance to Nevirapine or Efavirenz⁸. Moreover, individuals with prior exposure to antiretroviral drugs are up to three times more likely to have resistance to NNRTIs than those without exposure⁹. Studies have revealed that a sizeable number of infants under 18 months of age who were newly diagnosed with HIV exhibited high rates of drug-resistant strains^10,11. Shockingly, almost half of them had a - (NNRTI)-resistant strain even before starting the treatment. These findings underscore the need to expedite studies to design innovative HIV therapy.

The development of drug-resistant strains and undesirable side effects from prolonged use have unavoidably posed challenges to the clinical use of NNRTIs¹². Initiating combination antiretroviral therapy (cART) extends the life expectancy of people living with HIV despite the potential for adverse side effects. These side effects encompass the risk of developing non-communicable diseases, including lipodystrophy, hyperlipidemia, reduced bone mineral density, elevated blood glucose level leading to type 2 diabetes mellitus, hypertension, an increased risk of stroke, and issues related to obesity¹³. Early diagnosis and prompt access to suitable medical care in the initial phase of HIV infection offer considerable advantages both from a clinical and public perspective. Timely initiation of ART and prophylaxis against opportunistic infections leads to a notable reduction in HIV-related illness and mortality.

The use of ART can also contribute to a decline in the potential for onward transmission of HIV by reducing the level of circulating HIV ribonucleic acid. Furthermore, addressing the treatment of other sexually transmitted diseases and co-infections can likewise diminish the probability of further HIV infection¹⁴. NNRTIs are amongst the optional treatment classes of inhibitors that interact with RT by binding to the allosteric region or site of HIV. This type of inhibition is commonly known as an uncompetitive inhibitor because the NNRTI does not bind at the active site of the substrate but rather on the outside. The effect alters the conformation of the substrate binding site, hindering the standard substrate from binding and leading to premature chain termination. Due to their reduced toxicity compared to NRTIs, simple structure, improved bioavailability over protease inhibitors, and excellent selectivity, NNRTIs have become the most attractive HIV inhibitors¹⁵. Therefore, the synthesis and design of novel NNRTIs are crucial from a pharmacokinetic perspective^16,17.

NNRTIs are vital in HIV/AIDS treatment due to their strong efficacy and low toxicity. However, early NNRTIs like Nevirapine, Delavirdine, and Efavirenz encounter resistance from viral mutations in the NNRTI binding site¹⁸. These drugs, part of the diarylpyrimidine (DAPY) family, exhibit potent activity against various NNRTI strains, including those resistant to early NNRTIs¹⁹, likely owing to their molecular flexibility and higher resistance barrier against HIV-1^20,21. Despite their success, the high mutation rate in HIV-1 RT and the absence of intrinsic proofreading activity have led to new resistance profiles in patients using Etravirine and Rilpivirine^22,23. These viral mutations differ from each other, as do those associated with early NNRTI drugs²⁴. Over 50 structurally diverse classes of compounds have been identified as NNRTIs. Notably, six NNRTIs have obtained approval for the treatment of HIV-1. These approved agents encompass Nevirapine (NVP), Delavirdine (DLV), Efavirenz (EFV), Etravirine (ETR), Rilpivirine (RPV), and Doravirine (DOR). Figure 1 shows the chemical structures of these six approved NNRTI drugs²⁵.

In a recent study²⁶, DFT calculations and molecular docking suggest that lovastatin and simvastatin show potential as anti-coronavirus agents. Virtual screening identified five new FDA-approved candidate molecules similar to the efavirenz scaffold, with superior binding affinity in the active pocket of the COVID-19 main protease.

Soltan et al.²⁷ carried out a similar study on the identification of novel molecules for improved binding capacity to HIV RT by employing a fragment-based strategy with FDA-approved drugs to design chemical derivatives. They specifically leveraged the structures of delavirdine, efavirenz, etravirine, and rilpivirine as foundational scaffolds. The drug-likeness of these derivatives was assessed through Swiss-ADME, followed by docking them into related crystal structures. The study concluded with the selection of compounds exhibiting superior binding affinities compared to their parent scaffolds, particularly noting a more pronounced improvement in derivatives designed from the second-generation NNRTIs, etravirine and rilpivirine. For instance, derivatives RPV01 and RPV15 showed significant enhancements in docked energy values over rilpivirine, indicating potential utility in targeting both wild-type and mutant forms of HIV RT.

Murugesan and coworkers²⁸ conducted research propose various medicinal chemistry approaches to enhance the efficacy and minimize resistance in HIV treatment. The study employed molecular hybridization, bioisosteric replacement and high-throughput screening. In their study, they managed to identify novel NNRTI scaffolds with high potency against both wild-type and drug-resistance strains of HIV. They developed DAPY derivatives that exhibit excellent selectivity and low toxicity, with some compounds showing effective inhibition at nanomolar concentrations.

Lately, using computational tools alongside in-silico research has gained popularity due to its practical analysis of a drug's quantum chemical characteristics²⁹. In this study, computer-aided drug design, density functional theory, qualitative structure-activity relationship, and molecular dynamics were applied to discover potent NNRTIs.

1. Computational details—protein preparation

1. Click on the Window icon on the monitor screen and click All Apps. Scroll down to navigate to the Schrodinger folder, open the folder, and click on the Maestro icon shown in Figure 2B, and select open as shown in Figure 2B to launch the software.
2. Retrieve the Protein structure of choice by navigating the File tab button in the software. From the short menu that pops up, select Get PDB, as shown in Figure 3A, and enter the PDB code of choice in the text box as shown in Figure 3B. Click on the download button, and the selected PDB file will be displayed in the project window.
3. Alternatively, download the protein of choice into the local computer from the Protein Data Bank by entering the Protein database identity code (PDB ID) in the search box and click download to download the PDB file into the local computer. Navigate to the File tab and select the Import Structures option. From the Import interface, locate the downloaded PDB file as shown in Figure 3C and select then Import button as indicated in Figure 3D.
   NOTE: The protein structure will be open as a 3D structure in a separate window, as shown in Figure 4.
4. Navigate to the top right of the software, select the task option, and type protein preparation in the search bar. Click Protein Preparation Workflow in the display on the right side, as shown in Figure 5A.
5. From the protein preparation Workflow window that pops up, write the Job name as the filename to be saved and click the green Run button in the bottom right corner, as shown in Figure 5B.
6. Monitor the job as it runs by clicking the Jobs button on the top right corner, as shown in Figure 5C.
7. Select and right-click on the prepared protein and select split ligand, as shown in Figure 5D. Choose to split into ligands, water, and others. This is to have the protein components as independent entries in the workspace navigator

2. Ligand preparation

1. Download the desired chemical compounds from the PubChem database³⁰ by typing the name of the compound of choice in the search bar. Go through the structures and select 3-dimensional (3D) structures. Click Download on the top right window to download the structure coordinates as a structured data file (SDF). If a 3D structure is unavailable, download the 2D structure and use other tools to generate a 3D structure from a 2D structure.
2. Click the File tab in Schrodinger and select Import Structures, as shown in Figure 6A. Navigate to the file location where the structures are saved in SDF file format to load the compounds to be prepared.
3. Select Task in the top right corner of the Schrodinger software. Type LigPrep in the search bar and select LigPrep in the right window on the left, as shown in Figure 6B.
4. Select Use Structures from to select files from the Workspace or Project table. Select the preferred options from the LigPrep window, save the file in the local computer, and click run to submit the job for ligand preparation as shown in Figure 6C.

3. Geometry and Optimization of ligands

1. Open the software³¹ for geometry optimization of the downloaded structures. Navigate to the File tab (Figure 7A) and select Open to choose the downloaded SDF file from the PubChem database.
   NOTE: The file will load in the main window. Click the other purple window to visualize the same chemical compounds. Every implemented setting will show the result in the purple window.
2. Navigate to the Calculate tab and select the Gaussian Calculation Setup Tab. Navigate to the Job type tab shown in Figure 7B and choose Optimization or Opt+Freq. 
   NOTE: Depending on the (size) number of atoms or compounds, optimize first, followed by frequency. If the compound is smaller, perform Opt+Freq. The bigger the molecule, the more time it takes to perform both Opt+Freq compared to Optimization followed by Frequency.
3. Navigate to the Method tab and select the quantum chemistry methods and then the Kohn–Sham global-hybrid exchange-correlation density functional of choice, basis set, charge, and spin of choice from the dropdown arrows in each section.
4. Navigate to the Title and specify a name for the compound under investigation.
5. Navigate to the Link 0 tab and specify the Memory limit and Shared processors of choice. Untick the Full path boxes as shown in Figure 7C.
6. Click the Edit button at the bottom to save the Gaussian input file, as shown in Figure 3C. Save the file in the preferred location with a file name of choice as a Gaussian job file (gjf).
   NOTE: After saving the Gaussian input file, a pop-up window will appear showing the file's contents in Notepad. Edit the file's contents, including options like changing the charge and functional that were unavailable in the settings options. The prepared Gaussian job file will be used as the input file for submission to run the optimization and frequency calculations via a local computer.
   Subsequently, geometric optimization of these structures was conducted using the MN15-L functional³² and 6-31++G(d,p) basis set³³. However, should issues arise in the running of the optimization and frequency process, one can submit the work with the aid of an HPC.

4. Receptor grid generation

1. Navigate to Tasks and select receptor grid generation³⁴. The receptor grid generation interface shown in Figure 8A identifies the protein's active site where the core-crystal ligand is bound. Click on Pick to identify the ligand and check if a co-crystallized ligand is present from the pop-up top notification shown in Figure 8B.
2. Select the ligands and/or residues from the workspace, and the selected compounds of interest will show a blue highlight in the workspace.
3. Select the settings tab in the receptor grid generation panel to set the grid box size. The default size of the grid box is 10 Å x 10Å x 10Å. Modify the box dimensions in the Grid Box tab by inputting the desired values directly or manually resizing the box within the workspace.
   NOTE: Alternatively, set additional parameters from the Advanced tab as shown in Figure 8C if necessary. Review all settings to ensure accuracy.
4. Click Run to begin grid generation. Once complete, save the generated grid files for subsequent docking simulations. A notification message will show when the job is completed, as shown in Figure 8D.

5. Molecular docking

1. Load the protein and the prepared ligands by navigating to Tasks, Docking³⁵, and Ligand Docking (glid docking), as shown in Figure 9A.
2. Load the grid file from step 4.1 above and select the ligands from the workspace using the Use ligand from the option in Figure 9B.
3. Choose the preferred docking precision method from the settings tab shown in Figure 9 C (default docking precision is Standard Precision (SP)).
4. Set the force field to OPLS4 for accurate modeling of molecular interactions. Set constraints (e.g., hydrogen bonds) in the Constraints tab.
5. Review all settings and save the docking job or file. Click RUN to start the docking process.
6. A notification indicates the job has been completed, and the Project Table is opened for docking results. Examine the binding poses, scores, and interactions.

6. 2D-QSAR model generation

1. Go to the KNIME community hub webpage and search for AutoQsar in the search bar. Select the second AutoQsar entry that appears shown in Figure 10A.
2. Click on the download workflow icon and select download workflow in the pop-up menu shown in Figure 10B. Save the workflow to the local computer.
3. Ensure that KNIME³⁶ is installed, then open it from the local computer. Compile a spreadsheet with the canonical smiles, name of the structure, activity/IC₅₀ values (in micromolar), and -log(activity/IC) or any chemical descriptor of choice. Save the file in a Comma-separated value (CSV) format on the local computer.
4. Navigate to File and import the AutoQSAR workflow³⁷. From the pop-up window, click Import KNIME Workflow to select the AutoQSAR workflow downloaded, as shown in Figure 10C and select Open | Next | Finish.
5. The AutoQsar workflow will show on the KNIME Explorer window on the top left. Double-click on the workflow to bring the workflow into the main window.
6. Double-click on the molecular reader (to MAE) icon, and from the Settings tab, click on Add file(s) to add the worksheet with the data input data or parameters.
7. Select the Memory Policy tab and click on Write Tables to disc | OK.
8. Right-click on Extract Properties button and select Configure, as shown in Figure 10D.
9. Select and filter molecular or chemical properties of choice from the exclude window into the Include window. Click Apply.
10. Right-click on AutoQSAR Build Model and click on Configure. Input the QSAR model name into the pop-up dialog.
11. Click on the setting button and select the property to be fit from those included. Choose the random training set value (e.g., 85%:15%) and several models to keep. Click Apply.
12. Right-click on the bottom Molecule Reader and select Configure. Load the test ligands from a prepared CSV Excel file containing the ligands to be tested.
13. Right-click on Extract Properties and set the settings to your choice.
14. Right-click on the top Molecule Reader and select execute. 
    NOTE: An orange dot indicates that the process has commenced, and a green light will show the successful execution of the process. The next step will carry on until all the nodes have been successfully executed.
15. .Execute the second Molecule Reader for Test ligands.
16. Save the zip folder with the results from the test and training set after the successful execution of all the nodes.
    NOTE: Select the best-performing model based on cross-validation results such as SD –standard deviation, R2 – training set correlation between actual and predicted activity values and Q2 – test set actual and predicted activity correlation.
17. Evaluate the workflow performance by assessing the model's performance using a scorer or statistical nodes.
18. To interpret the results, identify the main or key descriptors using the feature importance. Visualize the results as a scatter plot or bar chart to determine the performance and descriptor relationship. Save the workflow and export the result data, predictions, and visualizations as a scatter plot and/or CSV files.

7. Enumeration

1. Navigate to Task and search for Ligand Designer, as shown in Figure 11A.
2. Select a pair of the docked protein and the ligand complex from the workspace navigator. Click Analyze Workspace from the Ligand designer window. To generate and evaluate new ligands, select Isostere scanning from the list of workflows that appear, as shown in Figure 11B implies the growing method that extends the ligand by adding fragments to existing parts of the molecule.
3. After setting the enumeration option, click enumerate from the Isostere Scanning notification window that appears. Repeat step 6.2 for the same protein and a different ligand until all ligands have run the same process. Review the enumerated ligands from the Projects Table.
4. Save the structures by exporting the designed ligands.
   NOTE: The enumeration method generated a set of docking scores when the simulation was completed.
5. Individually select the enumerated compounds with the more negative docking scores and redock them using the docking procedure outlined in section 4.

8. HOMO and LUMO

1. Open the software and upload the optimized ligand in Gaussian job file format.
2. Navigate to tools and select MO editor in red and green dots, as shown in Figure 12A.
   NOTE: The enumeration method generated a set of docking scores when the simulation was completed.
3. In the MO editor panel under method, load the FChk file by clicking on load Mos from the existing Chk or FChk file, as shown in Figure 12B.
4. Click on the visualize tab | update. Wait for ~10 s for the current surface to appear, as shown in Figure 12D.
5. Click on one of the two tick boxes next to the highlighted yellow figures to select either HOMO or LUMO³⁸ surface to display in the window next to it.
6. Right-click on the purple background and select File. Click on save the image file to save the image of the current surface, as shown in Figure 12E.
7. Right-click on the purple background and select View. Choose the display format to change the background in the General tab, surface transparency in the Surface tab, font size and colour in the Text tab, image quality and the surface's preferred layout under the Molecule tab, as shown in Figure 12F.
8. Alternatively, generate a cube file and upload the FChk file in GaussView. Navigate to the Results tab and choose Surfaces/Contours. Click on cube actions and load the cube. Click on Surface Actions and choose a new surface. Repeat step 8.7 to edit the surface.
   NOTE: The smaller the energy gap between the HOMO and LUMO³⁹ (the difference between LUMO and HOMO), the more reactive a molecule is expected to be. Calculate the energy gap (E_gap) using Eq 1 for each molecule.
   

9. Molecular dynamics simulation—System preparation and minimization

1. Ensure that the Schrodinger Suite has been installed in the local computer and load the protein-ligand complex structures in the workspace⁴⁰.
2. Click on the Task button and choose Desmond System Builder⁴¹. In the system builder panel, select the Solvation tab and choose the Predefined solvent model as shown in Figure 13A, Box shape as shown in Figure 13B, and box size calculation method as shown in Figure 13C⁴², that suit the protein-ligand complex.
3. Select the ions tab and click recalculate to neutralize the system by adding counterions and setting the desired salt concentration.
4. Choose the OPLS^43,44 of choice as the force field.
5. Name the job appropriately and save the job file on the local computer. Click run to submit the job for preparation.
   NOTE: Ensure that the job name is saved. Use a detailed name.
6. Equilibration and production
7. View the project in the workspace After the system preparation. Select the protein-ligand complex from the Workspace Navigator, navigate the task, and choose molecular dynamics (Desmond).
8. Load the ligand-protein complex from the workspace in the molecular dynamics panel. Select the desired simulation timeline from the simulation tab. Select NPT as the ensemble class⁴⁵.
9. Name the job appropriately and write out the job. Click on Close to exit the molecular dynamics window.
10. Submit the written-out job for molecular dynamics preparation via a local terminal. Upon completion, open the completed job and continue the simulation time from the initial set simulation timeline until the desired simulation time⁴⁵, for example, 100 ns, 200 ns.
    NOTE: Repeat step 9 for the other protein-ligand complexes.
    1. Open the completed job in Schrodinger Maestro software and merge the different simulation timeframes if separate jobs were run. Navigate to the TASK button and select the simulation interaction diagram. Load the merged files or the single file to visualize the trajectory.
11. Generate a report that includes the visualizations and other detailed output results.

10. Molecular mechanics with generalized Born and surface area (MM-GBSA)

1. Open the trajectory file and play the trajectory. Visualize where the protein-ligand complex is equilibrated and note the number of frames. Submit the job via the terminal for processing.
2. Repeat section 9 (prepare the protein molecular dynamics simulation) for the other protein-ligand complex.
3. Catenate/view the output file contents to analyze the results generated. Read the ΔG Average to obtain the binding free energy of the protein-ligand complex.
4. Download the CSV file to visualize the contributions of different intramolecular molecules.
5. To compute the free binding energy of the complex, take note of the various thermodynamics and desolvation parameters, including the binding energy (ΔG_bind), Columbic solvation model (ΔG_bind Coulumb), nonpolar solvation term (ΔG_bind Lipo), hydrogen-bonding correction (ΔG_bind Hbond), covalent binding (ΔG_bind Covalent), π-π packing correction (ΔG_bind Packing), Generalized Born electrostatic solvation energy (ΔG_bind sol GB), and van der Waals interaction (ΔG_bind vdW).
   1. Derive the final ΔGbind by averaging the ΔG_bind values determined for each snapshot within the MD simulation, as shown in Eq 2.

Receptor grid generation and molecular docking

Maestro's receptor grid generation tool was used to properly characterize the binding site for the subsequent docking. The co-crystallized ligand was used to define the grid. SP precision Glide setting for molecular docking was used. The LigPrep tool in Schrödinger Maestro was used to prepare the ligands for docking using the OPLS4 force field. For Molecular dynamics simulations, the OPLS4 force field was used in Desmond. Force fields are fundamental to classical molecular simulations, and their accuracy is critical for the quality of protein-ligand binding simulations in drug discovery. For OPLS4, charge and parameter assignment were completed using Schrödinger Maestro. The application of OPLS4 parameters resulted in significant improvements in both energetic and geometric comparisons compared to the default parameters of OPLS2005.

This entailed docking specific ligands known to have an affinity for the HIV-1 target protein, allowing for a thorough analysis of the interactions between ligand molecules and receptor residues. Table 1 represents a summary of compound classification for QSAR modeling.

Upon docking HBY561, the docking protocol was assessed by comparing the redocked ligand with the one found in the active site of the crystalline protein 1HQU. The docked and crystallized HBY561 structures are provided in Supplemental File 1 (Supplemental Figure S1 and Supplemental Figure S2). To evaluate the similarity between docked poses and reference structures, a root-mean-square deviation (RMSD) value of less than 2.0 Å is widely regarded as a criterion for reliable docking results. This threshold indicates that the predicted structure aligns closely with the experimental data. In this study, the ligands demonstrated an RMSD of 1.27 Å when comparing the reference structure to the docked pose, as illustrated in red in Supplemental Figure S3. This showed that the docking protocol was sufficient for this work, and as a result, all ligands were docked using the same settings. Examining the docking scores presented in Table 2, Efavirenz and Etravirine have shown the most favorable scores at -10.432 eV and -9.647 eV relative to the docking score of the co-crystallized ligand HBY561 (-9.242 eV). Supplemental FigureS4 from Supplemental File 1 shows the ligand interaction diagrams between the original HIV-1 protein and Crystal ligand, Etravirine, and Efavirenz.

Hydrogen bonding, π-π stacking, and hydrophobic interactions are the primary forces contributing to binding. A specific instance of hydrogen bonding emerged between HBY561 and the designated protein 1HQU, explicitly involving the amino acid residue LYS101. This bonding pattern mirrored observations made with the ligands Efavirenz and Etravirine, as shown in Figure 14.

Moreover, hydrophobic interactions were essential for binding in several protein locations involving HBY561, Etravirine, and Efavirenz, in addition to the intermolecular forces, hydrogen bonding and π-π stacking. π-π stacking was observed between TYR318 and the aromatic ring in Efavirenz. Both hydrogen bonding and π-π stacking were essential in maintaining the binding connections between the ligands and the protein. Ligand Interaction diagrams of HBY_561, Nevirapine, Doravirine, Efavirenz, and Etravirine are presented in Supplemental File 1-SupplementalFigureS5.

These intermolecular forces influence protein-ligand interactions and are crucial for developing medicinal drugs that reduce AMR in HIV-1. Their role in enhancing binding affinity, specificity, and mechanism of action aids in designing drugs that can effectively target and inhibit the virus, thereby addressing the growing concern of antimicrobial resistance in the context of HIV-1 treatment.

2D-QSAR data set preparation

The training and testing phases involved 94 compounds. These compounds were categorized into four classes, each representing ligands associated with the specific protein tested. The training process utilized the KNIME AutoQSAR workflow with activity, HOMO, and LUMO were selected as the three descriptors for this investigation.

2D-QSAR generation

All 94 molecules had their activity values determined through experimental data (Supplemental Table S1). The results generated from the QSAR model can be found in Table 3. Class 1 compounds in Supplemental Table S2 were used for QSAR modeling, showing the added Ar-groups and their respective activity values. The results from the four classes indicate that the highest score of 0.8223, R² of 0.815, and Q ² of 0.8182 were achieved, corresponding to Class 1. This aligns with the previous criteria of aiming for R² close to 1 and Q² greater than 0.7⁴⁶. Consequently, we selected class 1 for the training of our QSAR model. While models for classes 3 and 4 demonstrated an excellent correlation value R² of 0.8172 and 0.6673, respectively, they did not match the performance of Class 1.

Cross-correlation was conducted to further validate the stability of the proposed QSAR model according to the criterion that the difference between the Q² score 0.8223 and R² should be less than or equal to 0.3⁴⁷. Our proposed model for class 1 has a difference of 0.0038. The scatter plot depicting the observed activity versus the predicted activity is provided in Figure 15.

HOMO-LUMO energy gap

The determination of the energy gap between the lowest unoccupied molecular orbital and the highest occupied molecular orbital, commonly known as the HOMO-LUMO energy gap, plays a crucial role in characterizing a molecule's chemical reactivity and kinetic stability in the context of the six NNRTI compounds. The frontier molecular orbitals play a pivotal role in facilitating charge transfer interactions with the binding site of the HIV protein. The confirmation of the energy minimum was ensured by examining vibrational frequencies and confirming the absence of negative or imaginary frequencies; subsequently, the HOMO and LUMO values were obtained for each energy minimum. A higher HOMO value signifies a molecule's proficiency as an electron donor, while a lower value suggests it acts as a weak electron acceptor. Supplemental Figure S6 represents the HOMO_LUMO energy gaps for the six optimized NNRTIs. Furthermore, a reduced energy gap between the HOMO and LUMO energy levels strongly influences intermolecular charge transfer interactions that occur between the studied molecules due to the strong electron-accepting ability and bioactivity of the molecules⁴⁸.

The trend in the energy gap values, as presented in Table 4, follows a descending order: Efavirenz > Etravirine > HBY-561 > Nevirapine > Delavirdine > Doravirine > Rilpivirine. The substantial energy gap observed for Efavirenz and Etravirine signifies that the analyses of docking scores reveal a correlation between bioactivity and the HOMO-LUMO gap. Notably, the antiviral potential increases with larger HOMO-LUMO gap values. It indicates not only the stability of the compounds but also their potential to form stable interactions with the receptor. The HOMO-LUMO gap plays a significant role in understanding the bioactivity of molecules, particularly in the context of HIV-1 drug design.

Enumeration

Various tools have been created for the enumeration of virtual libraries. Tools used for enumeration include Schrödinger. It relies on the core hopping method, where libraries are created by substituting one or several attachments on a core structure with fragments from reagent compounds⁴⁹. The enumeration tool in Maestro v13.1 was used to add custom side groups or atoms to each of the six NNRTIs. The new compounds were also used to predict activities. There was an improvement in the expected activity values of the enumerated molecules compared to the initially optimized NNRTI molecules, as can be seen in Table 5.

The improvement shown in the activity values of the Enumerated ligands was due to the enumeration process performed as the added custom R-group influenced the forces of interaction of the novel proposed compound and the original protein. The enumerated compounds were prepared for quantum mechanics by optimizing these molecules and calculating their vibrational frequencies. Their energy gaps were calculated and compared to the energy gaps of the NNRTIs. The overall observation, as shown in Table 6, indicates that the enumerated compounds are more stable than their optimized counterparts.

In Table 6, it was observed that the energy gap of the enumerated compounds compared to the optimized compounds showed a similar trend. This indicates that the chemical properties of the enumerated molecules remained unchanged, regardless of any conformational rotation. Thus, they were able to maintain important intermolecular forces with the amino acid residues of the protein.

Before performing the enumeration process for the NNRTIs as outlined in protocol section 6, the observed output results had their initial docking scores from the enumeration process, represented in Table 7 under the enumerated docking score column. As explained in protocol section 4, the molecular docking process was carried out to validate the proposed docking score for the enumerated compounds. It could be observed that after redocking the enumerated compounds, the new docking scores improved, as shown in the 'redocked enumerated ligands' column. The redocked scores for the enumerated compounds were compared to the docking scores of the original NNRTIs represented under the 'original docking score' column in Table 7. It could be observed that for the enumerated compounds, the docking scores for HBY_561, Etravirine, Efavirenz, and Doravirine were better than those of their equivalent optimized compounds. However, Delavirdine possesses the same docking score as enumerated and optimized Delavirdine.

Molecular dynamics

Three hydrogen bonds are formed by HBY 561 with LYS101, the highly electronegative N atoms, and the OH. The crystal ligand, or third molecule, establishes two hydrogen bonds simultaneously with sulfur and hydrogen and a third hydrogen link with GLU138. Importantly, π-π stacking and hydrophobic interactions significantly contribute to the intermolecular forces involved. Furthermore, the presence of additional hydrogen bonds is crucial for molecular dynamics, particularly reflected in the resulting root mean square deviation (RMSD) graphs. In all, three hydrogen bonds are observed to be formed by Efavirenz with the very electronegative nitrogen atom, the oxygen atom on the other non-aromatic ring, and the benzene ring and TYR318 between the ligand highly electronegative N at the central ring. A second hydrogen bond between N from the ring and the OH and LY101 is seen. Etravirine shows three hydrogen bonds with LYS101. Doravirine forms a hydrogen bond with GLU138. Nevirapine shows two hydrogen bonds with LY101.

In this instance, the amino acid residue interaction shared by all the ligands is LYS101. Even though their structures differ, they all interact with the same amino acid residue. MD simulations were carried out with the parameters listed in protocol section 2.8 to ascertain how well or poorly each ligand (NNRTI and the listed NNRTI) binds to the active site of 1HQU. The ligand interaction depicted in Figure 16 indicates robust hydrogen bonds between the amino acid LY101 of the protein and the molecules HBY 561, Nevirapine, Efavirenz, and Etravirine. As can be seen in the comparison in Table 7, these strong interactions account for each compound's high docking scores.

To assess the binding efficacy of each ligand, including NNRTI and enumerated NNRTI at the active site of 1HQU, molecular dynamics simulations (MDS) were conducted. Specifically, the four enumerated compounds that demonstrated better docking scores than the original optimized combos were selected. These chosen compounds were subjected to MDS as a validation method to examine and observe the reaction of each molecule with the HIV-1 protein over a specific period, considering the interatomic interactions in the presence of a ligand.

Before starting MDS for the recently enumerated glycans, it was essential to confirm the suitability of the simulation protocol for our system. To achieve this, the initial step involved running MDS of the free protein 1HQU. No ligand was present within the protein's active site (Figure 17A and Supplemental Figure S7). Up to ~60 ns, there are considerable fluctuations in the protein structure, producing Cα RMSD shifts of up to 4.5 Å; after that, the protein seems to stabilize with an RMSD fluctuation of ~3.5 Å up to 200 ns. This stabilization convinced us that the MD protocol would be appropriate for our protein-ligand complexes, which will be analyzed in the following section.

The 200 ns MDS of Etravirine and enumerated Etravirine (Figure 17B,C) showed root mean square deviation (RMSD) fluctuations of Etravirine close to 5.0 Å and equilibration of 4.5 Å. The enumerated Etravirine showed RMSD fluctuations of 4.5 Å and equilibration of 3.5 Å. This stabilization indicates that enumerated Etravirine can be a potential NNRTI ligand for treating HIV/AIDS. A trajectory with an RMSD less than 5 Å signifies a robust binding effect between the active site's protein and ligand. This observation was held for all the previously mentioned compounds, except for Nevirapine and Doravirine, among the enumerated compounds (Supplemental File 1: Supplemental Figure S8, Supplemental Figure S9, Supplemental Figure S10, and Supplemental Figure S11).

Figure 18 and Figure 19 further analyze the binding between Etravirine, enumerated Etravirine, and the protein. The analyzed data include a histogram of the Interaction, ligand-protein, and protein-ligand contacts. The interaction contacts histogram for each respective ligand bound to the protein is directly related to the corresponding interaction forces between the residues of the protein's amino acid and the ligand. The high abundance of LYS101 is very distinct for Etravirine and enumerated Etravirine, and a visible thick orange band was observed. A faintly discernible light orange band positioned at the lower extremity of the enumerated Etravirine chart was observed in correlation with TYR181. This correspondence indicates the existence of two intermolecular forces of attraction between GLU138. This positive observation for the ligands and their enumerated forms were compared to analyze a better ligand as a potential NNRTI compound. Based on the results provided, the enumerated Etravirine possesses the potential to be used for the treatment of HIV/AIDS. 

Molecular mechanics with generalized Born and surface area (MM-GBSA) calculations

In this study, the primary source of energetic input towards the binding free energy, ΔG_bind, was the contribution of the van der Waals, ΔG_VdW, interaction. Enumerated Etravirine has a higher ΔG_VdW of -66.146 kcal/mol than its known NNRTI counterpart with a ΔG_VdW of -64.669 kcal/mol. The higher ΔG_Hbond value of -2.541 kcal/mol enumerated Etravirine indicates the significant contribution of the hydrogen forces of attraction between the ligand and the protein. The contributions of the ΔG_Coulomb and ΔG_Covalent for enumerated etravirine (-17.976 and 2.807 kcal/mol) were far greater than those of the known NNRTI counterpart, which had -11.196 and 2.491, respectively.

The observed results during the binding of Etravirine and the enumerated Etravirine with the 1HQU protein are somewhat consistent. Enumerated Etravirine was found to be better than its counterpart, Etravirine, due to two additional hydrogen bonds. The study also revealed that enumerated Etravirine was preferred for binding to the identified pocket. The more negative ΔG_bind (-89.684 kcal/mol) for enumerated Etravirine compared to that of Etravirine (-80.551 kcal/mol) shows that enumerated Etravirine is a good inhibitor of HIV-1 RT.

Figure 1: Chemical structures of six approved United States Food and Drug Administration non-nucleotide reverse transcriptase inhibitors for human immunodeficiency virus-1. NVP = Nevirapine; DLV = Delavirdine; EFV = Efavirenz; ETV = Etravirine; RPV = Rilpivirine; DOR = Doravirine. Please click here to view a larger version of this figure.

Figure 2: Opening the Maestro Schrodinger application on Windows in the local computer. (A) Navigating to the Maestro Schrodinger Application in the local computer. (B) How to open and run the Maestro Schrodinger Application in the local computer. Please click here to view a larger version of this figure.

Figure 3: Importing a PDB file structure from the local computer to the project window in Schrödinger. (A) Import structure function in Maestro Schrodinger. (B) PDB ID text box. (C) Downloaded PDB file in the local computer. (D) Import button to allow the inserted PDB ID file to be imported. Please click here to view a larger version of this figure.

Figure 4: Structure of PDB file imported into the Schrodinger project window. Please click here to view a larger version of this figure.

Figure 5: Protein preparation workflow. (A) Protein Preparation Workflow search interface. (B) Saving the job file name and initiating the protein preparation process. (C) Monitoring window for running jobs. (D) Breaking down the ligand into its individual components. Please click here to view a larger version of this figure.

Figure 6: Ligand preparation workflow. (A) Importing structures from the local computer to the Protein preparation Schrodinger project window. (B) Searching for Ligand preparation process. (C) Ligand preparation workflow window. Please click here to view a larger version of this figure.

Figure 7: Geometry and optimization workflow. (A) GaussView menu window for geometry optimization. (B) Available Job Types in the Calculate Tab of GaussView. (C) Options available under the Link0 tab in GaussView. Please click here to view a larger version of this figure.

Figure 8: Glide gride generation and molecular docking workflow. (A) Receptor grid generation workflow interface. (B) Pop-up notification for selecting an atom within the ligand. (C) Receptor advanced settings. (D) Job completion notification. Please click here to view a larger version of this figure.

Figure 9: Glide ligand docking. (A) Interface for glide ligand docking Search. (B) Ligand Docking Interface. (C) Precision settings for Glide Docking. Please click here to view a larger version of this figure.

Figure 10: KNIME QSAR preparation workflow. (A) Searching for AutoQSAR node on KNIME community hub webpage. (B) Download button to obtain the AutoQSAR KNIME Node. (C) Importing the downloaded AutoQSAR KNIME Workflow. (D) Configuration settings for ligands when building a QSAR model. Please click here to view a larger version of this figure.

Figure 11: Enumeration of ligands using Ligand Designer in Maestro Schrodinger. (A) Searching for Ligand Designer options in Schrodinger. (B) List of workflows for conducting the enumeration process. Please click here to view a larger version of this figure.

Figure 12: HOMO-LUMO generation workflow. (A) To access the molecular editor options in the Tools tab in GaussView. (B) Loading an existing Chk or FChk file to generate frontier molecular orbitals. (C) Illustration of the HOMO and LUMO frontier orbitals in GaussView. (D) Visualize window to display the HOMO and LUMO frontier orbitals. (E) Saving the frontier orbitals of the HOMO and LUMO. (F) Display Format Interface in GaussView. Please click here to view a larger version of this figure.

Figure 13: Molecular dynamics preparation, setup, preparation, and execution workflow. (A) Desmond System Builder: Solvation options to determine rigid water models. (B) Desmond System Builder Boundary options to determine Box shape. (C) Desmond System Builder: Method for Calculating Box Size. Please click here to view a larger version of this figure.

Figure 14: Ligand interaction diagrams between 1HQU and 3 top docked ligands. Forces of interaction between protein (1HQU) and (A) the crystal ligand (HBY561), (B) Efavirenz, and (C) Etravirine. These forces influence protein-ligand interactions and are crucial for developing medicinal drugs that reduce antimicrobial resistance in HIV-1. Their role in enhancing binding affinity, specificity, and mechanism of action aids in designing drugs that can effectively target and inhibit the virus, thereby addressing the growing concern of antimicrobial resistance in the context of HIV-1 treatment. Please click here to view a larger version of this figure.

Figure 15: A scatter plot shows the observed activity versus the predicted activity for class 1 of the QSAR model. The graph represents the fitting between class 1 as the training set and the NNRTI compounds as the test set to give a predictive activity value. Abbreviations: NNRTI = non-nucleotide reverse transcriptase inhibitors; QSAR = quantitative structure-activity relationship. Please click here to view a larger version of this figure.

Figure 16: Ligand interaction diagrams. Forces of interactions between the protein and (A) enumerated crystal ligand HBY_561, (B) enumerated Nevirapine, (C) enumerated Doravirine, (D) enumerated Efavirenz, and (E) enumerated Etravirine. Please click here to view a larger version of this figure.

Figure 17: Molecular dynamics simulation interaction diagram of the free protein Etravirine and Enumerated Etravirine. (A) Molecular dynamics interaction diagram of the free protein. (B) Molecular dynamics interaction diagram of Etravirine. (C) Molecular dynamics interaction diagram of Enumerated Etravirine. Please click here to view a larger version of this figure.

Figure 18: Interaction contacts histogram between Etravirine and the protein. (A) Protein-Ligand Contacts timeline for Etravirine. (B) The timeline of the protein-ligand interactions over time, including H-bonds, hydrophobic, ionic, and water bridge contacts. (C) A schematic showing detailed interactions between ligand atoms and protein residues. Only interactions occurring more than 30% of the simulation time (0.00 to 200 ns) are displayed. Please click here to view a larger version of this figure.

Figure 19: Interaction contacts histogram between enumerated Etravirine and the protein. (A) Protein-Ligand Contacts timeline for enumerated Etravirine. (B) The timeline of the protein-ligand interactions over time, including H-bonds, hydrophobic, ionic, and water bridge contacts. (C) A schematic showing detailed interactions between ligand atoms and protein residues. Only interactions occurring more than 30% of the simulation time (0.00 to 200 ns) are displayed. Please click here to view a larger version of this figure.

Table 1: Summary of compound classification criteria for QSAR modeling. A dataset of 94 dihydrofuro[3,4-d] pyrimidine compounds was synthesized by Kang and colleagues⁵⁰ for targeting various HIV strains, namely, NL4-3 wild-type HIV-1, IIIB WT HIV-1, RES056 NNRTI-resistant strain and ROD HIV-2 strain. These pyrimidine derivatives were grouped into four classes according to their target protein for preparation to perform our QSAR training. The table summarizes how these molecules were grouped into four classes according to their experimental EC₅₀ values, which represent the potency of a drug, operationalized as the concentration at which the drug exerts 50% of its maximal effect. Abbreviations: NNRTI = non-nucleotide reverse transcriptase inhibitors; QSAR = quantitative structure-activity relationship.

Table 2: Docking scores of six optimized NNRTIs and the HIV-1 protein. Docking scores are for the six NNRTIs under investigation. The more negative score indicates good binding effectiveness between the ligand and the protein. Efavirenz and Etravirine have shown the most favorable scores at -10.432 eV and -9.647 eV relative to the docking score of the co-crystallized ligand HBY561 (-9.242). Potential ligands were those with a more negative docking score, less than -9.242 eV.

Table 3: 2D-QSAR model statistical parameters. The table presents the standard deviation, the training set correlation between actual and predicted activity values (R²), and the test set actual and predicted activity correlation high score for each class (Q²). The higher score (R²) represents the class.

Table 4: HOMO-LUMO energy gap of the optimized NNRTIs. The table shows the results of the HOMO-LUMO energy gaps obtained after optimization of the six NNRTIs.

Table 5: Predicted activity scores of optimized NNRTIs versus corresponding enumerated NNRTIs. The table compares the predictive activity scores between optimized and enumerated ligands. The higher the activity scores, the better the compound as a potential drug lead.

Table 6: Comparison of the predicted energy gaps of the enumerated NNRTIs and the original energy gap of the optimized NNRTIs. The table shows a comparison of the HOMO-LUMO energy gap between the optimized and enumerated ligands and the energy gap differences between them.

Table 7: Docking score comparison between the original and enumerated NNRTIs. The table shows the comparison between enumerated docking scores, which are the docking scores after enumeration. The original docking scores are the docking scores of the optimized NNRTIs. The redocked enumerated scores are the docking scores of the ligands that were enumerated. Since the enumeration process gives out a predictive docking score, ligands had to be redocked with the same method as the optimized ligands. The added groups are that ones added to the ligands during the enumeration process.

Table 8: MMGBSA reXts of selected ligands on 1HQU. The table represents the Molecular Mechanics with Generalized Born and surface area calculations, which show an average binding free energy (ΔGbind) of the protein-ligand complex. The table shows a comparison between the originally optimized ligands that show good docking scores relative to the crystal ligands and their enumerated counterparts. The chosen compound with a higher docking score and a good molecular dynamic simulation RMSD fluctuation at equilibration should also satisfy the MMGBSA calculations by showing the most negativebinding free energy of (ΔG_bind). In this case, enumerated Etravirine satisfies both computational parameters.

Supplemental File 1: Other rXts obtained in this study. Please click here to download this file.

The DFT²⁶ method was used to analyze the electronic properties and stability of various pyrimidine derivatives, which is not the case in similar studies. The use of DFT allows for a deeper understanding of the molecular interactions at the quantum level, enhancing the design process of NNRTIs. In this study, we selected six NNRTIs with unknown activity values. We retrieved their molecular structures from the PubChem database and conducted geometric optimization using the Minnesota 15 Local (MN15-L)³² functional and 6-31++G (d,p)² basis set in the Gaussian 16 revision C01³³ quantum mechanics software package. To set up simulation input files, we used GaussView 6.0²⁶. The quantum structures were then docked into the active site of the HIV-1 protein.

The study utilized a QSAR approach⁵¹, which is tailored to predict the biological activity of NNRTIs based on their chemical structure. This model is generated using a dataset of 94 dihydrofuro[3,4-d] pyrimidine derivatives, enabling the identification of key functional groups that influence antiviral activity. Incorporating a robust QSAR model is a crucial advancement, as it helps filter compounds early in the design process, enhancing the efficacy of drug discovery.

The predicted activities of the new compounds were considered, and the docking scores were validated. The advanced molecular docking technique explored the binding interactions between NNRTIs and the HIV-1 reverse transcriptase enzyme. This was complemented by molecular dynamics simulations over an extended period (200 ns), providing insight into the stability and behavior of the drug-protein complexes under physiological conditions. Molecular dynamics⁵² validated the findings from docking studies. By simulating the behavior of drug-enzyme complexes over time, we can assess the dynamic stability and interactions that may not be captured in static docking studies⁵³. This approach provides a more realistic assessment of how NNRTIs may perform in biological systems. The simulations showed that Etravirine produced RMSD fluctuations of approximately 4.5 Å, while the enumerated Etravirine gave RMSD fluctuations of 3.5 Å. There were various similarities between Etravirine and enumerated Etravirine, with the enumerated compound possessing enhanced amino acid interactions within the protein's active site. The lower RMSD, improved amino acid interactions, and highest binding free energy all suggest that the enumerated Etravirine could serve as a viable alternative for HIV/AIDS treatment.

The use of enumeration techniques⁵⁴ to generate stereoisomers and perform isostere scanning is a noteworthy aspect of this research. This method allows researchers to explore a wider chemical space for potential NNRTIs by systematically modifying existing structures, which can lead to the discovery of compounds with improved efficacy and reduced resistance.

The research integrates HOMO-LUMO analysis derived from quantum mechanical calculations to evaluate the reactivity and stability of the compounds. This aspect is often overlooked in similar studies⁵⁵, yet it provides valuable insights into electronic properties that can influence biological activity

Using the MMGBSA method in this study to estimate binding free energies of the enumerated NNRTIs as potential HIV-1 reverse transcriptase (RT) inhibitors present a novel aspect that enhances the significance and potential impact of this work on HIV treatment development.

The application of MMGBSA in this study provides a more accurate estimation of the binding free energies for the enumerated Etravirine compared to its counterpart. By calculating binding free energy values, the study quantitatively assesses the binding affinity of enumerated Etravirine, which was considerably higher (9.133 kcal/mol) than that of standard Etravirine. This level of detail in binding energy calculations allows for a more nuanced understanding of how structural modifications can influence the efficacy of NNRTIs, which is often lacking in similar studies that may rely solely on qualitative assessments.

The reported ΔG_Covalent value of 1.477 kcal/mol further indicates the strength of the interaction for the enumerated compound. The comparative approach is innovative as it not only identifies a promising candidate but also situates it within the broader context of existing therapies, providing a benchmark for future NNRTI development.

The comparison shown in Table 8 also indicates that enumerated Etravirine can be a potential compound to be used as an alternative ART as its binding free energy of (ΔG_bind) is the highest as compared to Etravirine, HBY56, Etravirine, and their enumerated counterparts.

Our study's results could significantly impact the development of novel therapeutic agents for HIV infection. Additionally, this approach may allow for identifying the most promising anti-HIV candidate. These findings may pave the way for the rational design of novel HIV-1 NNRTIs.

The incorporation of MMGBSA in conjunction with other computational methods (like molecular docking, QSAR modeling, and molecular dynamics simulations) enhances the robustness of the findings. This integrated approach allows for a comprehensive evaluation of the compounds, from structural optimization to dynamic behavior in a biological context. Such a methodology is relatively novel in NNRTI research, where studies often focus on isolated methods without a holistic view of the drug design process. The research highlights several promising outcomes and an enhanced understanding of the molecular interactions between NNRTIs and the reverse transcriptase enzyme, which can inform future drug design efforts.

While other studies in the field of NNRTI development often focus on single computational methods or basic docking analyses, this research stands out by combining quantum chemical calculations with molecular dynamics and QSAR modelling, employing a systematic enumeration approach that expands the potential chemical space for NNRTIs, and utilizing a comprehensive in-silico framework that not only predicts binding affinities but also assesses the stability and dynamics of drug-enzyme interactions.

Overall, the novelty of this research lies in its multifaceted approach, leveraging advanced computational techniques to address the challenges of drug resistance in HIV treatment, thereby paving the way for the development of more effective NNRTIs. The findings underscore the potential for developing more effective NNRTIs that can overcome the limitations of current therapies.

Some future applications of the approach described in this protocol include the following.

Integration with machine learning

Future studies could integrate MMGBSA with machine learning algorithms to improve the predictive power of binding free energy estimates. By training models on large datasets, researchers could enhance the accuracy and reliability of predictions, enabling better identification of promising NNRTI candidates.

Application in fragment-based drug design

MMGBSA can be effectively utilized in fragment-based drug design, where small molecular fragments are optimized to improve binding affinity. This approach could lead to the identification of novel NNRTIs with enhanced potency and selectivity against HIV-1.

Exploration of drug resistance mechanisms

The technique can be applied to study the binding interactions of NNRTIs with mutant strains of HIV-1. By comparing binding free energies of NNRTIs against both wild-type and resistant variants, researchers can gain insights into the mechanisms of drug resistance and inform the design of next-generation inhibitors.

Evaluation of pharmacokinetic properties

MMGBSA can also be employed to evaluate the pharmacokinetic properties of NNRTIs, such as solubility and permeability. By understanding how these properties correlate with binding affinity, researchers can optimize drug candidates for better therapeutic profiles.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The authors would like to thank the Centre for High-Performance Computing (CHPC) for providing computational resources and the Department of Chemical Science at the University of Johannesburg.