The conservation of FtsZ protein was experimentally determined by performing BLAST analysis of whole UNIPROT database sequences. The default parameters like E-Threshold as 10, and the Auto-BLOSSUM 62 matrix were used. A pairwise sequence alignment of the two FtsZ sequences from S. epidermidis and S. aureus was performed using Clustal omega to identify the regions of similarity which indicates the structural, functional, and evolutionary relationship between the two sequences.

The 3D structure of the S. epidermidis FtsZ protein was predicted using the Swiss-Model server because the detailed structural information of the crystallized structure is unavailable in PDB as the literature report of its PDB entry (4M8I) is not published. The protein sequence was retrieved from the Uniprot database (Uniprot ID: Q5HQ06). The template was selected based on the parameters observed in the BLAST findings by mainly focusing on the sequence similarity, resolution, and experimental technique used to determine the structure. (Choudhary et al., 2020). The target protein’s predicted 3D structure was validated using the Ramachandran (RC) Plot. The 3D protein structure was visualized using Maestro 13.1. After that, the built model was pre-processed using the protein preparation wizard of Schrodinger-suite 2022, which refine the protein structure for docking by setting bond orders, filling in missing loops, adding hydrogens, and deleting water molecules that are more than 3Å distances from the protein (Sahayarayan et al., 2021). After the H-bond assignment, H-bonds were optimized using PROPKA tool. The binding sites of the FtsZ protein were anticipated using the Site Map tool of the Schrodinger-suite 2022 which predicted best five sites for a given protein entry.

The receptor grid was generated using the “Receptor grid generation” panel from the glide module of Schrodinger software by preserving the grid’s default settings and size. A receptor grid was generated for both the nucleotide and inter-domain binding sites.

The literature review obtained about 138 reported natural, semi-synthetic, and synthetic FtsZ inhibitors of various other pathogenic bacterial species. The experimentally proved, 138 natural, semi-synthetic, and synthetic FtsZ inhibitors of various other pathogenic bacterial species were obtained (Tripathy and Sahu, 2019). Among these ligands, the chemical structures of a few were obtained from PubChem, and the remainder were sketched using maestro 13.1’s 2D Sketcher. To prepare ligands, these structures were loaded into Schrodinger Workspace. The Ligprep module of Schrodinger-suite 2022 was used to prepare the ligands by optimizing their geometrical features and generating ionization states for the compounds to achieve the necessary pH of 7.0 ± 2.0.

Using the “ligand docking” panel of Schrodinger software, the prepared ligands were docked against both the nucleotide and inter-domain binding sites of the FtsZ protein. In order to generate hits from the ligand dataset, ligand docking was first carried out using the high throughput virtual screening (HTVS) approach with the precision mode set to HTVS, followed by standard precision (SP), and extra precision (XP). While taking the docking score into account, the Epik state penalties of the ligands were modified. (Kapusta et al., 2021). The docking validation is accomplished by redocking co-crystal ligands to their specific binding site of the receptor protein (El-Far et al., 2020). Dacomitinib (S2727), used as a control drug to assess the binding affinities and free energy of the proposed inhibitors, is a promising FtsZ inhibitor that Du et al. identified using in vivo and in vitro bioassays. while PC190723 (Elsen et al., 2012) served as the inter-domain binding site control drug. Instead of suppressing FtsZ filament assembly and condensation, PC190723 (difluoro-benzamide derivative) induces it (Andreu et al., 2010), causing FtsZ to assemble into delocalized cellular foci as opposed to the Z-ring (Adams et al., 2011).

The binding free energies of the top 10 docked complexes (nucleotide-binding site and inter-domain binding site) were determined using the Prime MM/GBSA module of Schrodinger-suite (Muthumanickam et al., 2022; Ramachandran et al., 2022). The equation used for calculating free energy is as follows Δ G bind = G complex − ( G protein + G ligand ) The G_complex indicates complex energy, G_protein indicates receptor energy, and G_ligand indicates the unbound ligand energy.

