Elucidation of Teicoplanin Interactions with Drug Targets Related to COVID-19

Teicoplanin is a glycopeptide antibiotic effective against several bacterial infections, has exhibited promising therapeutic efficiency against COVID-19 in vitro, and the rationale for its use in COVID-19 is yet to be recognized. Hence, in this study a number of molecular modeling techniques were employed to decrypt the mechanistic insight of teicoplanin interaction with several COVID-19 drug targets. Initially, molecular docking was employed to study the teicoplanin interaction with twenty-five SARS-CoV-2 structural and non-structural proteins which was followed by molecular mechanics/generalized Born surface area (MM/GBSA) computation for binding energy predictions of top ten models from each target. Amongst all macromolecular targets, the N-terminal domain of the nucleocapsid protein displayed the strongest affinity with teicoplanin showing binding energies of −7.4 and −102.13 kcal/mol, in docking and Prime MM/GBSA, respectively. Thermodynamic stability of the teicoplanin-nucleocapsid protein was further probed by molecular dynamics simulations of protein–ligand complex as well as unbounded protein in 100 ns trajectories. Post-simulation MM-GBSA computation of 50 frames extracted from simulated trajectories estimated an average binding energy of −62.52 ± 12.22 kcal/mol. In addition, conformational state of protein in complex with docked teicoplanin displayed stable root-mean-square deviation/fluctuation. In conclusion, computational investigation of the potential targets of COVID-19 and their interaction mechanism with teicoplanin can guide the design of novel therapeutic armamentarium for the treatment of SARS-CoV-2 infection. However, additional studies are warranted to establish the clinical use or relapses, if any, of teicoplanin in the therapeutic management of COVID-19 patients.


Introduction
The outbreak of COVID-19 pandemic in China, caused by severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) primarily spreads through close contact amongst people by sneezing, coughing, or by communicating verbally. The primary symptoms of the disease are fever, cough, fatigue, shortness of breath, and loss of smell. The complications of this disease include viral pneumonias, respiratory distress, and hypoxia [1]. The world is suffering due to lack of effective medicine. Therefore, there is an urgent need to figure out efficient medicine through drug discovery efforts [2,3].
Drug repurposing strategy is often used for recognizing novel usages of approved or investigational drugs because of several advantages such as lower risk of failure, reduced time for drug development, and lower investment needs [4]. Therefore, numerous antimalarials, antibacterials, antiparasitic agents, and antivirals are currently prevalent in pre-clinical as well as clinical investigations aimed at developing COVID-19 treatment by a drug repurposing approach [5]. In particular, a widely available FDA-approved antibiotic, teicoplanin (Figure 1), is among the molecule of interest as probable COVID-19 medicine. It belongs to the glycopeptide class of antibiotic having low toxicity profile in humans and hence routinely used in clinical practice for the treatment of bacterial infections. Interestingly, it has demonstrated antiviral efficacy against several kinds of viruses such as Computer-aided drug design methods are extensively employed in drug design and discovery projects owing to several advantages such as rapid development process and reduced cost [10][11][12]. In particular, molecular docking coupled with molecular dynamics simulation studies are intended to decipher the mechanism of binding interactions at the molecular levels [13]. Rapid mechanistic insight is vital for understanding structure-activity relationship and leading to the optimization of the design and discovery of potential Computer-aided drug design methods are extensively employed in drug design and discovery projects owing to several advantages such as rapid development process and reduced cost [10][11][12]. In particular, molecular docking coupled with molecular dynamics simulation studies are intended to decipher the mechanism of binding interactions at the molecular levels [13]. Rapid mechanistic insight is vital for understanding structureactivity relationship and leading to the optimization of the design and discovery of potential molecules [14][15][16]. In this study, several computational techniques such as molecular Antibiotics 2021, 10, 856 3 of 16 docking, molecular mechanics/generalized Born surface area (MM-GBSA), and molecular dynamics simulation were exploited to inspect the binding interactions between teicoplanin and potential drug targets associated with SARS-CoV-2. The study is envisioned to assist in finding potential leads and accelerating drug development process for the treatment of novel coronavirus, COVID-19.

Experimental Section
The methodology adopted in this study has been outlined as flowchart in Figure 2.
Antibiotics 2021, 10, x FOR PEER REVIEW 3 of 16 molecules [14][15][16]. In this study, several computational techniques such as molecular docking, molecular mechanics/generalized Born surface area (MM-GBSA), and molecular dynamics simulation were exploited to inspect the binding interactions between teicoplanin and potential drug targets associated with SARS-CoV-2. The study is envisioned to assist in finding potential leads and accelerating drug development process for the treatment of novel coronavirus, COVID-19.

