Initial screening of compounds identified by AtomNet in MIC and MBC experiments.

The AtomNet technology identified 74 compounds that could potentially interact with Ng-LdcA and 84 compounds that could potentially interact with Ng-LtgD. Initially, all the compounds were screened at a concentration of 50 µM in the standard broth microdilution assay against N. gonorrhoeae strain P9-17 (Fig. 1). We identified 26 compounds for Ng-LdcA (Fig. 1A) and 13 compounds for Ng-LtgD (Fig. 1B) that reduced the optical density readings by greater than 50% and these were chosen for further study. These compounds were then serially diluted in the broth microdilution assay and tested against strain P9-17 to determine the MIC50/90 values, shown in Supplementary Table 1 for Ng-LdcA and Supplementary Table 2 for Ng-LtgD. The individual MIC titration curves are shown in Supplementary Fig. 1 for the Ng-LdcA compounds and in Supplementary Fig. 2 for the Ng-LtgD compounds. Next, the MBC50/90 values were determined for only the compounds with the lowest MIC values, i.e. those most likely to show highest bactericidal effects (highlighted in Supplementary Tables 1, 2; Supplementary Fig. 3). This reductionist screening led to the final selection of the most active compounds, with three Ng-LdcA compounds and three Ng-LtgD compounds showing the highest activities against P9-17: the structures and properties of these compounds are summarised in Table 1. Compound 16 was the most active against Ng-LdcA, with a MIC50 value < 1.56 µM, MIC > 90 values of 3.125–6.25 µM and MBC50/90 values between 0.195 and 0.39 µM. In general, the Ng-LdcA compounds showed higher activity than the compounds directed against Ng-LtgD, of which compound 45 had MIC50 values of 1.56–3.125 µM, MIC > 90 values of 6.25–12.5 µM and MBC50/90 values between 3.125 and 6.25 µM (Table 1).

Fig. 1

Screening of compounds in a standard Minimum Inhibitory Concentration (MIC) assay. Ng-LdcA (n = 74) and Ng-LtgD (n = 84) compounds were diluted and tested at a single concentration of 50 µM in wells containing ~ 105 CFU of N. gonorrhoeae strain P9-17. Controls were bacteria alone, GC broth alone (control in the graphs), GC broth with DMSO (10µL/well) alone, and positive control was bacteria treated with ceftriaxone. DMSO alone has no effect on bacterial growth. Optical density was measured after 24 h incubation. The data are shown as the % reduction in optical density compared to the control bacteria alone. Data are from a representative experiment done twice. The red lines denote the 50% cut offline for selecting compounds for further analyses

Table 1 Summary of the structures, properties and MIC and MBC values for the selected Ng-LdcA and Ng-LtgD compoundsComputational modelling of Ng-LdcA and Ng-LtgD and analyses of their interactions with the hit compounds

Having identified the three Ng-LdcA (-16, -37 and -69) and three Ng-LtgD (-45, -52 and -69) compounds from the initial biological screen, we then did a deep computational study to model their interactions with their target enzymes.

Protein modelling and validation

In the absence of crystal structures, homology modelling has become the primary method for obtaining a three-dimensional (3D) representation of protein targets in recent years. Since no crystal structure of Ng-LdcA and Ng-LtgD was available, it was determined using homology modelling with AlphaFold. From the structural PDB BLAST results, the Ng-LdcA sequence showed the best alignment with the native Escherichia coli L,D-carboxypeptidase A, LdcA (PDB Id: 5Z01) [28] with a query coverage of 99% and sequence similarity of 41%. The Ng-LtgD sequence showed the best alignment with the native Escherichia coli Lytic transglycosylase D, LtgD (PDB Id: 1D0K) [29] with a query coverage of 88% and sequence similarity of 37%. These sequence similarity findings agreed with the initial Atomwise models. The observed gaps between the template and query proteins were found to be 1%. The remaining 11% of the sequence, which corresponded to regions at the C- terminal of the Ng-LtgD protein, was not covered by the alignment and was excluded from the final generated model of our protein. This portion of the sequence, including the C-terminal region (MEKRKILPLAICLAALSACTAMEARTPRANEAQAPRADEMKK ESRPAFDAAAVPVSDSGFAAN), was excluded from the final model as it did not include residues relevant to the catalytic domain or active site. The sequence alignment, highlighting the modelled regions and indicating the excluded parts, is provided in Supplementary Fig. 4.