The Qikprop module of Schrodinger’s suite 2022 was used to forecast the pharmacokinetic features, also known as absorption, distribution, metabolism, and excretion, of the top 10 hit compounds (Guo et al., 2016). This module provides information on the drug-like characteristics of the given ligands, such as molecular weight (mol MW), the number of hydrogen bond acceptors and donors (accptHB), solubility (QPlogS), Octanol/water Partition coefficient (QPlogPo/w), the percentage of oral absorption, apparent Caco-2 cell permeability (QPPCaco), and the brain/blood partition coefficient (QPlogBB).

ProTox—II was used to calculate the toxicity of the hit compounds. This programme delivers data based on chemical similarities, fingerprint propensities, etc. The ProTox-II employs machine learning models to determine the toxicity class, LD 50 values, organ toxicity, and toxicity endpoints, including hepatotoxicity, carcinogenicity, mutagenicity, immunogenicity, etc. Using a system that is universally accepted, toxicity classifications for chemicals are defined (GHS) (Banerjee et al., 2018). LD50 values are given in [mg/kg]:

Class I: If ingested, deadly (LD50 ≤ 5)

The molecular dynamics (MD) simulation of the top docking scored protein-ligand complex was performed using the desmond module of Schrodinger suite 2022 for 100 ns. The simulation system was built using a system builder task where a simple point charge water (SPC) model was used as a solvent system, and periodic boundary condition (PBC) was applied using an orthorhombic boundary box. The system was neutralized by adding counter ions like Na⁺ and Cl⁻. The OPLS4 force field was used for energy minimization of molecular dynamics system. For the simulation of the build system it was loaded on to workspace, molecular dynamics simulation task was performed for 100 ns with a trajectory recording interval of 100 ps while the system was equilibrated with NPT ensemble system with pressure 1.01 bar and temperature 300 k, respectively. Finally, RMSD, RMSF, protein-ligand interaction and ligand properties were used for analyzing molecular dynamics simulation results (Kaushik et al., 2018).

Using the Prime MM/GBSA module of the Schrodinger, the post MM/GBSA analysis was carried out for the complexes at the 100 ns time frame of the molecular dynamics simulation to determine the binding free energies.

With the aid of BLAST, the entire Uniprot database was aligned using the reference sequence of the FtsZ protein from S. epidermidis (Uniprot ID: Q5HQ06). This produced 250 hits, with the sequence identities ranging from 100% to 65.9%, indicating that they are the best matches or homologs present in different bacterial species such as Staphylococcus species, Macrococcus caseolyticus, Abyssicoccus albus, Bacillus sp. (Supplementary Figure S5) The FtsZ protein of S. epidermidis has shown a high sequence similarity about 92.39% identity with that of its kin i.e., S. aureus (Figure 3).

Using the Swiss model server, FtsZ protein homology modeling was carried out by choosing a template, PDB ID: 4M8I (Organism: Staphylococcus epidermidis RP62A) with 100% sequence identity and a Global Model Quality Estimation (GMQE) of 0.81%. The target protein’s 3D model was validated by the Ramachandran plot using the Swiss Model’s structure evaluation. (Figure 4). The 3D protein structure is visualized using Maestro 13.1). The active site for the nucleotide binding region was determined by using a co-crystallized ligand present in the template that we used for homology modeling followed by supplying the X, Y, and Z coordinates as -18.5, −9.81, and 19.98 Å, respectively, and retaining the other parameters at their default values, the receptor grid was generated at the nucleotide-binding site of FtsZ. The inter-domain binding site was predicted using the sitemap tool of Schrodinger-suite 2022. The top-ranked site showed a site score of interdomain binding site, volume = 116.620 Å³, hydrophilic score = 0.781, and hydrophobic score = 0.746. This predicted site contained Gln192, Gln195, Gly196, Asp199, Leu200, Val203, Leu209, Ile228, Leu261, Asn263, Ile264, Thr265, Val297, Asn299, Leu302, Val307, Thr309, and Ile311 residues. Similarly, the site coordinates X = −0.31, Y = 13.26, and Z = 24.61Å were provided to generate a receptor grid at the inter-domain binding site with the default settings for the remaining parameters.

