Background. Aliarcobacter butzleri is a Gram-negative, curved or spiral-shaped, microaerophilic bacterium and causes human infections, specifically diarrhea, fever, and sepsis. The research objective of this study was to employ computer-aided drug design techniques to identify potential natural product inhibitors of a vital enzyme in this bacterium. The pyrimidine biosynthesis pathway in its core genome fraction is crucial for its survival and presents a potential target for novel therapeutics. Hence, novel small molecule inhibitors were identified (from traditional Chinese medicinal (TCM) compound library) against it, which may be used for possible curbing of infection by A. butzleri. Methods. A comprehensive subtractive genomics approach was utilized to identify a key enzyme (orotidine--phosphate decarboxylase) cluster conserved in the core genome fraction of A. butzleri. It was selected for inhibitor screening due to its vital role in pyrimidine biosynthesis. TCM library ( compounds) was screened against it using pharmacophore model based on orotidylic acid (control), and the obtained lead-like molecules were subjected to structural docking using AutoDock Vina. The top-scoring compounds, ZINC70454134, ZINC85632684, and ZINC85632721, underwent further scrutiny via a combination of physiological-based pharmacokinetics, toxicity assessment, and atomic-scale dynamics simulations (100 ns). Results. Among the screened compounds, ZINC70454134 displayed the most favorable characteristics in terms of binding, stability, absorption, and safety parameters. Overall, traditional Chinese medicine (TCM) compounds exhibited high bioavailability, but in diseased states (cirrhosis, renal impairment, and steatosis), there was a significant decrease in absorption, Cmax, and AUC of the compounds compared to the healthy state. Furthermore, MD simulation demonstrated that the ODCase-ZINC70454134 complex had a superior overall binding affinity, supported by PCA proportion of variance and eigenvalue rank analysis. These favorable characteristics underscore its potential as a promising drug candidate. Conclusion. The computer-aided drug design approach employed for this study helped expedite the discovery of antibacterial compounds against A. butzleri, offering a cost-effective and efficient approach to address infection by it. It is recommended that ZINC70454134 should be considered for further experimental analysis due to its indication as a potential therapeutic agent for combating A. butzleri infections. This study provides valuable insights into the molecular basis of biophysical inhibition of A. butzleri through TCM compounds.

1. Introduction

Aliarcobacter butzleri was first described in 1991 and is found in water, soil, and various food sources [1]. It is considered an emerging foodborne pathogen and has been associated with a wide range of gastrointestinal diseases [2]. The pathogenesis of A. butzleri is not fully understood, but it is thought to involve a range of virulence factors, including adhesins, invasins, and toxins [3]. The organism is known to colonize the intestinal tract and has been shown to adhere to and invade intestinal epithelial cells [4]. It produces a variety of toxins, including cytotoxins, enterotoxins, and hemolysins, which can cause damage to host cells and tissues [5]. It can survive in harsh environments, such as low pH, high salt concentrations, and low oxygen levels, which makes it difficult to control and prevent infection [6]. The organism is known to be resistant to a variety of antibiotics, including macrolides, fluoroquinolones, and tetracyclines [7]. The diagnosis of A. butzleri infections is usually based on the isolation of the bacteria from clinical samples and the identification by laboratory methods such as microscopy, culture, and molecular biology [8]. Diagnosis of A. butzleri infection can be challenging, as the organism is not easily detectable using standard microbiological methods. However, molecular methods, such as PCR and DNA sequencing, can be used to identify the organism in clinical and environmental samples [9].

Whole genome of A. butzleri has been sequenced from various sources, like cattle [10], milk [11], poultry [12], shellfish [13], and diseased humans [14]. Pan-genome on 32 strains of this bacterium has previously been attempted and shown to have a hypervariable accessory genome [15]. Additionally, virulence and immune response inciting genes have been identified. Pan-genomics, which is the study of the entire genetic repertoire of a group of organisms, has previously revealed the existence of conserved gene clusters that are present in multiple strains of a species and are important for drug design [16]. By targeting such gene clusters, a broad-spectrum inhibitor may be obtained that could be effective against multiple strains of a pathogenic bacterial specie, including strains that have developed resistance to existing treatments. Previously, efflux pumps in the core genome have been identified as drug target in Yersinia sp., using pan-genomics, and used for identifying inhibitors from traditional medicinal compounds [17]. In another study, L-lysine biosynthesis pathway protein was targeted using a subtractive core genomics approach and phytochemical inhibitors were identified against Orientia sp. [18].