The 3D structure of Ng-LdcA was modelled based on the sequence alignment with E. coli-LdcA PDB Id: 5Z01 (Fig. 2) and the 3D structure of Ng-LtgD was modelled based on the sequence alignment with E. coli-LtgD PDB Id: 1D0K (Fig. 2). The modelled structure for Ng-LdcA was superimposed on the 5Z01 template structure and the RMSD value was calculated at 1.70 Å and the modelled structure for Ng-LtgD was superimposed on the 1D0K template structure and the RMSD value was calculated at 1.20 Å. These calculations indicated significant backbone similarity between the modelled protein and the template structure. Both the modelled structures were then validated with the PDBsum server and the Ramachandran plot (Supplementary Fig. 5), which revealed for Ng-LdcA that 90% of residues were in the allowed region, 9.6% in the favoured region and 0.4% was present in disallowed regions, whereas for Ng-LtgD, 81.7% of residues were in the allowed region, 14% in the favoured region, 2.6% in generously allowed region and 1.7% was present in disallowed regions. These values suggested good overall stereo-chemical quality and stability of the models. ERRAT analysis for the overall model quality factor, indicated that both generated models were of extremely high-quality and reliability with respect to 3D structure. The models’ backbone conformation, non-bonded interactions and energy scores were well within the range of the high-quality models. In the case of Ng-LdcA, two homologous chains A and B constituted the dimeric protein, whereas Ng-LtgD was a monomer.

Fig. 2

Homology-based models of Ng-LdcA and Ng-LtgD. The E. coli 5Z01 template-based model of Ng-LdcA and the E. coli 1D0K template-based model of Ng-LtgD were generated by AlphaFold. The Ng-LdcA homo-dimeric protein consists of two identical chains (Chain A in green and Chain B in red), whereas Ng-LtgD is a monomer (in green)

Active site analysis and validation

The identification of potential binding site(s) for the ligands provides an insight into the active site regions on the proteins and gives an idea of the interacting residues as well as the crucial interactions between the proteins and ligands. Identifying the binding site is essential for structure-based virtual screening of compound libraries. To access the binding cavity in the modelled Ng-LdcA protein, the Ng-LdcA moiety mostly contacted amino acid residues in the deep cleft between the enzyme’s two domains. As shown in Fig. 3 the residues from Chain A: Tyr267, Arg268, Arg271, Tyr302, Asp303, and from Chain B: Gly69, Phe70, Glu72, Arg103, Arg133, Gly135, Tyr136, Ser165, Asn236, Ser238, Val239, Asp261, Val262, Glu264, formed the main region of the binding site and were conserved across multiple Neisseria species and homologs in other bacteria. To access the binding cavity in the modelled Ng-LtgD protein, as shown in Fig. 3, the residues from the protein chain of Met101, Ile157, Glu158, Asn160, Asn164, Arg184, Tyr213, Ala214, Gln221, Phe222, Met223, Ser226, Tyr256, Gln340, Tyr341, Asn342, His343, Tyr347, formed the main region of the binding site and were similarly conserved. The conserved nature of active site residues across different species gives us confidence in the relative accuracy of both the models [55,56,57]. The identified active sites were confirmed with the help of previously reported literature for Ng-LdcA and Ng-LtgD protein complexes and validated by repeating the targeted docking and creating the binding cavity using the same residues.

Fig. 3

Identification of binding pockets with DeepFold. The binding pocket (active site) predicted by DeepFold is shown for Ng-LdcA and Ng-LtgD. The active site regions show the conserved amino acid residues predicted to be involved in protein–ligand binding and conserving the architecture of the binding cavity. See text for description of the amino acids involved

