Int J App Pharm, Vol 15, Issue 5, 2023, 171-180Original Article

THE POTENTIAL OF INDONESIAN MARINE NATURAL PRODUCT WITH DUAL TARGETING ACTIVITY THROUGH SARS-COV-2 3CLPRO AND PLPRO: AN IN SILICO STUDIES

CIKAL FIARSI NAHIR1, MASTERIA YUNOVILSA PUTRA2,3, JOKO TRI WIBOWO2, VANNAJAN SANGHIRAN LEE4, ARRY YANUAR*1,3

1Laboratory of Biomedical Computation and Drug Design, Faculty of Pharmacy, Universitas Indonesia, Depok, Jawa Barat-16424, Indonesia. 2Research Center for Vaccine and Drug, National Research and Innovation Agency, Bogor, Jawa Barat, Indonesia. 3National Metabolomics Collaborative Research Center, Faculty of Pharmacy, Universitas Indonesia, Kampus UI Depok, West Java-16424, Indonesia. 4Department of Chemistry, Faculty of Science, University of Malaya, Kuala Lumpur, Malaysia
*Corresponding author: Arry Yanuar; *Email: arry.yanuar@ui.ac.id

Received: 25 May 2023, Revised and Accepted: 10 Jul 2023

ABSTRACT

Objective: This research was conducted to find potential candidate compounds from one hundred thirty-seven Indonesian marine natural products capable of preventing SARS-CoV-2 with a computational approach.

Methods: The physicochemical properties and Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) profile of compounds were predicted using ADMETLab. The candidate compounds were filtered using AutodockVina. Molecular docking was carried out using AutoDockTools on the SARS-CoV-2 3-Chymotrypsin-like protease (3CLpro) and Papain-like protease (PLpro) that is essential for the SARS-CoV-2 life cycle. Also, AMBER22 was used to perform molecular dynamics simulations in this study.

Results: Based on molecular docking results, Pre-Neo-Kaluamine has good activity against 3CLpro with a bond energy value of-10.35 kcal/mol. Cortistatin F showed excellent binding activity on PLpro, with energy value results of-10.62 kcal/mol. Acanthomanzamine C has dual targeting activity and interacts well with protein 3CLpro and PLpro with binding energy values ranging from 10 kcal/mol to 14 kcal/mol.

Conclusion: The molecular docking results were corroborated by molecular dynamics simulation results and showed good stability of the candidate ligands, and we found that there were three potential compounds as protease inhibitors of SARS-CoV-2 including Pre-Neo-Kaluamine for 3CLpro, Cortistatin F for PLpro, and Acanthomanzamine C which had dual targeting activity against both proteases.

Keywords: Marine natural product, In silico, SARS-CoV-2, 3CLpro, PLpro, Virtual screening, ADMET prediction, Molecular docking, Molecular dynamic

© 2023 The Authors. Published by Innovare Academic Sciences Pvt Ltd. This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/)
DOI: https://dx.doi.org/10.22159/ijap.2023v15i5.48416. Journal homepage: https://innovareacademics.in/journals/index.php/ijap

INTRODUCTION

At the end of 2019, there was an outbreak caused by a new variant of the coronavirus known as SARS-CoV-2, which caused a COVID-19 pandemic in several countries around the world, and it became severe infections that increased rapidly and attacked the respiratory system causing pneumonia [1]. Several drugs are used to treat Covid-19, including corticosteroids, lopinavir-ritonavir, hydroxychloroquine, ivermectin, convalescent plasma, Janus-kinase inhibitors, and remdesivir [2]. Besides the therapeutic effects of these drugs, many side effects were also found, especially abnormalities in the immune system, respiratory disorders, signs of inflammation or liver cell damage, and side effects on the neurological system [3]. Based on the latest recommendations from WHO, the combination of nirmatrelvir and ritonavir has shown promising results in Covid-19 patients by inhibiting 3CLprotease from SARS-CoV-2. Still, it has side effects, such as severe drug interactions and contraindications to enzyme-inducing drugs. At first, the shortage of drugs and vaccines designed explicitly for Covid-19 treatment presented a significant obstacle. Therefore, we had to turn to repurposing antiviral drugs and utilizing convalescent plasma as the sole feasible solution [4]. Coronavirus has a protease enzyme system that helps its survival process: Papain-like Protease (PLpro) and 3 Chymotrypsin-like Protease (3CLpro). The PLpro enzyme processes polyprotein chains derived from viral RNA translation, producing functional viral proteins [5]. The 3CLpro enzyme has an essential role in viral genomic replication, so inhibiting this enzyme will significantly affect the blocking of the transcription and replication of coronavirus RNA [6]. Due to several side effects, discovering potential new compounds, such as Covid-19 therapy with fewer side effects, is necessary. Dual targeting of these two proteases will result in more efficient Covid-19 treatment. The discovery of a new drug certainly takes a long time, and it costs a lot to prove a compound has efficacy after passing a series of tests. With the rapid development of technology, Computational techniques are made based on algorithms to produce predictive outputs on the relationship of chemical structure, physicochemical properties, and biological activities from the compound. In computational studies, computer-assisted software simulates the interaction between drug molecules and specific targets. This technique, known as molecular docking, is commonly used in drug discovery and design. It helps identify the best binding model for a ligand to a protein, screen an extensive library of compounds, rank them based on their binding affinity, and propose structural hypotheses about how the ligand inhibits the target [7].

The drug discovery and development process has been accelerating, especially during the pandemic, to prevent and cure Covid-19. In the urgency of this pandemic, natural compounds have been chosen as sources to obtain potential compounds against coronavirus and become an alternative in developing and rapidly discovering medicinal compounds. Natural compounds are reported to have many biological activities that can effectively treat disease and play an essential role in finding and developing innovative drugs during a pandemic [5, 8]. Indonesian marine invertebrates produce potential secondary metabolites. Some of them are sponges, tunicates, and soft corals. The compounds produced have various activities, and one of the activities is antiviral [9, 10]. These secondary metabolites from Indonesian marine invertebrates can be explored and investigated further regarding their biological activity and potential as an antiviral for Covid-19 treatment. Several types of marine invertebrates were also reported to inhibit important proteases of SARS-CoV-2 based on studies [11]. In addition to considering the side effects of drugs that are already available, the target for discovering new drug compounds is also quite important to support the success of drug action on Covid-19. This research will focus on finding and identifying the dual targeting activity of secondary metabolites of Indonesian marine invertebrates against both SARS-CoV-2 protease targets, PLpro and 3CLpro through in silico analysis such as virtual screening with molecular docking method, physicochemical properties prediction, ADMET prediction, and molecular dynamic simulation.