Phytochemicals, which are naturally occurring compounds found in plants, have been found to have antibacterial properties and may be used as an alternative or alongside antibiotics for synergistic impact in some resistant isolates [19]. Phytochemicals may be better than antibiotics in curbing bacteria due to their broader spectrum of activity and lower risk of resistance development due to their multiple mechanisms of action [20]. Thus, bacteria may be less likely to develop resistance to them compared to antibiotics, which typically target a single bacterial pathway. Some phytochemicals have been found to have synergistic effects when combined with antibiotics, which can enhance the effectiveness of antibiotics. Antibiotics can have side effects such as gastrointestinal disturbances, allergic reactions, and drug interactions [21]. Phytochemicals from traditionally used medicinal plants are generally considered safer with fewer side effects, due to their innocuous usage by humans since long times [22].

Traditional Chinese medicine (TCM) has a long history of use in the treatment of infectious diseases, and many TCM compounds have been shown to have antibacterial activity [23]. In recent years, virtual screening methods have been used to identify TCM compounds that may be effective in inhibiting bacteria [24]. Some examples of TCM compounds that have been identified using virtual screening methods include antibiofilm agents [25], purine synthesis inhibition [26], and amino acid synthesis inhibition [27] in pathogenic bacteria. In pursuit, the present study attempted to identify drug-like compounds from TCM library against ODCase enzyme of A. butzleri. Dynamics simulation was conducted to assess stability of these compounds while ADMET profiling was done for toxicity assessment. PBPK modeling was also carried out for determining possible bioavailability and other parameters of importance, required for effectiveness as a drug.

2. Material and Methods

The phytochemical inhibitor prioritization approach consisted of several steps outlined below:

2.1. Pan-genomics