Experimental Section
The methodology adopted in this study has been outlined as flowchart in Figure 2. An outline of the adopted methodology in this study. B3LYP: Becke's three-parameter hybrid model, Lee-Yang-Parr correlation functional method; DFT: density functional theory; MM/GBSA: molecular mechanics/generalized Born surface area; ns: nano seconds; RMSD: root-mean-square deviation; RMSF: root-mean-square fluctuation.

Ligand Preparation
The molecular structure of teicoplanin was retrieved from PubChem database as two-dimensional coordinate in sdf format and converted to its three-dimensional conformation by means of Open Babel program [17]. Jaguar v10.9 of Schrodinger Suites 2020-3 was used for density functional theory (DFT)-based optimization by Becke's three-parameter hybrid model, the Lee-Yang-Parr (B3LYP) method at the level of 6-31G [18,19]. The optimized structure of teicoplanin is presented in Figure 3. Prepare_ligand4.py module of MGL Tools 1.5.6 [20] was used for the preparation of ligand in pdbqt format after merging all non-polar hydrogens, defining rotatable bonds/torsion tree and adding Gasteiger charges.

Ligand Preparation
The molecular structure of teicoplanin was retrieved from PubChem database as twodimensional coordinate in sdf format and converted to its three-dimensional conformation by means of Open Babel program [17]. Jaguar v10.9 of Schrodinger Suites 2020-3 was used for density functional theory (DFT)-based optimization by Becke's three-parameter hybrid model, the Lee-Yang-Parr (B3LYP) method at the level of 6-31G [18,19]. The optimized structure of teicoplanin is presented in Figure 3. Prepare_ligand4.py module of MGL Tools 1.5.6 [20] was used for the preparation of ligand in pdbqt format after merging all non-polar hydrogens, defining rotatable bonds/torsion tree and adding Gasteiger charges. Teicoplanin structure optimized by density functional theory (DFT) with Becke's three-parameter hybrid model, the Lee-Yang-Parr (B3LYP) correlation functional method at 6-31G level (Jaguar v10.9 software of Schrodinger Suites 2020-3 was used). Color code: carbon-dark gray, hydrogen-light gray, nitrogen-blue, oxygen-red, and chlorinegreen.

Molecular Docking Simulation
The grid box of 30 × 30 × 30 size in x, y, and z directions with 1-Å spacing was constructed around the active site residues as demarcated by the bound co-crystallized ligands (see Table 1 for grid points used for each receptor). However, blind docking, covering the entire macromolecular target, was opted if binding site information was unavailable. AutoDock Vina 1.1.2 [23] was used for molecular docking employing default parameters for flexible ligand and rigid protein. Moreover, the exhaustiveness parameter was set to 12. Upon successful completion of docking simulation, ten best poses were individually analyzed for intermolecular interactions using PyMol 2.4.1 [22] and Biovia Discovery Studio Visualizer 2020 [16,21,24] and further subjected to MM/GBSA computations in the next step.

Prime MM-GBSA Calculations
MM-GBSA technique are frequently used to rationalize the findings of molecular docking and virtual screening. Prime MM-GBSA module of Schrödinger Suite 2020-3 [25][26][27] was used to compute binding energy and relevant parameters as presented in Table 2. The Prime MM-GBSA protocol amalgamates OPLS molecular mechanics energies, a VSGB solvation model for polar solvation (G SGB ), and a nonpolar solvation expression (G NP ) involving nonpolar solvent-accessible surface area (SASA) and van der Waals interactions [28]. The binding-free energy (∆G bind ) of docked teicoplanin in complex with respective proteins was calculated using the following equation [29]: where ∆E MM is the difference in energy between the complex structure and the sum of the energies of the protein with and without teicoplanin, ∆G solv is the difference in the GBSA solvation energy of the teicoplanin-protein complex and the sum of the solvation energies for the teicoplanin-bound and unbound protein, and ∆G SA is the difference in the energy of surface area for the teicoplanin-protein complex and the sum of the surface area energies for the ligand and un-complexed protein.