MATERIALS AND METHODS

Protein and ligand structure preparation

The crystal structures of SARS-CoV-2 main protease 3CLpro [12, 13] with PDB ID: 7RFW and the crystal structure of SARS-CoV-2 PLpro with PDB ID: 7TZJ [14, 15] were obtained from rcsb.org the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank [16]. The protein structure was refined by filling in the missing residue using the Modeller program [17, 18]. Meanwhile, the structures of the 137 studied ligands were drawn using MarvinSketch 22.19, 2022, ChemAxon (http://www.chemaxon.com). In addition, co-crystal from 3CLpro (nirmatrelvir) and co-crystal from PLpro (3k, also known as N-[(3-fluorophenyl)methyl]-1-[(1R)-1-naphthalen-1-ylethyl]piperi dine-4-carboxamide) were used as reference ligands which were separated from the protein structure using AutoDockTools [19] were saved in pdb format.

Molecular docking method

In this study, a molecular docking study was carried out using the Autodock Vina [20] accessed from PyRx [21] to screen 137 studied ligands and AutoDockTools to simulate molecular docking. Water molecules are removed before molecular docking, and polar hydrogen atoms are added so that the docking results are not affected by other particles in the protein's surrounding environment. After optimizing the protein and ligands, a grid box uses 40 x 40 x 40 pts for docking with Autodock. The docking process with Autodock Vina uses a grid box with dimensions of 15 x 15 x 15 Å which is the result of the conversion of pts (multiplied by 0.375) with the center grid coordinates adjusted to the position of the co-crystal ligand of each protein, x = 9.9 y = 1.217 z = 20.552 for 3CLpro and x =-3.122 y = 4.041 z =-41.027 for PLpro. For the visualization, BIOVIA Discovery Studio Visualizer [22] software was used to provide information about the complex binding site from protein and the active residues from the compound.

Physicochemical and ADMET prediction

Prediction of the properties and characteristics of the studied ligands was carried out to see the physicochemical properties and how the studied ligands react to Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET). The prediction of physicochemical and ADMET was using ADMETLab 2.0 [23]. Twenty studied ligands from the virtual screening results were converted into SMILES code using Marvin Sketch, and the SMILES code for the co-crystal ligands was obtained from the protein data bank. The web server-based ADMETlab 2.0 was accessed from admetmesh.scbdd.com/service/evaluation. Lipinski’s rule of five (Molecular Weight<500 Da, Log P<5, Hydrogen Bond Acceptor Number<10, Hydrogen Bond Donor Number<5) [24] was used as a reference parameter for having drug-likeness property. Topological Polar Surface Area<140 Å2, Absorption (Human Intestinal Absorption>30%), Distribution (Volume of Distribution 0.04–20 L/kg), Metabolism activity from CYP3A4 and CYP2D6 enzymes, Excretion (T1/2>3h), and Toxicity (Carcinogenicity and AMES Toxicity) those are used as an assessment parameter to predict the drug-likeness of the studied compounds.

Molecular dynamic simulation

The stability of the interaction between ligand and protein can be predicted by conducting molecular dynamics simulations. The simulation process was carried out using AmberTools22 [25]. Before the simulation, The macromolecule files obtained from molecular docking were protonated by adding hydrogen atoms using reduce, an AmberTools subprogram. The ligands from molecular docking were protonated by adding hydrogen atoms (make explicit) using OpenBabel [26]. The ligand was loaded with AM1-BCC and converted to mol2 format using the Antechamber subprogram. The tLeap from AmberTools program created topology and coordinates to check the protein and ligands, combine protein and ligands into a complex, add ionic charges to the complex (Na+or Cl-), and create a solvation complex. UCSF Chimera minimized the protein's energy with default parameters [27]. Molecular Modeling Tool Kit (MMTK) provides energy minimization of molecular models, a Chimera unit. Amber parameters were used for standard residues and Antechamber parameters for non-standard residues. Minimization The minimization stage was set at 1000 steepest descent and used 100 conjugate gradients. The energy minimization process was carried out in two sets with the help of the pmemd. cuda sub-program, the file with the. rst format was utilized in the next stage, heating at 300 K, and energy equilibration was performed in 150-200 picoseconds to determine the system's stability [28]. The production stage was carried out for 100 ns with the pmemd. cuda subprogram. The equilibration stage was analyzed on temperature, density, and total energy parameters. In contrast, production yield analysis was carried out on several parameters such as Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and Molecular Mechanics Poisson-Boltzmann Surface Area/Generalized Born Surface Area (MM-PBSA/MM-GBSA) for the binding free energy from the complex between ligand and protein. The RMSD and RMSF graphics were visualized using QTgrace.

RESULTS AND DISCUSSION

Virtual screening

Table 1: 20 ligands with the best binding affinity when docked with 3CLpro

Ligand Binding affinity (kcal/mol) RMSD (Å)
Acanthomanzamine A -11.51 0.88
Acanthomanzamine B -10.33 2.89
Acanthomanzamine D -10.3 4.1
Mollamide B -9.95 0.76
Naamidine H -9.55 1.38
Pre-neo-kauluamine -9.52 0.7
Acantholactam -9.41 6.12
Lobozoanthamine -9.3 0.7
Epi-tetradehydrohalicyclamine B -9.29 2.77
Naamidine I -9.17 0.92
Jaspamide Q -9.15 1.21
Acanthomanzamine E -9.03 1.86
Sagitol C -9.01 4.17
Chloromethylhalicyclamine B -8.98 1.34
Acanthomanzamine C -8.68 1.04
Crambescidin -8.39 1.43
Sarcofuranocembrenolide A -7.91 0.9
Tetradehydrohalicyclamine B -7.8 1.78
Chloroscabrolide A -7.77 0.98
Sarcofuranocembrenolide B -7.51 2.54