Molecular docking was performed sequentially using the Schrodinger’s Glide module in three modes: HTVS, SP, and XP. It is crucial to know the binding affinity since docking was used to explore molecules that might inhibit the FtsZ protein. This research article discusses the outcomes of glide-XP docking since the extra precision (XP) mode generates precise results and uses the chemscore scoring program to evaluate the docked complex. The top 10 ligand structures are displayed in Table 1, along with their docking data and receptor-ligand interactions.

Out of 138 screened molecules, Compound A (2R,3S,4R,5R)-5-(2-amino-8-chloro-6-oxo-1,6-dihydro-9H-purin-9-yl)-3,4-dihydroxytetrahydrofuran-2-yl)methyl triphosphate) i. e, chloro derivative of GTP had the highest docking score of −9.549 kcal/mol compared to Compound I ((2R,3S,4S,5R)-5-(2-amino-6-oxo-1,6-dihydro-9H-purin-9-yl)-3,4-dihydroxytetrahydrofuran-2-yl)methyl diphosphate i.e., GDP molecule, a co-crystallized ligand which was found in the template protein, with a docking score of −8.886 kcal/mol. Following analysis of the 2D interaction diagram of the top four ligand-protein complex, it was found that while the ring structure of the GTP derivative produced a Pi-Pi stacking interaction with the aromatic amino acid Phenylalanine (Phe) 183 of active site, phosphate groups interacted with the majority of amino acids by forging hydrogen bonds Figure 5. (Figure 5A). Compound B produced hydrogen bonds and pi-cation interactions with the amino acids Gly104, Gly108, Thr109, Gly110, Asn166, Asp187, and Arg143, respectively (Figure 5B). It was observed that the compound C have shown similar interactions as that of compound A. (Figure 5C). The Compound D has interacted with Gly21, Arg29, Ala71, Gly72, Ala73, Met105, Gly108, Thr109, Gly110, Thr133, and Arg143 residues via hydrogen bonding (Figure 5D). whereas the compound E formed hydrogen bonding interaction with Gly21, Gly22, Ala71, Ala73, Met105, Gly108, Thr109, and Arg143 residues (Figure 5E). in the nucleotide binding site. The hydrogen bond interactions with Gly108 was found common in both GTP and non GTP derivatives which can be assumed as a crucial interaction for inhibiting FtsZ activity. The control drug i.e., Dacomitinib has interacted with Gly22, Thr109, Arg143 via hydrogen bonding and Salt Bridge with Glu139 of FtsZ’s nucleotide binding site. The docking results were validated by redocking the co-crystallized ligand (GDP) in the template protein (PDB—4M8I) to its nucleotide binding site. The co-crystal ligands’ original and docked confirmations were compared, and the computed root means square deviation (RMSD) between them was less than 2 Å i.e. 0.4302 Å (Figure 6).

Small molecules bind to the region between the N and C-terminal guanosine triphosphate (GTP) binding subdomains i.e., the inter-domain binding site of FtsZ, reducing its capacity to function allosterically, which eventually prevents bacterial division. However, the lack of adequate chemical tools to develop a binding screen against this region has hindered the search for FtsZ antibacterial inhibitors. (Huecas et al., 2021).

At the FtsZ protein’s inter-domain binding site, ligand docking was conducted using the ligand molecules taken from the literature. The top 10 outputs’ binding affinities range was −7.243 to −4.290 kcal/mol (Table 2). According to the findings of this study, the ligands interact with the active residues Asp199, Glu301, Asn263, Thr265, Glu192, Thr309, and Gly196 through interactions that entail both salt bridges and hydrogen bonds.