Molecular Dynamics Simulation
The best ranked conformation of teicoplanin in complex with N-protein-N-domain furnished by MM-GBSA experiments was further assessed for their thermodynamic behavior and stability by employing molecular dynamics (MD) simulation. Desmond 6.1 program integrated with Maestro (Schrödinger, Inc. LLC, New York, NY, USA) was used for MD simulation studies [30,31]. In addition, apo form of the N-protein-N-domain was also simulated in order to understand the conformational changes upon teicoplanin binding. The teicoplanin-protein complex was placed into an orthorhombic box of 10 × 10 × 10 Å size, filled with 7210 water molecules by means of simple point charge (SPC) model implemented in system setup protocol. An OPLS3 force field was applied for the MD computations [32]. The system was neutralized using 20 and 26 Na + and Cl − ions, respectively, with a salt concentration of 0.15 M that represents the physiological concentration of monovalent ions. An isothermal-isobaric (NPT) ensemble was utilized with temper-ature and pressure attuned to 300 K and 1.01325 bar, respectively. A simulation time of 100 ns was adjusted, whereas trajectories were saved at every 100 ps. A cut-off radius of 9.0 Å was used for short-range van der Waals and Coulomb interactions. Nose-Hoover thermostat [33] and Martyna-Tobias-Klein [34] methods were employed for maintaining the system temperature and pressure, respectively. In order to integrate the equations of motion, an RESPA integrator was used with an inner time step of 2.0 fs for bonded as well as non-bonded interactions within the short-range cut-off [35]. The system was minimized and equilibrated with the default protocols of the Desmond. Simulation event analysis, simulation quality analysis, and simulation interaction diagram protocols of the Desmond package was exercised to analyze the trajectory files.

Post-Simulation MM-GBSA Analysis
Post-simulation MM-GBSA analysis was performed by using the thermal_MMGBSA.py script of the Prime/Desmond module of the Schrodinger suite 2020-3 [25,26]. From each MD trajectory, every tenth frame was extracted from the last 50 ns of simulated trajectory, averaging over 50 frames, for binding free energy calculations of teicoplanin in complex with N-protein-N-domain. The Prime MM-GBSA method uses the rule of additivity wherein total binding free energy (kcal/mol) represents a summation of individual energy modules like coulombic, covalent, hydrogen bond, van der Waals, self-contact, lipophilic, solvation, and π-π stackings of ligand and protein [36].

Validation of Docking Protocol
The co-crystallized ligand of main protease, N3, was redocked into the inhibitor binding cavity of the enzyme and docked conformation of the ligand was compared with that of the crystal structure in order to validate the docking protocol implemented in AutoDock Vina. The root-mean-square deviation (RMSD) was less than 2 Å between the docked and X-ray crystal structure conformations of N3, which authenticated the scoring functions implemented in the docking program (data not shown). According to the reported protocols, a successful docking computation should furnish the RMSD of ≤2.0 Å [10,14]. Therefore, AutoDock Vina used in this study was deemed reliable for studying teicoplanin interaction with potential SARS-CoV-2 targets.

Molecular Docking of Teicoplanin with Potential Targets of SARS-CoV-2
Molecular docking is a computational method of accelerating the early stages of drug discovery through unravelling the precise mechanism of intermolecular interactions of potential drug candidates by comparing their energetic compatibility and molecular shape with the target proteins [11,37]. With the aim of inspecting the interaction of teicoplanin with putative drug targets of SARS-CoV-2, docking studies was performed on the receptors/targets constituting both non-structural and structural proteins with prospective to be used in the drug discovery field [38,39]. Docking-predicted binding energies of teicoplanin with twenty-five COVID-19 drug targets are given in Table 1, whereas comprehensive intermolecular interactions have been listed in Table S1 in the Supplementary Information.