Virtual screening is a computational approach by screening many compounds to obtain candidate compounds with potential activity against a target protein [29]. One hundred thirty-seven compounds from the review by Izzati et al. and Nurrachma et al. [9, 10] were virtually screened. We selected twenty ligands with the best binding affinity values for each target protein, namely 3CLpro (table 1) and PLpro (table 2). Physicochemical properties and ADMET prediction were performed on screened compounds with RMSD values below 2 Å. Therefore, several compounds with RMSD values more than 2 Å on 3CLpro, i.e., Acantholactam (6.12 Å), Acanthomanzamine B (2.89 Å), Acanthomanzamine D (4.1 Å), Epi-tetradehydrohalicyclamine B (2.77 Å), Sagitol C (4.17 Å), and Sarcofuranocembrenolide B (2.54 Å) were not processed to the next stage. While ligands from table 2, i.e., Cortistatin L (3.84 Å), Cortistatin K (3.68 Å), and (-)-Leptoclinidamine B (2.76 Å) were not processed to the next stage.

Table 2: 20 ligands with the best binding affinity when docked with PLpro

Ligand Binding affinity (kcal/mol) RMSD (Å)
Acanthomanzamine A -11.51 0.88
Acanthomanzamine B -10.33 2.89
Acanthomanzamine D -10.3 4.1
Mollamide B -9.95 0.76
Naamidine H -9.55 1.38
Pre-neo-kauluamine -9.52 0.7
Acantholactam -9.41 6.12
Lobozoanthamine -9.3 0.7
Epi-tetradehydrohalicyclamine B -9.29 2.77
Naamidine I -9.17 0.92
Jaspamide Q -9.15 1.21
Acanthomanzamine E -9.03 1.86
Sagitol C -9.01 4.17
Chloromethylhalicyclamine B -8.98 1.34
Acanthomanzamine C -8.68 1.04
Crambescidin -8.39 1.43
Sarcofuranocembrenolide A -7.91 0.9
Tetradehydrohalicyclamine B -7.8 1.78
Chloroscabrolide A -7.77 0.98
Sarcofuranocembrenolide B -7.51 2.54

Table 3: Physicochemical and ADMET prediction results from the candidate of 3CLpro inhibitor

Ligand name MW (Da) LogP HA HD Lipinski

TPSA

(Å2)

 HIA VD (L/Kg) Metabolism T 1/2 Carcinogenicity AMES toxicity
Nirmatrelvir 499.2 1.5 9 3 Yes 131.4 >30% 0.66 CYP2D6 <3h No No
Acanthomanzamine A 545.4 2.2 6 4 Yes 79.2 <30% 4.72 CYP2D6 >3h No No
Acanthomanzamine C 562.3 3.8 6 2 Yes 72.46 >30% 2.97 CYP2D6 <3h No No
Acanthomanzamine E 576.4 5.5 5 1 No 44.39 >30% 2.83 CYP2D6 <3h No No
Chloromethylhalicyclamine B 431.3 1.9 2 0 Yes 3.2 <30% 0.9 CYP2D6 >3h No Yes
Chloroscabrolide A 398.1 1.7 7 0 Yes 91.43 >30% 1.03 CYP3A4 >3h Yes No
Crambescidin 346.3 3.7 5 2 Yes 50.67 >30% 1.18 CYP3A4 >3h Yes No
Jaspamide Q 630.3 4.8 10 4 Yes 140 >30% 0.35 CYP3A4 >3h No No
Lobozoanthamine 497.3 4.1 6 1 Yes 76.1 >30% 1.48 CYP3A4 <3h Yes No
Mollamide B 696.4 3.9 12 5 No 161.9 >30% 1.56 CYP3A4 <3h No No
Naamidine H 493.2 2.6 11 2 Yes 130.7 <30% 0.45 CYP3A4 >3h Yes No
Naamidine I 506.2 2.7 11 3 No 131.4 >30% 0.5 CYP3A4 >3h Yes Yes
Pre-neo-kaluamine 580.3 4.1 7 3 Yes 84.85 >30% 3.18 CYP3A4 <3h No No
Sarcofuranocembrenolide A 388.2 3.9 7 1 Yes 99.11 >30% 0.63 CYP3A4 <3h Yes No
Tetradehydrohalicyclamine B 379.3 3.9 2 0 Yes 7.12 <30% 1.48 CYP2D6 >3h Yes No

Table 4: Physicochemical and ADMET prediction results from the candidate of PLpro inhibitor

Ligand name MW (Da) LogP HA HD Lipinski

TPSA

(Å2)

 HIA VD (L/Kg) Metabolism T 1/2 Carcinogenicity AMES toxicity
3k 390.21 4.29 3 1 Accepted 32.34 >30% 1.56 CYP2D6 <3h No No
Cortistatin J 438.27 5.28 3 0 Accepted 25.36 >30% 3.36 CYP2D6 <3h No No
Cortistatin G 458.33 5.61 3 0 Accepted 25.36 >30% 3.29 CYP2D6 <3h No No
Nakijiquinone V 437.27 3.99 6 4 Accepted 101.7 >30% 1.13 CYP3A4 >3h No No
Cortistatin F 478.39 5.61 3 0 Accepted 15.71 >30% 5.13 CYP2D6 <3h No No
5-Benzoydemethyl-aaptamine 318.1 2.67 5 1 Accepted 71.78 >30% 0.52 CYP2D6 >3h Yes Yes
Cortistatin E 480.41 5.99 3 0 Accepted 15.71 >30% 3.18 CYP2D6 <3h No No
Cortistatin H 460.35 5.99 3 0 Accepted 25.36 >30% 3.04 CYP3A4 <3h No No
Variabine A 322 0.58 8 2 Accepted 118.6 >30% 0.36 CYP2D6 <3h Yes No
Ingenine D 331.14 2.67 6 2 Accepted 86.79 >30% 0.98 CYP3A4 >3h No Yes
Ingenine C 275.12 2.93 5 2 Accepted 69.72 >30% 1.14 CYP3A4 >3h No Yes
Preclathridine B 231.1 1.07 5 2 Accepted 63.03 >30% 1.99 CYP2D6 >3h Yes Yes
Halioxepine 414.28 5.42 4 3 Accepted 69.92 >30% 4.74 CYP2D6 <3h No No
Psammaplysin L 770.81 3.84 11 3 Rejected 133.5 <30% 1.61 CYP3A4 <3h Yes No
Acanthomanzamine C 562.33 3.88 6 2 Accepted 72.46 >30% 2.97 CYP3A4 <3h No No