Molecular docking

This was used to assess the binding ability of the three Ng-LdcA and three Ng-LtgD hit compounds to the target Ng-LdcA and Ng-LtgD proteins. The compounds were docked and modelled onto Ng-LdcA and Ng-LtgD using AutoDock Vina to evaluate their binding affinity and to investigate how precisely the ligands docked into the binding region at the protein surface. All the successfully docked compounds exhibited favourable binding interactions with Ng-LdcA and Ng-LtgD (Table 2) and the docked compounds are shown in Fig. 4 for Ng-LdcA and Fig. 5 for Ng-LtgD with their interactions and best docked pose selected based on the docking scores. These compounds showed good interactions with some of the substrate-binding amino acids, such as Gly69, Phe70, Glu72, Arg103, Arg133, Gly135, Tyr136, Ser165, Asn236, Ser238, Val239, Asp261, Val262, Glu264, Tyr267, Arg268, Arg271, Tyr302 and Asp303 for Ng-LdcA, and Met101, Ile157, Glu158, Asn160, Asn164, Arg184, Tyr213, Ala214, Gln221, Phe222, Met223, Ser226, Tyr256, Gln340, Tyr341, Asn342, His343 and Tyr347 for Ng-LtgD, and they fit well into the active site cavity of their respective proteins. The significant interactions of these compounds plotted by LigPlot+ are shown in Fig. 6 for Ng-LdcA and Fig. 7 for Ng-LtgD.

Table 2 Characteristics of hit compounds with their docking, free energy and binding affinity scoresFig. 4

Molecular docking studies for Ng-LdcA. The left-hand exploded diagram shows the binding site at the surface of modelled protein for Ng-LdcA. On the right-hand side are images of compounds Ng-LdcA-16, -37 and -69 docked in the binding cavity. All the docked ligands exhibited good docking scores and retained all the conserved residues in the protein–ligand binding interaction

Fig. 5

Molecular docking studies for Ng-LtgD. The left-hand exploded diagram shows the binding site at the surface of modelled protein for Ng-LtgD. On the right-hand side are images of compounds Ng-LtgD-45, -52 and -69 docked in the binding cavity. All the docked ligands exhibit good docking scores and retained all the conserved residues in the protein–ligand binding interaction

Fig. 6

LigPlot+ 2D interaction diagrams for Ng-LdcA. The 2D interaction diagrams of the compounds docked with modelled Ng-LdcA proteins, plotted using LigPlot+. The interactions of all the molecules obtained by molecular docking are shown for compounds Ng-LdcA-16, -37 and -69

Fig. 7

LigPlot+ 2D interaction diagrams for Ng-LtgD. The 2D interaction diagrams of the compounds docked with modelled Ng-LtgD proteins, plotted using LigPlot+. The interactions of all the molecules obtained by molecular docking are shown for compounds Ng-LtgD-45, -52 and -69

Molecular dynamics simulations (MDS) of compounds