Prime MM-GBSA Calculations of Docked Complexes
Though several docking experiments are routinely employed to emphasize the binding mode and the affinity of a ligand relative to the protein, lack of accurate scoring function is one of the drawbacks of these algorithms. Therefore, MM-GBSA computations are routinely employed as a post-docking analysis in order to avoid false negatives as well as false positives [40,41]. The top ten conformations of protein-teicoplanin-docked complexes of each target were subjected to molecular mechanics/generalized Born surface area (MM-GBSA) computation for estimation of more accurate free binding energies. Table 2 enlists the MM-GBSA computed energies of best conformation of teicoplanin in complex with each target. Several other energy components such as Coulomb energy, hydrogen-bonding correction, lipophilic energy, generalized Born electrostatic solvation energy, Van der Waals energy, and ligand strain energy were also computed. Moreover, an elaborated MM-GBSA calculation of all poses has been listed in Table S2 in the Supplementary Information. As presented in Figure 4, nucleocapsid protein N-terminal domain ranked top among all studied targets exhibiting ∆G bind = −102.13 kcal/mol followed by main protease and spike protein-RBD with binding energies of −97.55 and −95.76 kcal/mol, respectively. However, MM-GBSA binding energy of teicoplanin-main protease complex obtained from AutoDock 4.2 has been recently reported as -68 kcal/mol. Interestingly, RdRp-RTP site was rendered as least promising for teicoplanin in this study, demonstrating 75.8 kcal/mol binding energy. The most promising target also hosted ample non-bond interactions with teicoplanin within the active site of N-domain of the N-protein. Moreover, appreciable intermolecular forces were established with all other targets. It is evident that the N-domain of the N-protein hosts nucleic acid substrates during processing of the viral genome into a ribonucleoprotein unit. Therefore, targeting this protein by teicoplanin can interrupt crucial steps of the SARS-CoV-2 life cycle such as transcription and translation processes [42]. Several proteins such as Nsp12, RdRp-RNA site, and papain-like protease were ranked as moderately favorable targets, whereas variable affinity was noted against the rest of the targets. The chemical structure of teicoplanin represents a unique blend of polar groups in the form of several amino acid fragments, carbohydrate cores such as N-acetylglucosamine, β-D-glucosamine and mannose, and non-polar fragments such as fatty acyl side-chain and numerous aromatic rings. These moieties enabled teicoplanin to interact with diverse drug targets of COVID-19 which warrants further experimental validation and may be useful in drug design based on teicoplanin template. Intermolecular interactions observed after MM-GBSA computation of all studied targets are presented in Figures S1-S4 and Table S1 in the Supplementary Information.   Table 2 for target names against respective numbers). The analysis is based on molecular mechanics-generalized Born surface area (MM-GBSA) computed binding energy (ΔGBind) in kcal/mol.

Molecular Dynamics Simulation Studies
Dynamic and thermodynamics parameters of living systems under specific conditions of physiological environments can be estimated by the application of molecular dynamics (MD) simulation, a widely employed computer-aided drug design technique [11,13,14,43]. Therefore, the best docked pose of teicoplanin in complex with SARS-CoV-  Table 2 for target names against respective numbers). The analysis is based on molecular mechanics-generalized Born surface area (MM-GBSA) computed binding energy (∆G Bind ) in kcal/mol.