Physicochemical properties and ADMET prediction

In the drug development and discovery process, drug properties such as pharmacokinetics and pharmacodynamics are crucial parameters for the efficacy of a drug [30]. These properties can be predicted with several online-based or offline-based programs using a computational approach [31]. Lipinski’s rule of five (MW<500 Da, Log P<5, HA<10, HD<5) is used as a reference parameter for having such an excellent drug-likeness property [24]. The result is shown in table 3 and table 4. In this study, ligands that rejected Lipinski’s rule of five and were predicted as toxic compounds were not processed further. Human Intestinal Absorption (HIA) was used to indicate a drug's absorption properties. Compounds showing less than 30% absorbance value are considered poorly absorbed [32].

Compounds with absorbance values above 30% can be adequately absorbed in the human intestine. Poorly absorbed compounds were not continued in the next stage of this study. The Topological Polar Surface Area (TPSA) parameter indicates the amount of surface contribution of a compound tabulated based on its polar fragment. Based on Veber's rule, the optimal TPSA value of a combination is in the range of 0 to 140 Å232. Compounds with TPSA values above 140 Å2 did not proceed to the next stage in this study. The volume of distribution (VD) shows the volume of the drug distributed in the blood plasma. The value of VD was calculated from the amount of drug in the body divided by the blood plasma concentration [33]. The optimum VD values are in the 0.04-20 L/Kg [32]. The predicted metabolic parameters used in this study were CYP3A4 and CYP2D6. CYP3A4 is an enzyme that metabolizes 30%-50% of marketed drugs. Similar to CYP3A4, CYP2D6 is responsible for the metabolism of more than 20% of drug biotransformation [34, 35]. The half-life of a drug involves the clearance and volume of distribution of a drug compound to estimate its excretion properties. Compounds with a half-life of more than 3 h indicate that the compound has a long half-life. Compounds that were predicted to be carcinogenic were not selected for further analysis. In addition to the carcinogenicity toxicology parameter, the AMES parameter was also used as the endpoint for predicting the toxicity of the test compound. The mutagenic effect is related to the carcinogenicity of a compound because this effect affects the level of mutagenicity in humans. Compounds predicted to have mutagenic properties are not processed at the following testing stage [32].

Molecular docking simulation

Validation of molecular docking

The root mean square deviation (RMSD) was applied to confirm the validation of the docking method. The validation of a docking method can be seen from the RMSD value resulting from the re-docking performed on the co-crystal ligand of the protein. The RMSD value states that the docking method's validity is at a value less than or equal to 2 Å. The smaller the RMSD value, it is assumed that the position of the co-crystal ligand before docking doesn’t move significantly so that the docking method can be carried out on the next studied ligand [36, 37]. Re-docking of the co-crystal ligand 3k from PLpro and nirmatrelvir from 3CLpro showed RMSD values of 1.082 Å and 0.490 Å, respectively (fig. 1). It is shown that the docking method was valid for docking simulation for the ligand and the protein.

Fig. 1: Overlay of co-crystal before docking (red) and re-docking co-crystal (yellow) position: (A) Co-crystal (3k) ligand from PLpro; (B) Co-crystal (nirmatrelvir) ligand from 3CLpro

Binding energy and residue interaction

After predicting the physicochemical properties and ADMET profile of the ligands, the selected ligands were docked to their respective receptors to observe further the interaction of the bond forms between the ligands and amino acid residues from the proteins. Table 5 shows the binding energy, Ki, and bond interaction between the amino acid residues from 3CLpro and the ligands.

Table 5: Binding interaction between ligands with 3CLpro from molecular docking simulation

Ligand

Binding

affinity

(kcal/mol)

RMSD

(Å)

Ki

(nM)

 Residue interaction
Hydrogen bond Electrostatic bond
Nirmatrelvir -8.95 1.29 276.21 Phe140; Gln189; Glu166; Pro168; Thr190; Leu167; Cys145; His163. Met165; Met49; His41;
Pre-neo-kaluamine -10.35 0.72 25.85 Thr190; Glu166; Asn142 Met49; Pr168
Acanthomanzamine C -9.52 1.07 105.48 Glu166 His41; Leu27; Cys145
Jaspamide Q -9.15 1.21 196.07 Met49; Glu166; Thr190; His41; Asn142; Asp187; Tyr54; Pro52; His164; Arg188; Leu167; Gln192; Ala191; Met165; Pro168; Met49;

Binding energy from the studied ligands was more excellent than the reference ligand nirmatrelvir (-8.95 kcal/mol). The binding energy value of pre-neo-kaluamine, Acanthomanzamine C, and Jaspamide Q were-10.35,-9.52, and-9.15 kcal/mol, respectively. The inhibition constant (Ki) obtained from docking results means less Ki value means much great inhibition activity from the ligand. Nirmatrelvir showed activity as a peptidomimetic inhibitor of 3CLpro [38]. The nitrile function of nirmatrelvir’s structure is bound to the active catalytic site (cysteine and histidine residue) from 3CLpro [39]. The hydrogen from the carboxamide group of nirmatrelvir formed several hydrogen bonds between the Glu166, Cys145, His163, and Thr190 residues (fig. 2).

While the carbon atoms formed an alkyl bond with Met49 and a pi-alkyl bond with Met165 and His41, pre-neo-kaluamine showed several interactions with amino acid residues of 3CLpro, including the hydrogen bonds between the hydrogens on the NH and OH groups in the ligand with the amino acid residues Thr190 and Glu166 (fig. 3).

On the other side of the ligand, the hydrogen group on the CH molecule formed a hydrogen bond with Asn142. Other bonds between ligands and amino acids (Pro168, Met49) were pi-sigma, alkyl, and pi-alkyl Pre-neo-kaluamine did not bind to the catalytic site of 3CLpro (Cys145 and His41). Acanthomanza-mine C interacted with amino acid residues in the active catalytic site of 3CLpro (fig. 4). Acanthomanzamine C created a bond with a pi-pi T-Shaped between the CH molecule from the structure and Cys145 and a pi-cation bond with His41 and CH molecule from the ligand structure, which is a critical amino acid in the catalytic site of 3CLpro.