To gain insights into the stability and dynamic properties of the protein–ligand complexes, explicit solvent MDS were done. MDS provided detailed insight into protein–ligand interactions in motion, contributing to their stable bound conformation and visualizing the effect of ligand binding on protein conformational changes. MDS for up to 50 ns were done for control (protein alone) as well as for the compounds docked with protein. The trajectories were also analysed for all the protein–ligand docked complexes after completion of the simulation run. The time evolution of the Root Means Square Deviation (RMSD) during MDS is used to monitor protein stability. The distributional probability of RMSD up to 50 ns trajectories is shown in Fig. 8. The mean RMSD values of Ng-LdcA-control complex, Ng-LdcA-16, Ng-LdcA-37 and Ng-LdcA-69 complex were ~ 0.2–0.4 nm, which shows the structural stability throughout the MD run [58, 59], where the RMSD values of Ng-LdcA-16 were found to be more stable (Fig. 8A, top panel). We also calculated the radius of gyration (Rg) value, a parameter directly associated with the overall conformational changes in the structure of the enzyme upon ligand binding, to further validate our results. It also revealed the stability, compactness, and folding behaviour of the structure. We calculated the Rg values of all the selected compounds and the reference complex to determine their compactness. The average Rg values for the Ng-LdcA-control complex, Ng-LdcA-16, Ng-LdcA-37 and Ng-LdcA-69 complex were below 2.6 nm, were stable after 15 ns implying increased compactness, and improved binding (Fig. 8B, top panel). In addition to RMSD and Rg, the number of hydrogen bonds generated between protein and ligand throughout the simulation duration was determined. Figure 8D (top panel) displays the graphs of these H-bonds formed between the protein and the corresponding ligand throughout the 50 ns simulation run. The average number of H-bonds formed between Ng-LdcA-16 is highest followed by Ng-LdcA-37. We have also analysed an important parameter, which is SASA (Solvent-Accessible Surface Area), which is an approximate surface area of a biomolecule that is accessible to a solvent with respect to the simulation time. Figure 8C (top panel) indicates that for Ng-LdcA-16 and Ng-LdcA-control complex the solvent accessible surface was optimally acquired and within the acceptable range [60, 61].

Fig. 8

Plots to investigate the energy deviation, conformation stability and surface area accessible during simulation for Ng-LdcA and Ng-LtgD proteins and ligands in bound state with the protein. A represents RMSD, B represents radius of gyration, C represents Solvent Accessible Surface Area (SASA), and D represents number of hydrogen bonds. Black colour shows control, whereas red, green and blue colour shows 1st, 2nd and 3rd ligands respectively. Ng-LdcA 1st ligand = compound -16, 2nd = compound -37, 3rd = compound -69. Ng-LtgD 1st ligand = compound -45, 2nd = compound -52, 3rd = compound -69

Similarly, in the case of Ng-LtgD, the distributional probability of RMSD up to 50 ns trajectories is shown in Fig. 8 (bottom panel). The mean RMSD values of Ng-LtgD-control complex, Ng-LtgD-45, Ng-LtgD-52 and Ng-LtgD-69 complex were originated from around 0.2 nm and ranges up-to 0.6 nm, which shows the structural stability throughout the MD run, where the RMSD values of Ng-LtgD-69 were found to be more stable (Fig. 8A, bottom panel). We also calculated the Rg values of all the selected compounds and the reference complex to determine their compactness. The average Rg values for the Ng-LtgD-control complex, Ng-LtgD-45, Ng-LtgD-52 and Ng-LtgD-69 complex were below 2.2 nm, were stable after 10 ns implying increased compactness, and improved binding (Fig. 8B, bottom panel). Figure 8D (bottom panel) displays the graphs of the H-bonds formed between the protein and the corresponding ligand throughout the 50 ns simulation run. The average number of H-bonds formed between Ng-LtgD-69 is highest followed by Ng-LtgD-52. In addition, Fig. 8C (bottom panel) indicates that for Ng-LtgD-69 and Ng-LtgD-control complex the SASA was optimally acquired and within the acceptable range (50–200 nm2) [60, 61].

Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA)

Conformations from the last 10 ns MDS run were used to gain more insight into the structural dynamics of protein–ligand complexes. To validate the MDS findings, free energies of docked complexes from all the potential inhibitor small molecules along with control were calculated using MM-PBSA with the g_mmpbsa tool (Fig. 9). All the complexes exhibited negative binding energy in the MM-PBSA threshold. As shown, in case of Ng-LdcA, out of the three compounds, Ng-LdcA-16 exhibited a better net binding energy score which may be comparable with the binding energy of control. The binding energy score for Ng-LdcA-16 as depicted by MM-PBSA, exhibited better protein–ligand binding compared to two other compounds. Similarly, for Ng-LtgD, out of the three compounds, Ng-LtgD-69 exhibited a better net binding energy score, which may be comparable with the binding energy of control. The binding energy score for Ng-LtgD-69 depicted by MM-PBSA exhibited better protein–ligand binding, compared to the two other compounds.