The two-dimensional interaction diagram of the ligand-protein complex at this binding site is shown in (Figure 7). While Glu301 formed hydrogen connections with Compounds 1 and 2, Asp199 altered the salt bridge (Figures 7A,B). Compound 3 (chlorogenic acid) formed hydrogen bonds with Asn263, Thr265, and Asp199 (Figure 7C). Compound 4 has shown only hydrogen bond interaction with Asp199 (Figure 7D). Compound 5 has interacted with Arg191, Glu192, Asn263 residues via hydrogen bonding (Figure 7E). The interaction with Asp199 via hydrogen bonding or through salt bridge was found to be common among all the hits. This indicates that the interaction with Asp199 may be essential for the compounds to exhibit FtsZ inhibitory activity.

The binding free energies of the top five docked complexes were calculated using the MM/GBSA approach. The results revealed that the binding free energy of FtsZ’s nucleotide-binding site with compounds A, B, C, D and E were found to be −41.92, −58.09, −65.68, −34.80, and −26.90 kcal/mol, respectively (Table 3). According to the total binding free energies, the hit molecules may form a potent interaction within the binding region of the selected target, thereby inhibiting enzyme activity.

The binding free energies of the top five docked complexes were calculated using the MM/GBSA approach. The results revealed that the binding free energy of FtsZ’s nucleotide-binding site with compounds 1, 2, 3, 4, and 5 was found to be −49.14, −48.19, −19.67, −28.92, and −40.51 kcal/mol, respectively (Table 4). These hits may generate a strong interaction within the binding region of the chosen target, according to the total binding free energies, which would impede enzyme activity.

ADME parameters were predicted for the top 5 hits of the nucleotide-binding site (Table 5). It was noticed that although most of the successful compounds were predominantly nucleotide analogues, which have shown excellent results, their expected oral absorption percentage in humans was low. The top 10 hit molecules were anticipated to be non-toxic compounds, as determined by Protox II (Table 6).

Following docking to the inter-domain binding site of FtsZ, the top 5 hits were predicted for their ADME properties (Table 7). The majority of the compounds have shown encouraging findings, and it was found that the projected values fall within the acceptable range. Although none of the compounds shows organ toxicity, Protox II predicted toxicity endpoints like immunotoxicity. (Table 8).

To investigate the real-time dynamics and conformational stability of a protein upon binding to a specific ligand, 100 ns of MD simulations were performed on the five best-docked compounds, (A, B, C, D, and E) at nucleotide binding site and (compound 1, 2, 3, 4, and 5) at inter-domain binding site. Our analysis of simulation interaction diagrams (SIDs) for the 100-ns SPC water model-based simulations provided a better understanding of Protein RMSD, Ligand RMSD, Protein RMSF, Ligand RMSF, Protein-Ligand contacts, and ligand characteristics were analyzed, and a simulation interaction diagram was constructed to assess the stabilities of the protein-ligand complexes.