Fig. 2: Interaction between nirmatrelvir and 3CLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Fig. 3: Interaction between pre-neo-kaluamine and 3CLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Fig. 4: Interaction between Acanthomanzamine C and 3CLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Fig. 5: Interaction between Jaspamide Q and 3CLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Other bonds were a hydrogen bond with Glu166 and a pi-Alkyl bond with Leu27 residue. Leu27, Residue Cys145, and His41 formed bonds with the benzene structure from the ligand. The benzene ring of a com-pound significantly increases the inhibitory properties due to the bulkiness of the benzene structure, which allows for many bonds to occur with the amino acids of the target protein [40]. Jaspamide Q interacted with amino acid residues in the catalytic active site of 3CLpro (fig. 5).

Created a conventional hydrogen bond with a Met49 that interacts with NH molecules from Jaspamide Q, Glu166 interacts with NH and CO molecules simultaneously, and Thr190 interacts with CO molecules from the ligand. Also, several van der Waals bonds were formed between Asn142, Asp187, Tyr54, Pro52, His164, Arg188, Leu167, Gln192, Ala191, and Met165. His41 residue forms a Pi-cation bond with a benzene structure from the ligand. Besides that, Met165 and Pro168 form Pi-Alkyl bonds. Table 6 shows the binding energy, Ki, and bond interaction between the amino acid residues from PLpro and the ligands 3k as a reference ligand has a-9.47 kcal/mol binding energy value.

The other ligands showed greater values than the reference ligand, i.e.,-9.82 kcal/mol (cortistatin J),-9.63 kcal/mol (acanthomanzamine C), and-10.62 kcal/mol (cortistatin F). The inhibition constant (Ki) of the ligands cortistatin J, acanthomanzamine C, cortistatin F, and 3k were 63.39, 16.41, 105.98, and 113,67 nM, respectively. Tyr268 and Gln269 are essential amino acid residues at PLpro catalytic sites where small molecules bind. Tyr268 and Gln269 are the main sites of the flexible loop, which can generate several backbone and side chain conformations [41]. The interaction shown in fig. 6 involved 3k (co-crystal ligand of PLpro) forming hydrogen bonds with the amino acid residue Tyr268, an amino acid residue located at the catalytic site of PLpro. The hydrogen bond comes from the interaction between NH molecules in the 3k ligand structure. In addition, bonds were formed in the naphthalene structure of the 3k ligand, which were pi-sulfur bonds with Met208 residues, pi-anion bonds with Asp164, Pi-sigma Bonds with Tyr264, Amide-Pi Stacked bonds with Ala246, Pi-Alkyl bonds with Pro247, and Pro248.

Cortistatin J formed several binding interactions with amino acid residues in the 3CLpro protein (fig. 7). Cortistatin J formed a hydrogen bond with the amino acid Leu162. Other bonds were Alkyl bonds with amino acids Pro248, Pro247, Tyr264, and Tyr268 which bind to the region of the ring-like ligand structure, and pi-pi T-Shaped bond between the ring structure of the ligand and the amino acid residue Tyr268, which is the amino acid residue on the active site of PLpro. Compounds containing bonds with the Pi character group can provide a variety of beneficial steric and electronic properties [42].

Table 6: Binding interaction between ligands with PLpro from molecular docking simulation

Ligand

Binding affinity

(kcal/mol)

RMSD

(Å)

Ki

(nM)

 Residue interaction
Hydrogen bond Electrostatic bond
Cortistatin F -10.62 1.68 16.41 Glu167 Leu162; Pro248; Pro247; Tyr264; Tyr273; Tyr268
Acanthomanzamine C -9.63 1.96 105.98 Tyr268; Gln269 Leu162; Tyr264; Asp164
Cortistatin J -9.82 0.32 63.39 Leu162 Tyr264; Tyr268; Pro248; Pro247
3k -9.47 1.48 113.67 Tyr268 Tyr264; Pro248; Pro247; Ala246; Met208; Asp164

Fig. 6: Interaction between 3k and PLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Fig. 7: Interaction between Cortistatin J and PLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Cortistatin F created hydrogen bonds with Glu167 amino acid residues (fig. 8). Other formed bonds were Alkyl and Pi-Alkyl bonds with amino acid residues Leu162, Pro248, Tyr268, Tyr273, and Pro247. In addition, the pi-sigma bond was formed between the Tyr264 amino acid residue and the CH molecule in the structure of the Cortistatin F.

Fig. 8: Interaction between cortistatin F and PLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Acanthomanzamine C formed hydrogen bonds with residues Gln269 and Tyr268 that interact with the NH group of Acanthomanzamine C (fig. 9). These two residues are key residues that occupy the active site of the PLpro receptor. Other bonds were pi-pi T-shaped with Tyr264 residues, Pi-Anion bonds with Asp164 residues, and Pi-Alkyl bonds with Leu162 residue.

Fig. 9: Interaction between Acanthomanzamine C and PLpro after docking molecular simulation: (A) 2-Dimension visual interaction (B) 3-Dimension visual interaction

Molecular dynamic simulation

Molecular dynamics simulations provide a stable interaction between the behavior of proteins and other small molecules in detail so the interaction stability can be appropriately described [43]. With the accuracy and accessibility of simulation, followed by experimental structural data, it can be used well and help practical work that requires a lot of time and effort [44]. The molecular dynamics simulation stages were carried out to observe conformational changes formed between ligands and macromolecules (proteins). The process of molecular dynamics was broadly divided into two phases, namely the equilibration process, where the state of the ligand was analyzed at a constant temperature of 300 K, as well as the total energy and density. The production process was carried out for 100 ns by analyzing RMSD, RMSF, the total energy of MM-PBSA/GBSA, and the hydrogen bond.

Analysis of system equilibration

Temperature contained in the two protease systems was gradually increased until 300 K at the end of the equilibrium simulation, which lasted up to 140 ps. Temperature is one of the analytical parameters in the equilibrium simulation to check whether the temperature contained in the system can maintain its stability. The total energy is one of the analytical parameters in checking the equilibrium system before starting the production stage. The total energy is expected to remain stable, and no significant fluctuations will occur. The following process matches the density value, which is expected to be stable at approximately 1.00 g/ml. The density was stabilized at about 1.00 g/ml for both systems, representing the density of water as a solvent.