Pan-genome of A. butzleri was accessed from ProPan database (https://ngdc.cncb.ac.cn/propan/pananalysis/28197; accessed 16 Feb 2023). The database integrates information on pan-genomes (the collection of all genes present in a set of genomes from a given species or population) and gene families from over 2000 prokaryotic species [28]. The database allows users to browse the genomes of individual prokaryotic species and visualize the distribution of gene families and other genomic features. It is designed to provide a comprehensive overview of the evolution and diversity of prokaryotic genomes and to facilitate comparative genomics and functional annotation. Nucleotide diversity analysis and COG function categories were obtained for gene clusters for A. butzleri. In total, 86.7 MB gene file, with a total 9870 gene clusters from 42 strains (listed in Table 1), was obtained. They were differentiated into core, unique, and pan-fraction gene clusters. Core clusters () were subjected to in-house subtractive genomics pipeline in bash, according to parameters described in the previous studies [29, 30]. Initially, redundant sequences were removed by CD-Hit. Remaining sequences were BLAST against DEG and CEG databases for essentiality analysis. Translated product of common genes from these databases was then BLAST against human proteome and, later, DrugBank. Obtained hits were labelled as drug targets.

2.2. Target Selection and Inhibitor Screening

Among the core gene clusters identified as drug targets, one depicting ODCase (cluster 1186) from pyrimidine metabolism pathway was selected for further downstream analysis based on its novelty, utility, and importance for the bacterial survival. 3D structure of the protein was modeled using SWISS-model [31] and assessed using ERRAT at SAVES server (https://saves.mbi.ucla.edu/) and Ramachandran plot through PDBSum server (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/). Ligand and protein structure was prepared before docking by assigning bond orders, protonating, and adding missing atoms and energy minimization [32]. Orotidylic acid, natural binder of ODCase [33], was taken as control for docking. TCM library of 36,000 molecules was filtered for drug-like candidates using Lipinski’s druggability filter and 11,431 compounds fulfilling the criteria were subjected to pharmacophore generation, according to the parameters described previously [34]. Screening against the obtained 77 compounds was then carried out using AutoDock Vina. Fpocket was used for pocket calculation, and pocket parameters were 3-6 Å radius with druggability score between 0 and 1 and volume of 100-2000 Å3. Box offset of 12 Å and pocket having a druggability score 0.8 and volume 366 were used for analysis. Only the pose with best binding affinity was retained. Interactions were mapped, and PRODIGY (PROtein binDIng enerGY prediction) server (https://bianca.science.uu.nl/prodigy/lig; accessed 9 July 2023) was used for calculating (kcal/mol) of the complexes.

2.3. ADMET Profiling

The PKCSM tool (https://biosig.lab.uq.edu.au/pkcsm/; accessed on 18 Feb, 2023) was used to conduct ADMET analysis. This tool utilizes graph-based machine learning and consists of 14 regression models and 16 classification models [35]. These models predict various pharmacokinetic or toxicity properties and categorize outcomes into two classes, respectively. The ADMET descriptors evaluated in this analysis include caco-2 permeability, intestinal absorption in humans, BBB permeability, cytochrome inhibition, renal enzyme substrate, total clearance, Rat LD50, and Ames toxicity [36]. The SMILE string was used as input for the evaluation of these parameters, and the obtained results were compared for top compounds.

2.4. Pharmacokinetics

To simulate the physiological-based pharmacokinetic (PBPK) parameters of a 100 mg oral tablet of the compound, a compartmental model was developed using GastroPlus software (version 9.8.2, Simulation Plus, Inc., Lancaster, PA, USA). The simulation was conducted on human subjects in a fasting prandial state, with two sets of population taken into account. The first set consisted of 100 healthy individuals with an age range of 20-60 years, measure  h, and a pH of 7.2. The second set consisted of 100 cirrhotic, 100 renally impaired, and 100 people having steatosis. Additionally, 100 pregnant women with an age range of 20-50 years were put in this group. Physiological measurements for various parameters were accounted (Supplementary Figure 1), and to ensure consistency with the previous studies [24, 37], the particle radius, dissolution model, and density were kept constant. To determine the pKa values of the compound, the ADMET profiler version 10 was utilized. A fixed first pass was assumed for the liver, and a separate jejunal and paracellular permeability model was included. A persistent electrical potential gradient was assumed for the length of the intestinal tract to enhance accuracy. The simulation included the calculation of the compound’s bioavailability, absorption, and plasma concentration to aid in determining optimal dosing. Mean of the values obtained for 100 runs of simulations for bioavailability, time to reach the maximum concentration in plasma, area under curve, etc., were taken and analyzed.

2.5. Dynamics Simulation

In order to gain a more comprehensive understanding of the interaction and stability of the top complexes, dynamics simulations were conducted using the Desmond software from Schrodinger LLC. To ensure the geometries were correct, OPLS3e force field was employed and energy minimization was performed. The TIP3P water solvation model was selected in an orthorhombic box for boundary [38]. The complexes were neutralized by adding Na+ ions and used a salt concentration of 0.15 M. The simulation lasted for 100 ns with a recording interval of 10 ps and an energy of 5. The standard NPT ensemble class was used with a pressure of 1.01325 bar and a temperature of 300 K. After the simulation, interaction evaluation was done, to gain more insights. The Bio3D module of the R package [39] was used to conduct post-MD analysis. To align the trajectory to reference frame 0, the fit.xyz() function was employed. The resulting data was then used to generate plots of RMSD and principal component analysis (PCA). PC1 and PC2 were examined for visualizing the conformational dynamics [39].

3. Results

3.1. Pan-genome Analysis

The pan-genome was open and comprised of 1271 core gene clusters (12.88%), 3556 dispensable (36.03%), and 5,043 unique (51.09%) gene clusters. Maximum nucleotide polymorphism amounted to 0.11 in core and 0.40 in pan-genome (Supplementary Figure 2A). Lowest number of accessory genes were in strain 16CS0821-2 () and highest in RMI (). Unique gene count was lowest () in three strains RM4018, L350, and L352 while highest in Ab_CR1132 () (Table 1).

Most of the core genes were involved in translation, ribosomal structure, and protein genesis, while most of the dispensable genes were involved in replication, recombination, and repair (Supplementary Figure 2B). Majority of these proteins had an unknown function.

3.2. Drug Target Selection

Targeting the core genome of a bacterial species for drug development can be advantageous for drug targeting due to its conservation of function, essentiality, and broader coverage [40]. Core fraction revealed 2123 nonredundant protein sequences. DEG and CEG homologs were 947 and 855, respectively. Common sequences from both these databases () and dissimilar to human proteome were 398. Among these, 163 were dissimilar to gut proteome, while 59 matched DrugBank sequences, and 49 of these were linked with various KEGG pathways (Supplementary Table 1). Pyrimidine metabolism-related pathway enzymes were only two (WP_014468039.1 and WP_198407793.1). Cluster number 1186 was selected, depicting ODCase (accession no. WP_198407793.1) for further analysis. ODCase is considered a promising target for the development of antimicrobial and antifungal agents as it is involved in the biosynthesis of pyrimidine nucleotides, which are a class of nucleotides that include cytosine, uracil, and thymine [41, 42]. These nucleotides are essential building blocks of DNA and RNA. Inhibition of ODCase can lead to the accumulation of toxic intermediates in the pyrimidine biosynthesis pathway, which can ultimately result in cell death [43].

3.2.1. Structure Modeling and Docking

Modeled structure (Figure 1(a)) consisted of one sheet (with eight strands of topology -1X -1X -1X -1X -1X -1X -1X), seven beta-alpha-beta units, two beta bulges, 12 helices, nine helix-helix interacs, 15 beta turns, and two gamma turns. Overall ERRAT quality factor was 99.07, while Ramachandran plot (Figure 1(b)) showed 91.8% residues in core allowed, 7.7% in allowed, 0.5% in generously allowed, and not even a single residue in disallowed regions.

Orotidylic acid was taken as control (Figure 2(a)) and used for generating seven pharmacophoric features (Figure 2(b)).

Docking of the ODCase with orotidyclic acid (Figure 2(c)) gave a score of -3.39 (Table 2) and made interactions with 11 residues of ODCase, including two acidic, three basic, and six hydrogen bond interactions (Figure 2(d)). The binding affinity scores of best binding top three TCM compounds were -4.30 kcal/mol (ZINC70454134; IUPAC name: (3S,8S,9S,10R,13S,14R,17R)-17-[(1R)-2-hydroxy-1-[(2R)-5-(methoxymethyl)-4-methyl-6-oxo-2,3-dihydropyran-2-yl]ethyl]-10,13-dimethyl-1-oxo-2,3,4,7,8,9,11,12,14,15,16,17-dodecahydrocyclopenta[a]phenanthrene-3-sulfonic acid), -4.2 kcal/mol (ZINC85632684; IUPAC name: (E,3R,7R)-1-(5-amino-2-oxa-4,6-diazabicyclo[8.4.0]tetradeca-1(10),5,11,13-tetraen-7-yn-12-yl)-7-(4-hydroxyphenyl)-7-pentoxyhept-5-ene-3-sulfonic acid), and -4.02 kcal/mol (ZINC85632721; IUPAC name: (E,3R,7R)-7-hydroxy-7-(4-hydroxyphenyl)-1-[4-(4-phenylpiperidin-4-yl)oxyphenyl]hept-5-ene-3-sulfonic acid), depicting stronger binding compared to the control. ZINC70454134 and ODCase complex (Figure 3(a)) made one acidic and four basic interactions (Figure 3(b)). Five interactions were conserved between control and ZINC70454134 binding with ODCase (Asp59, Lys61, Val154, Thr119, and Pro176). ZINC85632684 and ODCase complex (Figure 3(c)) made three acidic and three basic interactions (Figure 3(d)). Except for Arg154, all other interactions were conserved between control and ZINC85632684. ZINC85632721 and ODCase complex (Figure 3(e)) made 16 interactions (Figure 3(f)). Except for Ser6, the rest of the interactions were similar between ZINC85632721 and control with ODCase.

3.3. ADMET Profiling

Best hits (ZINC70454134, ZINC85632684, and ZINC85632721) were subjected to ADMET profiling, along with control orotidylic acid. Predicted values for various pharmacokinetic and toxicity parameters were obtained and compared (Table 3). Water solubility of prioritized TCM compounds decreased in the order of ZINC70454134, ZINC85632684, and ZINC85632721. The Caco2 permeability across the intestinal lining was highest for ZINC70454134, followed by ZINC85632721 and then ZINC85632684.

VDss represents the theoretical volume into which a drug is distributed in the body, assuming it is uniformly distributed in plasma and tissues at equilibrium [44]. A negative value for VDss suggests that the drug is mainly confined to plasma. Fraction unbound refers to the proportion of a drug that is unbound or free in plasma and available for distribution to tissues [45]. A higher value for fraction unbound suggests that a larger proportion of the drug is available to exert its pharmacological effects. ZINC70454134 was predicted to be more restricted to plasma. It also had higher clearance or elimination from the body compared to ZINC85632684 and ZINC85632721. All the TCM compounds were a substrate of just one class of CYP450 enzymes, i.e., CYP3A4. Ames toxicity was null for all compounds. ZINC70454134 showed maximum tolerated dose in humans and no skin sensitization. Based on the data in Table 3, it appears that ZINC70454134 has the most favorable toxicity profile among the four compounds. It has a higher oral rat chronic toxicity LOAEL value than ZINC85632684 and ZINC85632721 and does not have any hepatotoxicity.

3.4. PBPK Analysis

PBPK was simulated under two conditions, normal (including male and female without any other disease), steatosis, cirrhosis, renal impairment, and in pregnant women. ZINC70454134 had the highest (%) and maximum concentration in blood (Cmax) values among the three compounds, indicating a relatively high degree of bioavailability and exposure in both normal and pregnant condition (Table 4). Time to reach peak exposure was almost the same for all compounds in both conditions. All the three TCM compounds had different values for Fa (%) and FDp (%) under normal and pregnant conditions, but the variation was not drastic. Under normal conditions, TCM compounds had a relatively high bioavailability, indicating that a large fraction of the administered dose reached the systemic circulation, but this value was lower in pregnant condition. This indicates a need for dosage adjustments in pregnant females. However, the peak concentration in blood, and the hours to reach peak exposure did not have much difference. The AUC values indicated a high total exposure to the compound over time and were even higher in pregnant condition. These results suggest a relatively high bioavailability and exposure levels under both normal and pregnant conditions, but a higher exposure under pregnant conditions. The first-pass metabolism did not show a large variation, and AUC values were notably lower in diseased patients, indicating reduced total exposure to the drugs over time. Overall, cirrhosis, renal impairment, and steatosis led to a pronounced decrease in absorption, Cmax, and AUC. Understanding these differences is important for tailoring drug regimens to individual patients with these specific medical conditions.

3.5. Dynamics Simulation

Simulation of control bound with ODCase as well as top-scoring compounds ZINC70454134, ZINC85632684, and ZINC85632721 was conducted. 28,202 atoms of the control-ODCase, 28,266 atoms of the ODCase-ZINC70454134, 28,238 atoms of the ODCase-ZINC85632684, and 28,271 atoms of the ODCase-ZINC85632721 complex, with respective 8154, 8161, 8153, and 8165 water molecules, were simulated in a box.

On the average, the RMSD did not exceed 3 Å for control-ODCase (Figures 4(a) and 4(b)), while the region around residue 130, 160, and 180 showed fluctuations larger than the rest of the ODCase (Figure 4(c)). These regions comprised of helix and loop region. The RMSD of ZINC70454134 and ZINC85632721 did not exceed 2.5 Å, but the RMSD of the ZINC85632684-ODCase complex reached up to 3 Å during 40-60 ns interval (Figures 5(a)5(c)). However, it decreased and became stable after 70 ns and kept around 1.5 Å afterwards. Lys30 and Leu60 made hydrogen bonds while Asp59 made a water bridge with the control compound, throughout the simulation time (Supplementary Figure 3A). Apart from this, Lys30, Asp59, Lys61, Pro176, Arg179, and Lys189 retained contact with ligand throughout simulation (Supplementary Figure 3B). No hydrogen bond or other interaction was retained for whole time of the simulation for ZINC70454134 (Supplementary Figure 4A, 4B). Leu7 and Glu10 retained hydrogen bond with ZINC85632684 for entire simulation time, along with interaction with Asp8 for more than 30% of simulation time (Supplementary Figure 4C, 4D). ODCase made interactions like hydrogen bond and water bridges with ZINC85632721 but did not retain for the entire time of simulation (Supplementary Figure 4E). Asp64 made interaction for more than 30% of time (Supplementary Figure 4F).

PCA and relative proportion of variance were calculated for all complexes. PC1 accounted for 31, ~36, and ~33%, while PC2 accounted for approximately 11, 20, and 20% of the total variance for ZINC70454134, ZINC85632684, and ZINC85632721, respectively. Proportion of variance (Supplementary Figure 5) was the amount of variability in the data explained by each principal component. The eigenvalue rank, on the other hand, referred to the relative importance of each principal component in explaining the overall variability of the data. PC1 captured the largest proportion of the variance, indicating that it explains most of the variability in the data. This could be due to large-scale conformational changes in the protein-ligand complex or due to fluctuations in the binding site caused by the ligand. In this scenario, the second principal component captured a smaller proportion of the variance but still was important in explaining the overall variability of the data. A higher eigenvalue rank of PC2 vs. PC1 indicates that it captures a more subtle, but still significant, aspect of the protein-ligand conformational dynamics, and the variation was higher for all complexes. Thus, all complexes had overall best binding affinity, visible from simulation results. The PCA proportion of variance and eigenvalue rank also validate this. Both of these are related but they do not always perfectly align. Nevertheless, they are important in understanding the contribution of each principal component to the overall variability of the protein-ligand complex’s conformational dynamics. In summary, the compounds showed stable binding, with average RMSD less than 3 Å over the course of 100 ns. The simulation replicates of the complexes showed the same dimensions of the trajectories and, hence, RMSD (Supplementary Figure 6).

4. Discussion

A. butzleri is allied with Campylobacter spp. and is a zoonotic pathogen [11]. It can cause diarrhea, nondiarrheal gastrointestinal illness, and arcobacteriosis [46]. Treatment of A. butzleri infection typically involves the use of antibiotics, such as fluoroquinolones and macrolides, but resistance to these antibiotics is becoming increasingly common [2, 7]. Therefore, prevention of A. butzleri infection through improved hygiene and sanitation practices, as well as proper food handling and preparation, is essential to control the spread of this emerging pathogen. Additionally, new compounds that could serve as drug for curbing A. butzleri infection are needed. For this purpose, a pan-genomics-allied subtractive genomics approach was used to identify novel therapeutic target and screen TCM natural product inhibitors against it.

Pan-genomics is a field of study that focuses on comparing the complete genetic makeup of individuals from a given population or species [47]. Understanding the pan-genome of a bacterial species can provide insights into its evolution, virulence, and antibiotic resistance potential, which can inform the development of new diagnostic tools, vaccines, and treatments [48]. One key application of pan-genomics in bacterial drug target identification is targeting the core fraction due to conservation and similar structures across different strains [17]. The percentage of core genes for A. butzleri was just 12.88%, while the higher percentage of dispensable and unique genes in the pan-genome, i.e., 36.03% and 51.09%, respectively, is likely due to horizontal gene transfer and genetic drift. The core fraction of A. butzleri comprised the set of genes present in all strains of this species, most of which were essential for the survival and replication and involved in basic metabolic processes and other fundamental functions. Among these, focus was on pyrimidine biosynthesis pathway enzyme cluster in A. butzleri, as it represents a promising target for the development of new drugs [49]. This pathway is essential for bacterial growth due to production of the building blocks for DNA/RNA synthesis and is highly conserved in bacteria, making it an attractive target for developing broad-spectrum inhibitors [50]. ODCase was selected in this pathway for its novelty as it is not much explored against phytochemical based inhibitors. It is a member of the family of pyridoxal -phosphate- (PLP-) dependent enzyme [51] and catalyzes the decarboxylation of orotidine -phosphate (OMP) to uridine -monophosphate (UMP), which is an important precursor for the synthesis of pyrimidine nucleotides [52].

Several classes of ODCase inhibitors have previously been identified, including ginkgolide [53] and novel C6 uridine replacements for antimalarial activity [54]. Anticancer activity of 6-azido-5-fluoro and 5-fluoro-6-iodo derivatives [55] has also been reported. However, to the best of authors’ knowledge, TCM-based discovery of phytochemicals against the ODCase of any pathogenic bacteria has not been reported yet. Therefore, the goal of this study was to identify natural compounds from TCM library that may be used as lead compounds for the development of new drugs against ODCase of A. butzleri. Pharmacophore mapping and then molecular docking were performed for screening a compound quickly and efficiently. This method helped identify potential inhibitors before their synthesis or isolation. Best docked compounds were ZINC70454134, ZINC85632684, and ZINC85632721, depicting stronger binding compared to the control. Their promiscuity was assessed to see if they could interact with multiple biological targets or proteins. However, they were nonpromiscuous and are best fit for inhibition of ODCase only (Figure 6).

MD simulation can help identify conformational changes in ligand binding site and even cryptic or allosteric binding sites [56]. It also helps predict small molecule binding energies and allow for the introduction of protein flexibility before or after a docking protocol, refining the structure of protein-ligand complexes in various environments, and ranking complexes with more accurate binding energy calculations [57]. The observed fluctuations in conformation and stability over simulation time provide dynamic view of interactions in drug development. A compound that maintains stability and favorable interactions with the target protein in a cellular environment is more likely to exhibit the desired therapeutic effects due to its continued occupation of the protein site. Salt concentration and overall environment can also influence electrostatic behaviour of a drug binding with receptor and simulation mimics that, along with watery environment of a cell which impact dissolution of a drug. Results suggest that ODCase-ZINC70454134 complex has a better overall binding affinity ( Å) than ODCase-ZINC85632684 complex, which is supported by the PCA proportion of variance and eigenvalue rank as well (Supplementary Figure 5). Hence, the stable binding and favorable conformational changes observed in the simulations suggest that ZINC70454134 has potential as a drug candidate.

The favorable outcomes in the computer-aided drug design study against A. butzleri, particularly with a focus on ZINC70454134 as an inhibitor, can be attributed to its optimized binding, stability, absorption, and safety parameters, distinguishing it among the screened compounds. Its superior characteristics, as evidenced by its distinct bioavailability compared to other TCM compounds, highlight its potential as a promising drug candidate. The observed impact of diseases such as cirrhosis, renal impairment, and steatosis on decreased absorption, Cmax, and AUC of compounds also underscores the importance of considering physiological conditions in drug performance evaluation. Moreover, MD simulation results supported the superior binding affinity of ODCase-ZINC70454134 and boosted evidence for its potential efficacy against A. butzleri.

It is important to note that the development of new drugs is a long and complex process, and more research is needed to determine the efficacy and safety of potential drug targets. In silico ADME/Tox prediction and PBPK simulation were therefore carried out to test various aspects of the possible compound action in body and their safety. Varied profile for simulated compounds in normal vs. pregnant condition was observed and may be due to changes in the physiology of the gastrointestinal tract and/or blood vessels during pregnancy could affect the absorption and distribution of drugs. For example, the reduced gastric acid secretion and increased intestinal motility during pregnancy [58] could lead to altered drug absorption. Previously lower effective exposure and enhanced metabolism of anti HIV-integrases were noted in pregnant women, possibly due to enhanced hormone release [59]. Altered renal or hepatic function during pregnancy could also impact the elimination of drugs from the body, since the glomerular filtration rate increases during pregnancy and may lead to increased clearance of some drugs [60]. In Table 4, it can be seen that the different compounds have different patterns of change between normal and pregnant conditions. Here, TCM compounds had a higher AUC(0-inf) value under pregnant conditions, indicating higher overall exposure and retention for the drug. This may be due to increased absorption and/or reduced metabolism or elimination of the drug during pregnancy. On the other hand, compound ZINC85632684 has a shorter Tmax value under normal conditions, which suggests faster absorption and/or distribution of the drug in nonpregnant individuals. These differences highlight the importance of considering the effects of pregnancy on drug pharmacokinetics, especially for drugs that are used during pregnancy. Most of the studies and clinical trials exclude pregnant women from the study, and little information is available for comparing the impact of various drugs on normal cohort vs. pregnant women. The results of this study and, in general, in silico simulation can, therefore, be used to inform further experimental or clinical studies, as well as to make predictions about the likely effects of changing doses, dosing regimens, or physiological conditions. Furthermore, in patients with cirrhosis, renal impairment, and steatosis, there was a decrease in overall bioavailability and absorption rates compared to the general population. This suggests that a smaller fraction of the drug reaches the systemic circulation in individuals with these medical conditions. The first-pass metabolism showed relatively consistent patterns, but AUC values were notably lower in patients with cirrhosis, renal impairment, and steatosis, indicating reduced total exposure to the drugs over time in these populations. Understanding these variations is crucial for tailoring drug regimens to individuals with specific medical conditions. Higher exposure levels in pregnant women may necessitate dosage adjustments to ensure therapeutic efficacy while avoiding potential adverse effects. In patients with cirrhosis, renal impairment, and steatosis, the pronounced decrease in absorption, Cmax, and AUC underscores the importance of individualized treatment approaches to achieve optimal therapeutic outcomes while minimizing the risk of side effects. These pharmacokinetic insights provide valuable guidance for clinicians in optimizing dosing regimens based on the specific physiological conditions and medical histories of patients. However, simulation-based findings may not perfectly replicate the complex dynamics and interactions that occur in vivo. Thus, further experimental validation, such as testing in cell lines and animal models, is necessary to confirm the efficacy and safety of drugs and, in this case, ZINC70454134 as a potential therapeutic agent against A. butzleri.

Antibiotics are still the first line of defense against severe bacterial infections and are often necessary to treat life-threatening conditions. Moreover, the overuse and misuse of antibiotics have led to the emergence of antibiotic-resistant bacteria, which is a major public health concern [61]. Phytochemicals can offer an alternative or complementary approach to antibiotic therapy, particularly for less severe infections or as a prophylactic measure [62, 63]. Research has shown that phytochemicals can have a range of antibacterial effects, such as inhibiting bacterial growth or disrupting bacterial membranes [64, 65]. Medicinal chemistry and structure-activity relationship studies can be used to further modify the structure of phytochemicals to enhance their antibacterial activity and reduce any potential toxic effects [66, 67]. It is also important to note that the dosage, purity, and quality of phytochemicals can vary, which can affect their effectiveness and safety [68]. Therefore, extensive experimentation and validation are needed before using phytochemicals for the treatment of bacterial infections.

5. Conclusion

In conclusion, the pan-genome analysis of the studied A. butzleri strains revealed a dynamic genomic landscape, characterized by a significant proportion of dispensable and unique gene clusters. The core genome exhibited conservation, and ODCase enzyme was selected for computational-aided drug design. This target has been studied for other pathogenic organisms as well. Three TCM inhibitors were prioritized and ZINC70454134 depicted strongest binding, compared to the control. All compounds were predicted as safe, but PBPK results were different for normal vs. pregnant cohort, highlighting the importance of considering the effects of pregnancy on drug pharmacokinetics. Dynamics simulations confirmed the stability of the ODCase-ligand complexes, with RMSD values indicating consistent binding over the simulation period. Principal component analysis supported the overall stability and variance of the complex dynamics, further validating the suitability of the selected compounds. The selection of a compound as a potential drug candidate involves a comprehensive evaluation of its pharmacological, toxicological, and physicochemical properties. Key factors include the compound’s efficacy in the targeted biological process or disease, safety profile, pharmacokinetics, and bioavailability. ZINC70454134 was selected based on its compliance of these parameters. The in silico simulations, as carried out in this study, provide a useful starting point for further investigation and optimization of drug targets. Hence, this multifaceted study provides foundation for development of new natural product inhibitors to control the spread of A. butzleri infection and may be replicated for other pathogens. These findings lay groundwork for further experimental validation and drug development efforts against A. butzleri as the true effectiveness and safety of the compound can only be determined through rigorous experimentation. Therefore, further testing on cell lines and in model organisms is proposed.

Data Availability

All the data used or generated in this study is provided as accession number or relevant information as tables in the manuscript.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Authors’ Contributions

ZB, SMA, and YW conceived and designed the experiments. YW and IA provided resources. ZB, SMA, and AM performed the experiments and wrote the manuscript. IA, AM, SMA, and AM validated findings and edited the manuscript. SMA and YW edited the manuscript and supervised the project. AM, IA, and ZB curated the data, validated findings, and procured funding for the project. All authors read and approved the final version of the manuscript. All authors approved the manuscript for publication in the current form.


Ibrar Ahmed is a recipient of the National Research Foundation of Korea (NRF) Grant No. RS-2023-00223245. The authors would like to thank the Deanship of Scientific Research, Qassim University, for supporting this research.

Supplementary Materials

Supplementary Figure 1: physiology information for pharmacokinetic parameters for (A) normal, (B) cirrhosis, (C) renally impaired, (D) steatosis, and (E) pregnant. Supplementary Figure 2: (A) Nucleic acid variation of core and pan-genome of A. butzleri. (B) COG categorization of pan, core, and accessory/dispensable genome. Supplementary Figure 3: (A) hydrogen bonds, hydrophobic, ionic, and water bridge contacts for the control and ODCase-simulated complex. (B) Interactions that occur more than 30.0% of the simulation time for control-ODCase simulation in the trajectory (0.00 through 100.00 nsec). Supplementary Figure 4: (A) hydrogen bonds, hydrophobic, ionic, and water bridge contacts for the ZINC70454134 and ODCase-simulated complex. (B) Interactions that occur more than 30.0% of the simulation time for ZINC70454134-ODCase simulation in the trajectory (0.00 through 100.00 nsec). (C) Hydrogen bonds, hydrophobic, ionic, and water bridge contacts for the ZINC85632684 and ODCase-simulated complex. (D) Interactions that occur more than 30.0% of the simulation time for ZINC85632684-ODCase simulation in the trajectory (0.00 through 100.00 nsec). (E) Hydrogen bonds, hydrophobic, ionic, and water bridge contacts for the ZINC85632721 and ODCase-simulated complex. (F) Interactions that occur more than 30.0% of the simulation time for ZINC85632721-ODCase simulation in the trajectory (0.00 through 100.00 nsec). Supplementary Figure 5: (A) eigenvalue plot of ODCase-ZINC70454134 complex. (B) Eigenvalue plot of ODCase-ZINC85632684 complex. (C) Eigenvalue plot of ODCase-ZINC85632721 complex. Supplementary Figure 6: (A) RMSD plot of ODCase-control complex for 100 ns. Ligand trajectory is shown in red and receptor in green. R1 depicts replica 1. (B) RMSD plot of ODCase-control complex for 100 ns. R2 depicts replica 2. (C) RMSD plot of ODCase-ZINC70454134 complex for 100 ns. R1 depicts replica 1. (D) RMSD plot of ODCase-ZINC70454134 complex for 100 ns. R2 depicts replica 2. (E) RMSD plot of ODCase-ZINC85632684 complex for 100 ns. R1 depicts replica 1. (F) RMSD plot of ODCase-ZINC85632684 complex for 100 ns. R2 depicts replica 2. (G) RMSD plot of ODCase-ZINC85632721 complex for 100 ns. R1 depicts replica 1. (H) RMSD plot of ODCase-ZINC85632721 complex for 100 ns. R2 depicts replica 2. Supplementary Table 1: KEGG pathways of drug targets mined from the core genome of A. butzleri. Empty cells indicate no information. (Supplementary Materials)