The RMSD is generally used to calculate protein backbone (Cα, C, and N) deviation during the 100 ns simulation period. The FtsZ bound to compound A (Chloro derivative of GTP) showed a protein RMSD of 1.6–3.2 Å and ligand RMSD of 1.1–2.0 Å. Convergence was observed between the two plots around the initial 0–35 ns, which indicates the FtsZ-Chloro derivative of GTP complex’s stability. Additionally, for ligand, slight fluctuation in RMSD, was observed between 35–45 ns of trajectory; however, it remained consistent after that with RMSD of 0.6–1.8 Å Figure 8 (Figure 8A). In case of Compound B (naphthalene-1,3-diyl bis(3,4,5-trihydroxybenzoate)-FtsZ complex, the protein and ligand RMSD plots slightly fluctuated during the initial 0–20 ns and remained converged throughout the simulative period with RMSD of 1.8–5.6 Å which can be inferred as a stable complex (Figure 8B). The compound C-FtsZ complex MD trajectory of 100 ns indicates that the complex tends to be stable during simulation with a RMSD of less than 3 Å. This indicates that the FtsZ- Guanosine triphosphate (GTP) complex did not lead to much conformational changes in the dynamics environment (Figure 8C). The ligand RMSD of compound D bound to FtsZ was found to be in the range of 2.0–2.4 Å, and the protein’s RMSD lay in the range of 1.8–3.2 Å. Initially, the slight fluctuation was observed in the trajectory at 0–15 ns time frame followed by the convergence of both the RMSDs till 60 ns, a slight change in order was observed at 70–100 ns, which falls in the acceptable range (Figure 8D). This ultimately indicates the FtsZ-morpholine derivative of GTP complex’s stability. The FtsZ protein’s RMSD was 1.8–7.8 Å, whereas compound E’s RMSD was 1.1–8.0 Å (methylpiperazine derivative of GTP). Even though the RMSD plots of the ligand and the protein converged for a certain period, after 70 ns, the complex was entirely out of the binding pocket when the simulation trajectory was observed. As the ligand was found to show significant variations in their RMSD, compound E is considered an unstable molecule (Figure 8E).

The RMSD of the protein in the presence of compound 1 [1-(2-((amino(iminio)methyl)amino)ethyl)-3-(3-(tert-butyl)phenyl)-6,7-dimethoxyisoquinolin-2-ium] was found to be 1.8–10.5 Å whereas the ligand RMSD was in the range of 0.8–12 Å (Figure 9 (Figure 9A). Variation in both protein and ligand RMSDs was recorded until the end of the simulation, indicating that the ligand diffused away from the binding and complex is unstable. The compound 2 RMSD plot converged well with the FtsZ protein backbone for the 100 ns simulation period with a RMSD of 2.3–4.2 Å indicating ligand’s tight binding to the protein’s binding pocket and hence a stable complex (Figure 9B). In the case of compound 3 [chlorogenic acid], the RMSDs of the FtsZ protein and ligand were found fluctuating for the initial 0–10ns, but were found to show consistent stability with a ligand RMSD of less than 3.6–4.8 Å and protein RMSD was in the range of 2.9–6.2 Å (Figure 9C). The protein RMSD was found to be in the range of <3 Å in the case of both compounds 4 [(Z)-2-(2-(5-fluoro-2-methyl-1-(4-(methylthio)benzylidene)-1H-inden-3-yl)acetamido)-N,N-dimethylethan-1-aminium] and compound 5 [Curcumin], but the ligand RMSDs was found to be 0.9–10.5 Å and 0.5–10 Å respectively. Due to the large deviation among the protein and ligand RMSDs, the compound 4-FtsZ complex was found to be unstable throughout the 100ns simulation (Figure 9D) and when the simulation trajectory was observed, the ligand protruded away after 20 ns from the binding pocket indicating poor stability with FtsZ. The compound 5-FtsZ complex was unstable during the 100 ns simulation due to a substantial discrepancy between the protein and ligand RMSDs (Figure 9E).

Based on the RMSDs of complexes, the stable complexes were considered for further analysis of simulation parameters like RMSF of protein and ligand, protein-ligand contacts and changes in the ligand properties during simulation.

Each protein residue’s flexibility and mobility are represented by the RMSF value. Greater RMSF values suggest more flexibility throughout the MD simulation, while a lower RMSF value affects the system’s stability. The protein’s RMSF in the presence of compound A was found to be 0.5–2.0 Å Figure 10 (Figure 10A). Similarly, compound A did not show any major fluctuations in their RMSF, while element number 10 of the ligand (involved in making hydrogen bond interactions with the active site) exhibited small fluctuation, which is under the acceptable range, retains the stability of Compound A-FtsZ complex (Figure 10B). In the presence of compound B, the backbone residues demonstrated a substantially lower than permitted amount of fluctuation and was determined to be 0.8–2.4 Å Figure 11 (Figure 11A). The fluctuation was only observed at ligand’s element number 5, which is not engaged in making any contact with the active site (Figure 11B). According to the protein RMSF assessment results, the residues did not exhibit any observable flexibility in the presence of compound C, which ranged from 0.5 to 1.8 Å Figure 12 (Figure 12A). The ligand-protein complex may not result in any structural variation as the ligand’s RMSF was within the acceptable range (Figure 12B). In the presence of compound D, the protein’s RMSF was found to be 0.6 to 1.8 Å was the observed RMSF of protein which falls under the acceptable range and indicates the protein’s stability Figure 13 (Figure 13A). Variations in the ligand’s elements 36 and 37, which do not interact with the active site, may not lead to structural variability in the ligand-protein complex (Figure 13B).

The protein’s RMSF showed fewer fluctuations i.e., 1.2 and 2.4 Å in the presence of compound 2, indicating protein stability (Supplementary Figure S1A). There is slight fluctuation observed in the ligand’s RMSF among the elements 18 and 19, which are involved in forming hydrogen bonds with the active site in the ligand (Supplementary Figure S1B). This fluctuation may not affect the structural variability of the Compound 2-FtsZ complex as the observed RMSF was within the acceptable range. The RMSF of protein backbone residues in the inter-domain binding site vary from 1.9–2.0 Å, which is under the acceptable range (Supplementary Figure S2A), and the RMSF for Compound 3 was found to be 1.5–3.0 Å where the element numbers 5 and 6 of the ligand, which are not engaged in forming contacts with the active site, have no impact on the structural variation of the Compound 3-FtsZ complex (Supplementary Figure S2B).

During the simulation phase, the atomic-level interaction information is critical for predicting the binding affinity of the protein and ligand. During the 100 ns simulation, intermolecular interactions between protein and ligand molecules, such as hydrogen bonds, ionic interactions, hydrophobic contact, and the salt bridge, were thoroughly investigated for binding analysis. The compound A was in contact with 34 amino acid residues in the nucleotide-binding domain, forming crucial interactions like hydrogen bonds, water bridges, ionic bonds and hydrophobic bonds (Figure 10C). 32 amino acid residues of the nucleotide-binding site came into contact with compound B during the simulation, creating critical interactions like hydrogen bonds, water bridges, ionic bonds, and hydrophobic bonds, most of which are hydrogen bond interactions (Figure 11C). Throughout the simulation period, 30 amino acid residues in the nucleotide-binding site came into contact with compound C, forming crucial interactions such as hydrogen bonds, water bridges, and hydrophobic bonds (Figure 12C). 40 amino acid residues in the nucleotide-binding domain interacted with chemical D, forming important interactions such as hydrogen bonds, water bridges, ionic bonds, and hydrophobic bonds (Figure 13C).

During the course of a 100 ns simulation, 16 amino acid residues in the inter-domain binding site made contact with compound 2, forming significant interactions such as hydrogen bonds, water bridges, ionic bonds, and hydrophobic bonds, the majority of which are hydrogen bond interactions (Supplementary Figure S1C). Compound 3 generated significant interactions with 25 amino acid residues of FtsZ, including hydrogen bonds, water bridges, ionic bonds, and hydrophobic bonds, among others, the water bridges and hydrogen bonds making most of these interactions (Supplementary Figure S2C).

To assess the stability of the lead molecules in the nucleotide and inter-domain binding sites of the FtsZ, the five molecular characteristics of ligand (ligand RMSD, the radius of gyration [rGyr], Molecular surface area [MolSA], solvent accessible surface area [SASA], and polar surface area [PSA]) were also investigated.

As stated in Table 9, the molecular dynamics simulation was generated with a 100 ns frame time interval for the purpose of analyzing the binding free energy calculation following MM/GBSA. The findings revealed that, on average, complexes (Compound A, B, C, and D) at nucleotide binding sites had better binding free energies at post simulation period (−97.15, −57.11, −73.07, and −42.02 kcal/mol) compared to pre-MM/GBSA complexes, but the complexes (Compound 2 and 3) at inter-domain binding sites had comparable binding free energies (−45.49 and −45.38 kcal/mol) compared to pre-MM/GBSA free energies. As there is less variation among the pre and post MM/GBSA energies of complexes at the inter-domain site, they can be considered as strong binders (Table 10). It was observed that for the binding of compounds at both binding sites, the MM/GBSA dG bind coulomb, and MM/GBSA dG bind solv GB energies majorly contributed to the binding energies of the complexes.