Analysis of the production stage

After completing a thorough check of the system equilibrium, the production stage of molecular dynamics simulation was carried out. The production stage lasted for 100 ns, during which RMSD, RMSF, total energy of MM-PBSA/GBSA, and hydrogen bonds were analyzed. Fig. 10 shows the fluctuation of the RMSD (Root Mean Square Deviation) ligand value during the simulation process. All ligands fluctuated in the initial 1-10 ns to find a stable location. RMSD values of several ligands range from<2 Å to 4 Å.

Fig. 10: Trajectory analysis and visualization of the root mean square deviation (RMSD) results in a 100 ns simulation: (A) Candidate inhibitor with target protein 3CLpro (B) Candidate inhibitor with target protein PLpro

The protein 3CLpro, nirmatrelvir, acanthomanzamine C, and jaspamide Q achieved stability at 2 Å to 3 Å. As for the pre-neo-kaluamine, it fluctuated quite significantly between 3 Å to 3.5 Å. In the PLpro protein, the 3k ligand has better stability than other ligands. 3k was stabled at 2 Å. The cortistatin F and cortistatin J fluctuated until they stabilized at 2.5 Å to 3 Å. At the same time, the acanthomanzamine C ligand fluctuated significantly between 10 ns to 30 ns and reached stability at 3 Å. Although RMSD values seem pretty stable at the beginning of the simulation, the ligand fluctuated to seek conformational stability. Besides referring to the RMSD analysis, RMSF parameters also need to be considered to see fluctuations in the flexible amino acid residues of the target protein. The graph in fig. 11 shows the evaluation results of the target protein amino acid residues during the simulation process. Most of the amino acid residues of the protein did not experience significant fluctuations. It indicated no significant change in the protein structure during the simulation.

Fig. 11: Visualization and analysis results of root mean square fluctuation (RMSF) during the simulation for each residue number of (A) 3CLpro and (B) PLpro

Fig. 12: Number of hydrogen bonds from candidate ligand to target protein: (A) 3CLpro (B) PLpro during 100 ns simulation

The last parameter analyzed was the hydrogen bond between the candidate ligand and the target receptor. Hydrogen bonds are essential in molecular systems as a significant non-covalent structural force, mainly determined by electrostatic properties. Hydrogen bond formation occurs when a hydrogen atom is effectively shared between a donor atom and another acceptor atom. The general rule used to determine the probability of hydrogen bonding was based on FON, meaning hydrogen bonded to F, O, and N atoms could act as hydrogen bond donors, and F, O, and N atoms could act as hydrogen bond acceptors. Hydrogen bond analysis was carried out with the VMD program to process the topology and coordinates of candidate ligands so that hydrogen bonds could be analyzed. According to Bakker and Hubbard, several essential parameters in the shielding of hydrogen bonds were the donor and acceptor atoms being separated by a distance of 3.5 Å and the angle having to be more than 100 °. Even so, there were variations between the distance and the angle of each type of hydrogen bond that could form, but most studies used 120° as a parameter of the hydrogen bond angle [45]. The graph visualizes the number of hydrogen bonds formed in fig. 12.

The results showed a correlation between the calculation results of MMGBSA and MMPBSA derived from molecular dynamics simulations (table 7). The table shows that the Pre-neo-kaluamine has a higher binding energy value (-70.17 kcal/mol and-59.63 kcal/mol) than other ligands. Nirmatrelvir as a reference ligand showed binding energy values of-46.30 kcal/mol and-42.36 kcal/mol. In the PLpro protein, all ligands correlated with the previous molecular docking simulation results.

Table 7: MM-generalized born surface area and MM-poisson boltzmann surface area binding energy calculation results

Protein Ligand MMGBSA (kcal/mol) MMPBSA (kcal/mol)
Mean (± SEM) Mean (± SEM)
3CLpro Nirmatrelvir -46.30 (± 0.374) -42.36 (± 0.464)
Pre-neo-kaluamine -70.17 (± 0.378) -59.63 (± 0.524)
Jaspamide Q -23.36 (± 0.337) -21.29 (± 0.420)
Acanthomanzamine C -14.11 (± 0.276) -10.09 (± 0.464)
PLpro 3k -11.62 (± 0.409) -10.81 (± 0.391)
Cortistatin F -23.67 (± 0.502) -19.27 (± 0.540)
Cortistatin J -17.49 (± 0.213) -16.20 (± 0.244)
Acanthomanzamine C -13.29 (± 0.371) -13.58 (± 0.372)

Cortistatin F (with binding energy values-23.67 kcal/mol and-19.27 kcal/mol) has a better binding energy value than the others, even if it compares with the 3k as a reference ligand (with binding values-11.62 and-10.81 kcal/mol). Acanthomanzamine C interacts well with protein 3CLpro and PLpro, with binding energy values ranging from-10 kcal/mol to 14 kcal/mol.

CONCLUSION

Main protease (3CLpro) and papain-like protease (PLpro) are two critical non-structural proteins in the transcription process of SARS-CoV-2, so they were chosen as targets of the candidate compounds in this study. Several secondary metabolite compounds from Indonesian marine natural products have the potential to inhibit 3CLpro and PLpro. Based on the molecular docking simulation results between candidate compounds and 3CLpro protein, Pre-neo-kaluamine had the best inhibitory potential with values of-10.35 kcal/mol and Ki 25.85 nM. In the PLpro protein, Cortistatin F has the best binding energy value of-10.62 kcal/mol with a Ki of 16.41 nM. Also, we found that acanthomanzamine C is a ligand with dual targeting activity and shows a good interaction with protein 3CLpro and PLpro with binding energy values ranging from-10 kcal/mol-14 kcal/mol. All ligands showed better results than reference ligands. Molecular dynamics simulations also showed the stability of the ligand with calculated values that correlate with the molecular docking simulation results.

ACKNOWLEDGMENT

The computation in this work has been done using the facilities of MAHAMERU BRIN HPC, National Research and Innovation Agency of Indonesia (BRIN). Also, we would like to thank Alfrina Irene from the Faculty of Pharmacy University of Indonesia for providing us with the 3-Dimension structure of the Indonesian marine invertebrate's secondary metabolites used in this study.

FUNDING