Molecular Dynamics Simulation Studies
Dynamic and thermodynamics parameters of living systems under specific conditions of physiological environments can be estimated by the application of molecular dynamics (MD) simulation, a widely employed computer-aided drug design technique [11,13,14,43]. Therefore, the best docked pose of teicoplanin in complex with SARS-CoV-2 nucleocapsid protein N-terminal domain was subjected to MD simulation study in order to investigate the stability of the ligand-protein complex as well as main intermolecular interactions during the simulated trajectory. MM-GBSA-optimized pose bearing minimum binding energy was selected for MD simulation study. Desmond software was employed for the MD simulation of 100 ns in explicit solvent system. The resulting trajectories of the simulated complex was inspected for different standard simulation parameters such as backbone RMSDs for alpha-carbons, side chains, and heavy atoms. In addition, the rootmean-square fluctuations (RMSFs) of individual amino acid residue and intermolecular interactions involved were also evaluated. Figure 5A shows the RMSD of Cα atoms of apo protein whereas the complex of teicoplanin bound to N-protein-N-domain is presented as Figure 5B. The analysis of RMSD of teicoplanin-bound protein indicates that the simulated system fluctuated initially during 0-22 ns but acquired stability during rest of the simulated time. However, fluctuations can be clearly seen in the RMSD plot of the apo protein alone ( Figure 5A). Therefore, it seems logical to infer that upon teicoplanin binding, the protein maintains its conformational stability when simulated for 100 ns period. The conformational changes observed during simulated path has been visualized by the snapshots taken at the intervals of every 20 ns and presented in Figure 6.
The MM-GBSA method exploits molecular mechanics, the generalized Born solvation models and a solvent accessibility approach to estimate the free energies of binding based on snapshots obtained from MD simulations [40,44]. MM-GBSA computation was performed for the last 50 snapshots of 100 ns simulation, which estimated an average binding energy of −62.52 ± 12.22 kcal/mol.
Simulation interactions diagrams presented in Figures 8 and 9 during entire simulation time signifies a comprehensive intermolecular interaction profile of teicoplanin with nucleocapsid protein N-domain. The modus of interaction pattern of teicoplanin illustrates that only fewer docking-predicted main contacts were retained upon MD simulation time of 100 ns. Figure 8 presents intermolecular interactions observed in first and last frames of the teicoplanin bound to nucleocapsid protein N-domain. Residues participating in both hydrogen bonding and water bridges include Asn7, Asn8, Thr9, Ser11, Thr51, Arg53, Arg55, Tyr71, Pro111, and Asn114. In addition, Thr51, Arg53, Arg55, and Arg67 also contributed in ionic bond interaction. However, amino acid residues Ala10, Ala50, Tyr69, and Ala112 were noted to be important for hydrophobic interaction.    The local conformational alterations along SARS-CoV-2 nucleocapsid protein N-domain were investigated by analyzing the RMSF during simulation time. Terminal residues fluctuated most in both bound and unbound states. Maximum fluctuation was noted with Gly1 as 13.29 Å which was declined upon teicoplanin binding at 7.00 Å (Figure 7). Key residues involved in hydrogen bonding showed fluctuations at 3.74, 2.97, 1.76, 1.54, 4.12, 0.96, and 2.68 Å by Asn7, Asn8, Thr9, Ser11, Arg55, Tyr69, and Pro111, respectively. Simulation interactions diagrams presented in Figures 8 and 9 during entire simulation time signifies a comprehensive intermolecular interaction profile of teicoplanin with nucleocapsid protein N-domain. The modus of interaction pattern of teicoplanin illustrates that only fewer docking-predicted main contacts were retained upon MD simulation time of 100 ns. Figure 8 presents intermolecular interactions observed in first and last frames of the teicoplanin bound to nucleocapsid protein N-domain. Residues participating in both hydrogen bonding and water bridges include Asn7, Asn8, Thr9, Ser11, Thr51, Arg53, Arg55, Tyr71, Pro111, and Asn114. In addition, Thr51, Arg53, Arg55, and Arg67 also contributed in ionic bond interaction. However, amino acid residues Ala10, Ala50, Tyr69, and Ala112 were noted to be important for hydrophobic interaction.

Conclusions
By using computer-aided drug design techniques, the current study explained the intermolecular interaction of the antibacterial drug, teicoplanin, with potential drug targets of SARS-CoV-2. Molecular docking studies employing AutoDock Vina highlighted the importance of hydrophilic and hydrophobic interactions in supporting the teicoplanin Figure 9. Protein interactions with teicoplanin, monitored throughout the simulation trajectory. These interactions are clustered by type and summarized in bar diagram including hydrogen bonds (H-bonds), hydrophobic, ionic, and water bridges. Ala: alanine; Arg: arginine; Asn: asparagine; Asp: aspartic acid; Cys: cysteine; Gln: glutamine; Glu: glutamic acid; Gly: glycine; His: histidine; Leu: leucine; Met: methionine; Phe: phenylalanine; Pro: proline; Ser: serine; Thr: threonine; Tyr: tyrosine; Val: valine. The number of corresponding amino acid residue is according to the crystallographic information obtained from the protein data bank.

Conclusions
By using computer-aided drug design techniques, the current study explained the intermolecular interaction of the antibacterial drug, teicoplanin, with potential drug targets of SARS-CoV-2. Molecular docking studies employing AutoDock Vina highlighted the importance of hydrophilic and hydrophobic interactions in supporting the teicoplanin molecule inside the nucleocapsid protein N-domain binding cavity. MM-GBSA analysis and molecular dynamics simulation results not only reinforce the credibility of the docking results, but also authenticate the stability of the simulated system. This study is expected to assist lead optimization and design of COVID-19 drugs based on molecular skeleton of teicoplanin. However, the requirement of additional experimental and clinical validation is an obvious limitation of this study. Nevertheless, findings of this study justify further exploration of teicoplanin repurposing.  Table S1: Intermolecular interactions observed in MM/GBSA optimized complexes of teicoplanin and 25 potential COVID-19 targets (AutoDock Vina 1.1.2 was used for molecular docking); Table S2. Computation of MM/GBSA energies of top ten docked poses of teicoplanin in complex with 25 potential COVID-19 targets.

Funding:
The author gratefully acknowledges the Deanship of Scientific Research, Qassim University, Saudi Arabia for funding the publication of this research.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.