Fig. 9

Binding free energy calculations using MM-PBSA tool for the potential small molecule inhibitor compounds along with the control. Colour coding is represented in the figure. Ng-LdcA 1st ligand = compound -16, 2nd = compound -37, 3rd = compound -69. Ng-LtgD 1st ligand = compound -45, 2nd = compound -52, 3rd = compound -69

Cytotoxicity of the compounds for human cells in vitro

We tested the cytotoxicity of the three Ng-LdcA compounds and three Ng-LtgD compounds on human Chang conjunctival epithelial cells using a standard resazurin assay, with cells treated for 18 h with 50 μM of the compounds and measurements of cell death after the addition of resazurin recorded at 4 and 18 h (Fig. 10A). There was no significant cytotoxicity shown by compounds Ng-LdcA-16, -37 and -69, nor by Ng-LtgD-45 and -52 after 4 h of cell viability measurement (P > 0.05), whereas there was a 50% reduction in cell viability recorded with Ng-LtgD-69 (P < 0.05). However, cell viability appeared to be restored when the cells were examined after 18 h, most likely due to rapid division of surviving cells unaffected by the presence of any compound.

Fig. 10

A Cytotoxicity of compounds. Human Chang conjunctival cells were treated with 50 µM (final concentration) of Ng-LcdA-16, -37 and -69 and Ng-LtgD-45, -52 and -69 compounds and cytotoxicity was measured using a standard resazurin assay. Controls included untreated cells, DMSO alone, medium alone, cell lysis (i.e. induced death) and ceftriaxone. The columns represent the means, and the error bars the standard error of the means of three independent experiments. B Determination of time to kill for compound Ng-LdcA-16. Bacteria (105 CFU/well, n = 3) were treated with 50 µM (final concentration) of Ng-LdcA-16 and viable counts were made over time. Controls were bacteria alone, and bacteria treated with ceftriaxone (50 µM final concentration). Data are from one representative experiment of experiments done at least twice

Specificity of the compounds for killing N. gonorrhoeae

We tested the hypothesis that the three Ng-LdcA compounds and three Ng-LtgD compounds did not kill bacteria of other genera and thus would be specific for gonococci. The majority of the Ng-LdcA or Ng-LtgD compounds tested at 50 µM did not kill P. aeruginosa or a variety of staphylococcal strains (Supplementary Table 3). Ng-LtgD-69 recorded a MIC50 of 50 µM against S. aureus NCTC8325.4, but when examined on other staphylococcal species, none of the Ng-LtgD compounds were active (Supplementary Table 3). Against the commensal L. gasseri, Ng-LdcA-69 reported a MBC50 of 50 µM and Ng-LtgD-45 and Ng-LtgD-69 compounds recorded MBC50/90 values of 12.5 µM (Supplementary Table 3). The MIC and MBC titration curves for compounds tested against P. aeruginosa, Staphylococcus spp and L. gasseri are shown in Supplementary Figs. 6, 7 and 8.

Anti-gonococcal activity of compounds Ng-LdcA-16 and Ng-LtgD-45

Based on the initial MIC and MBC screening study and the computational modelling, we chose the best performing compounds Ng-LdcA-16 and Ng-LtgD-45 as exemplars for further biological studies. All the studies with Ng-LdcA-16 were done with compound provided by Atomwise (sourced from Enamine), whereas additional lots were necessary to complete the experiments with Ng-LtgD-45 (sourced from Chemspace and Enamine).