This research was funded by The Directorate of Research and Development, Universitas Indonesia, under Hibah PUTI 2022 (Grant No. NKB-069/UN2. RST/HKP.05.00/2022).

AUTHORS CONTRIBUTIONS

The manuscript was written through the contributions of all authors, and all authors have approved the final version.

CONFLICT OF INTERESTS

The authors declare no conflict of interest.

REFERENCES

  1. Chen N, Zhou M, Dong X, Qu J, Gong F, Han Y. Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study. Lancet. 2020 Feb 15;395(10223):507-13. doi: 10.1016/S0140-6736(20)30211-7, PMID 32007143.

  2. Sukumaran S, Sathianarayanan S. A review on covid-19 pandemic a global threat-current status and challenges and preventive strategies. Int J Appl Pharm. 2021 Sep 7:10-4.

  3. Mohammad Zadeh N, Mashinchi Asl NS, Forouharnejad K, Ghadimi K, Parsa S, Mohammadi S. Mechanism and adverse effects of Covid-19 drugs: a basic review. Int J Physiol Pathophysiol Pharmacol. 2021;13(4):102-9. PMID 34540130.

  4. Jain NK, Agrawal A, Kulkarni GT, Tailang M. Molecular docking study on phytoconstituents of traditional ayurvedic drug TULSI (ocimum sanctum linn.) against covid-19 MPRO enzyme: an in silico study. Int J Pharm Pharm Sci. 2022 Apr 1:44-50. doi: 10.22159/ijpps.2022v14i4.43181.

  5. Rajpoot S, Alagumuthu M, Baig MS. Dual targeting of 3CLpro and PLpro of SARS-CoV-2: a novel structure-based design approach to treat covid-19. Curr Res Struct Biol. 2021 Jan 1;3:9-18. doi: 10.1016/j.crstbi.2020.12.001, PMID 33319212.

  6. Tumskiy RS, Tumskaia AV. Multistep rational molecular design and combined docking for discovery of novel classes of inhibitors of SARS-CoV-2 main protease 3CLpro. Chem Phys Lett. 2021 Oct 1;780:138894. doi: 10.1016/j.cplett.2021.138894, PMID 34276059.

  7. Mourad T, Alahmad S. A computational study of ciprofloxacin metabolites and some natural compounds against resistant methicillin staphylococcus aureus (MRSA). Int J Pharm Pharm Sci. 2022 Aug 1;14(8):22-8. doi: 10.22159/ijpps.2022v14i8.44560.

  8. Ali S, Alam M, Khatoon F, Fatima U, Elasbali AM, Adnan M. Natural products can be used in therapeutic management of covid-19: probable mechanistic insights. Biomed Pharmacother. 2022 Mar 1;147:112658. doi: 10.1016/j.biopha.2022.112658, PMID 35066300.

  9. Izzati F, Warsito MF, Bayu A, Prasetyoputri A, Atikana A, Sukmarini L. Chemical diversity and biological activity of secondary metabolites isolated from Indonesian marine invertebrates. Molecules. 2021;26(7):1898. doi: 10.3390/molecules26071898, PMID 33801617.

  10. Nurrachma MY, Sakaraga D, Nugraha AY, Rahmawati SI, Bayu A, Sukmarini L. Cembranoids of soft corals: recent updates and their biological activities. Nat Prod Bioprospect. 2021 Jun 22;11(3):243-306. doi: 10.1007/s13659-021-00303-2, PMID 33890249.

  11. Syahputra G, Gustini N, Bustanussalam B, Hapsari Y, Sari M, Ardiansyah A. Molecular docking of secondary metabolites from Indonesian marine and terrestrial organisms targeting SARS-CoV-2 ACE-2, M pro, and PL pro receptors. Pharmacia. 2021;68(3):533-60. doi: 10.3897/pharmacia.68.e68432.

  12. Owen DR, Allerton CMN, Anderson AS, Aschenbrenner L, Avery M, Berritt S. An oral SARS-CoV-2 Mpro inhibitor clinical candidate for the treatment of covid-19. Science. 2021;374(6575):1586-93. doi: 10.1126/science.abl4784, PMID 34726479.

  13. RCSB PDB-7RFW: structure of SARS-CoV-2 main protease in complex with a covalent inhibitor. Available from: https://www.rcsb.org/structure/7RFW. [Last accessed on 17 Nov 2022]

  14. Calleja DJ, Kuchel N, Lu BGC, Birkinshaw RW, Klemm T, Doerflinger M. Insights into drug repurposing, as well as specificity and compound properties of piperidine-based SARS-CoV-2 PLpro inhibitors. Front Chem. 2022 Apr 12;10:861209. doi: 10.3389/fchem.2022.861209, PMID 35494659.

  15. RCSB PDB-7TZJ: SARS CoV-2 PLpro in complex with inhibitor 3k. Available from: https://www.rcsb.org/structure/7TZJ. [Last accessed on 17 Nov 2022]

  16. Burley SK, Bhikadiya C, Bi C, Bittrich S, Chen L, Crichlow GV. RCSB protein data bank: powerful new tools for exploring 3D structures of biological macromolecules for basic and applied research and education in fundamental biology, biomedicine, biotechnology, bioengineering and energy sciences. Nucleic Acids Res. 2021 Jan 8;49(D1):D437-51. doi: 10.1093/nar/gkaa1038, PMID 33211854.

  17. Webb B, Sali A. Comparative protein structure modeling using modeller. CP in Bioinformatics. 2016;54(1). doi: 10.1002/cpbi.3.ncbi.nlm.nih.gov/27322406.

  18. Fiser A, Sali A. ModLoop: automated modeling of loops in protein structures. Bioinformatics. 2003 Dec 12;19(18):2500-1. doi: 10.1093/bioinformatics/btg362, PMID 14668246.

  19. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS. AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J Comput Chem. 2009 Dec;30(16):2785-91. doi: 10.1002/jcc.21256, PMID 19399780.

  20. Eberhardt J, Santos Martins D, Tillack AF, Forli S. AutoDock Vina 1.2.0: New docking methods, expanded force field, and python bindings. J Chem Inf Model. 2021 Aug 23;61(8):3891-8. doi: 10.1021/acs.jcim.1c00203, PMID 34278794.

  21. Dallakyan S, Olson AJ. Small-molecule library screening by docking with PyRx. Methods Mol Biol. 2015;1263:243-50. doi: 10.1007/978-1-4939-2269-7_19, PMID 25618350.

  22. Biovia. Biovia discovery studio V.21.1.0. Dassault Syst. 2021.

  23. Xiong G, Wu Z, Yi J, Fu L, Yang Z, Hsieh C. ADMET lab 2.0: an integrated online platform for accurate and comprehensive predictions of ADMET properties. Nucleic Acids Res. 2021 Jul 2;49(W1):W5-W14. doi: 10.1093/nar/gkab255, PMID 33893803.

  24. Lipinski CA. Chapter 27. Bioisosterism in drug design. In: Annual reports in medicinal chemistry; 1986. p. 283-91. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0065774308611379. [Last accessed on 29 Apr 2021]

  25. Case DA, Aktulga HM, Belfon K, Ben-Shalom IY, Berryman JT, Brozell SR. San Fran: University of California; 2022.

  26. O’Boyle NM, Banck M, James CA, Morley C, Vandermeersch T, Hutchison GR. Open Babel: an open chemical toolbox. J Cheminform. 2011 Dec 7;3(1):33. doi: 10.1186/1758-2946-3-33, PMID 21982300.

  27. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC. UCSF Chimera-a visualization system for exploratory research and analysis. J Comput Chem. 2004 Oct;25(13):1605-12. doi: 10.1002/jcc.20084, PMID 15264254.

  28. AMBER advanced tutorial 33-constant pH and redox potential MD: Section 2. Available from: https://ambermd.org/tutorials/advanced/tutorial33/section2. htm. [Last accessed on 07 Jan 2023]

  29. Kondapuram SK, Sarvagalla S, Coumar MS. Docking-based virtual screening using PyRx tool: autophagy target Vps34 as a case study. In: Molecular docking for computer-aided drug design. Elsevier; 2021. p. 463-77.

  30. Daina A, Michielin O, Zoete V. Swiss ADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep. 2017;7(1):42717. doi: 10.1038/srep42717, PMID 28256516.

  31. Ma S, McGann M, Enyedy IJ. The influence of calculated physicochemical properties of compounds on their ADMET profiles. Bioorg Med Chem Lett. 2021 Mar;36:127825. doi: 10.1016/j.bmcl.2021.127825, PMID 33508464.

  32. Xiong G, Wu Z, Yi J, Fu L, Yang Z, Hsieh C. ADMET lab 2.0: an integrated online platform for accurate and comprehensive predictions of ADMET properties. Nucleic Acids Res. 2021 Jul 2;49(W1):W5-W14. doi: 10.1093/nar/gkab255, PMID 33893803.

  33. Alsanosi SMM, Skiffington C, Padmanabhan S. Pharmacokinetic pharmacogenomics. Handbook of pharmacogenomics and stratified medicine; 2014 Jan 1. p. 341-64.

  34. Urquhart BL, Nolin TD. Drug metabolism in chronic kidney disease. Chronic Ren Dis. 2019 Jan 1:1035-51.

  35. Vuppalanchi R. Metabolism of drugs and xenobiotics. Practical hepatic pathology: a diagnostic approach a volume in the pattern recognition series, expert consult; 2011 Feb 4. p. 45-52.

  36. Mardianingrum R, Endah SRN, Suhardiana E, Ruswanto R, Siswandono S. Docking and molecular dynamic study of isoniazid derivatives as anti-tuberculosis drug candidate. Chem Data Coll. 2021 Apr 1;32:100647. doi: 10.1016/j.cdc.2021.100647.

  37. Ruswanto R, Nofianti T, Mardianingrum R, Kesuma D, Siswandono. Design, molecular docking, and molecular dynamics of thiourea-iron (III) metal complexes as NUDT5 inhibitors for breast cancer treatment. Heliyon. 2022 Sep 1;8(9):e10694. doi: 10.1016/j.heliyon.2022.e10694, PMID 36177227.

  38. Vangeel L, Chiu W, De Jonghe S, Maes P, Slechten B, Raymenants J. Remdesivir, Molnupiravir and Nirmatrelvir remain active against SARS-CoV-2 omicron and other variants of concern. Antiviral Res. 2022 Feb 1;198:105252. doi: 10.1016/j.antiviral.2022.105252, PMID 35085683.

  39. Ullrich S, Ekanayake KB, Otting G, Nitsche C. Main protease mutants of SARS-CoV-2 variants remain susceptible to nirmatrelvir. Bioorg Med Chem Lett. 2022 Apr 15;62:128629. doi: 10.1016/j.bmcl.2022.128629, PMID 35182772.

  40. Das A, Dasgupta S, Pathak T. Crescent-shaped meta-substituted benzene derivatives as a new class of non-nucleoside ribonuclease a inhibitors. Bioorg Med Chem. 2022 Oct 1;71:116888. doi: 10.1016/j.bmc.2022.116888, PMID 35944385.

  41. Ismail MI, Ragab HM, Bekhit AA, Ibrahim TM. Targeting multiple conformations of SARS-CoV2 papain-like protease for drug repositioning: an in-silico study. Comput Biol Med. 2021 Apr 1;131:104295. doi: 10.1016/j.compbiomed.2021.104295, PMID 33662683.

  42. Jin MY, Zhen Q, Xiao D, Tao G, Xing X, Yu P. Engineered non-covalent π interactions as key elements for chiral recognition. Nat Commun. 2022;13(1):3276. doi: 10.1038/s41467-022-31026-8, PMID 35672365.

  43. Harvey MJ, De Fabritiis G. High-throughput molecular dynamics: the powerful new tool for drug discovery. Drug Discov Today. 2012 Oct 1;17(19-20):1059-62. doi: 10.1016/j.drudis.2012.03.017, PMID 22504137.

  44. Hollingsworth SA, Dror RO. Molecular dynamics simulation for all. Neuron. 2018 Sep 9;99(6):1129-43. doi: 10.1016/j.neuron.2018.08.011, PMID 30236283.

  45. Tan KP, Singh K, Hazra A, Madhusudhan MS. Peptide bond planarity constrains hydrogen bond geometry and influences secondary structure conformations. Curr Res Struct Biol. 2021 Jan 1;3:1-8. doi: 10.1016/j.crstbi.2020.11.002, PMID 34382009.