The ability of Ng-LdcA-16 to kill several gonococci from the FDA/CDC AR gonococcal biobank, which were reported to show increased resistance to ceftriaxone, was examined (Table 3; Supplementary Fig. 9). The reported MIC90 value for ceftriaxone for all these isolates was 0.19 µM, whereas the MIC90/MBC90 values for compound Ng-LdcA-16 were ~ 4–16 fold higher at between 0.78 and 3.125 µM. Interestingly, strain P9-17 was particularly sensitive to ceftriaxone, whereas the effects of Ng-LdcA-16 on this strain were like those observed on the FDA/CDC isolates. In time-to-kill assays, Ng-LdcA-16 killed 100% of gonococci by 6 h (Fig. 10B). For comparison, the antibiotic ceftriaxone was able to kill 100% of gonococci by 3 h (Fig. 10B).

Table 3 Activity of Ng-LdcA-16 against different gonococcal strains with reported resistance to ceftriaxone

In our study, the additional new lots of Ng-LtgD-45 from Chemspace and Enamine purchased to finish the remaining experiments (i.e. testing against other gonococcal strains, against other Neisseria spp. and time-to-kill assays), surprisingly, were not effective in the MIC and MBC assays against P9-17 gonococci growing in either supplemented GC broth or in a simple 1% (w/v) proteose peptone medium. This was not the case with any additional lots of Ng-LdcA-16. Thus, to test new batches of Ng-LtgD-45 against P9-17 and the gonococcal isolates from the FDA/CDC AR biobank, we used a bactericidal MBC assay previously used to test other compounds that were inactivated by complex growth media [53, 54]. In this established assay, we treated 105 CFU of gonococci with various doses of the compound for 1 h in a PBSB solution and then examined survivors with viable counting. The MBC50 and MBC90 values are shown in Table 4 and the titration curves in Supplementary Fig. 10. In this assay, Ng-LtD-45 was highly active against all the strains, with MBC50 values ranging from 0.02–0.39 μM and MBC > 90 values from 0.05–0.78 μM (Table 4). Again, P9-17 was highly sensitive to this compound with a MBC value of 0.006 μM and MBC > 90 value of 0.048 μM. Under the test conditions, the bactericidal effects occurred within the 1 h of the exposure of the bacteria to this compound.

Table 4 Activity of Ng-LdcA-45 against different gonococcal strains with reported resistance to ceftriaxone

Finally, it is possible that the Ng-LdcA-16 and Ng-LtgD-45 compounds could have activity against other Neisseria spp and so they were tested against reference strains N. meningitidis MC58 and N. lactamica strain Y92-1009. Prior to this, we examined the amino acid similarities of both enzymes between the three species by Clustal. There was ~ 89% similarity in the amino acid sequences of LdcA of the three species and 94% similarity for LtgD (Supplementary Fig. 11). In addition, the amino acids belonging to the binding sites interacting with the compounds were completely conserved between the species for both LdcA and LtgD (Supplementary Fig. 11, shown in green). Thus, our hypothesis was that the two compounds would be equally effective against N. lactamica and N. meningitidis. Ng-LdcA-16 showed highly significant efficacy against N. meningitidis MC58 with MIC 50/ > 90 values of 0.1 μM and against N. lactamica with MIC50 values of 1.56 µM and MIC > 90 values of ~ 25 µM, and MBC50/90 values of 0.78—1.56 µM (Fig. 11A, B). Using the viable count assay with Ng-LtgD-45, high activity was observed against meningococci, with a MBC > 90 value of 0.39 µM, and MBC50 values < 0.01 µM (Fig. 11C). Against N. lactamica, MBC50 and MBC > 90 values of 0.39 µM and 0.78 µM were recorded (Fig. 11D).

Fig. 11

Activity of Ng-LdcA-16 and Ng-LtgD-45 compounds against other Neisseria spp. A Meningococci and B N. lactamica (105 CFU/well, n = 3) were treated with various concentrations of compound Ng-LdcA-16 in the standard MIC and MBC assay, over 24 h. C Meningococci and D N. lactamica (105 CFU/well, n = 3) were treated with various concentrations of compound Ng-LtgD-45 in the MBC assay in PBSB for 1 h with viable counting. Symbols represent the mean and any error bars the standard error of the means from three independent experiments