*3.8. Molecular Docking of the CYP128 Model*

Non-covalent docking of menaquinone-9 (MK9), a known substrate of CYP128A1, and different azole compounds (Figure S2) was done with MOE's dock utility. Prior to this, all compounds were "washed" using MOE, i.e., we generated the most dominant protonation state of each compound at neutral pH, computed its atomic partial charges, and minimized the generated 3D structure. Docking was performed into the catalytic site of our CYP128A1 model, setting the placement method to "Triangle Matcher" and the scoring and rescoring methods to "London dG" and "GBVI/WSA dG", respectively. After docking, the ligand structures were further refined using the fixed receptor option. This refinement step entails energy minimization using the conventional Amber10:Extended Huckel Theory molecular mechanics force field to take electronic effects into account. For each compound, the top 10 complexes as identified from GBVI/WSA-dG scores were considered for further MD equilibration (see Section 3.10). In the case of MK9, since the molecule was predicted to undergo omega hydroxylation via the heme group, only poses with proper orientation (hydrophobic tail facing the heme group) were kept, resulting in six (out of 10) poses being selected.

#### *3.9. Density Functional Theory*

MD parameters—e.g., partial charges and force constants—for non-standard residues like the heme group are not provided by standard MD force fields. Hence, to run further equilibration of our docked complexes, it was necessary to generate those parameters via quantum calculations. Such calculations were performed using Gaussian 09 (g09) together with the Metal Center Parameter Builder (MCPB.py) available in the Amber16 package [55,56]. MCPB.py was applied to create the correct g09 input files by including the heme ferric cation and its nearby residues from our 3D structures. g09 was successively utilized for geometry optimization, force constant calculation and Merz-Kollman RESP charge calculation of the selected atoms. Every calculation was performed at the B3LYP/6-31G\* level of theory. Finally, the MCPB.py program was re-applied to fit RESP charges and generate parameters compatible with Amber's ff14SB force field. Note that a partial charge of 0.250*e* was

calculated for the ferric cation in the heme group. Keeping in mind that the partial charge of a metal ion is less than 2*e*, no matter how big its formal charge (oxidation state) is, our result looks reasonable.

#### *3.10. Molecular Dynamics*

While MCPB.py and g09 were used to get the correct force field parameters for the ferric core of the heme and nearby atoms, Amber's antechamber utility [55,56] was applied to generate MD parameters for MK9 and azole compounds. Amber's tleap program was applied to solvate all our docked complexes in TIP3P water and to generate the correct topology file to conduct MD simulations using the ff14SB force field. For each complex, a rectangular box with at least 10 Å distance between the box edges and the protein was considered, resulting in about 16,380 water molecules in each case. Na<sup>+</sup> and Cl<sup>−</sup> ions were added to approximate 0.15 M concentration as well as to neutralize the system. Minimization, NVT, NPT, and MD production runs were all performed with Amber's pmemd utility. Minimization of each structure was carried out in two phases, both using the steepest descent and conjugate gradient methods successively. Briefly, minimization was done in 10,000 minimizations steps on hydrogens and solvent atoms only, i.e., by restraining the protein–ligand complexes. Next, a 20,000-step minimization was run without restraints. The structures were then equilibrated in the NVT ensemble during 20 ps and in the NPT ensemble during 40 ps, setting the temperature to 298 K and the pressure to 1 bar. Finally, MD production was run to relax each complex for at least 20 ns. The stability of each complex was assessed by checking if the RMSD of both the protein and the ligand reached a plateau by the end of the simulation (see Figure S3). In general, we found 20 ns sufficient to equilibrate the complexes, although in some cases, an extra 5–10 ns may have been required to equilibrate the structure fully.

#### *3.11. Redocking of Compounds*

MK9 and the azoles compounds were re-docked on their corresponding set of equilibrated structures using the MOE's dock (six structures for MK9, 10 structures for each of the azole drugs). The same options as in the first docking step were used (see Section 3.8).

#### **4. Conclusions**

The availability of quite a large number of bacterial genome sequences and different bioinformatics tools gives us an opportunity to understand the role of different genes/proteins in bacterial communities at large rather than confining the results to a single bacterium. In this study, we utilized such information and tools to understand the CYP128 P450 family profiles in different mycobacterial species and understand its structure. The study revealed interesting aspects; for example, P450 was found to be highly conserved in the MTBC species causing lung disease in humans and other animals, indicating the important function of this enzyme during the latent phase of these organisms, as this P450 was found to be expressed during such stage [11,12]. Furthermore, only 12–34% of non-MTBC species that have this P450 strongly support its important role during the latent phase of MTBC species. Contrasting CYP128 subfamily profiles between MTBC and four other different mycobacterial categories revealed by this study suggests that CYP128 P450s may play different roles in different category species, as previously observed for the CYP53 family in fungi [25]. Molecular dynamic simulation of CYP128A1 interactions with different ligands revealed that this P450 has the highest binding affinity to its substrate compared to azole drugs. Binding of azole drugs to CYP128A1 protein models further shows the complex biology of mycobacterial species, as azole drugs have been found to be effective against *Mycobacterium tuberculosis* [18–20], indicating the expression of *CYP128A1* in the latent phase [11,12], and possible binding of azole drugs during experimentation did not resulted in a hypervirulent bacterium. This may be due to the inhibition of essential P450s such as CYP121A1 and CYP125A1, leading to the death of *M. tuberculosis* H37Rv. One of the interesting observations of this study is that more than one copy of this gene was present in some species and one of these genes was found not to be part of the gene cluster. This poses a fascinating evolutionary question on the

need for having more than one copy of this gene. Research on the CYP128 P450s that are not part of gene clusters would be interesting in discerning other roles of CYP128 P450s in mycobacterial species. It is also important that *CYP128* gene clusters were found to have *CYP124A1* or both *CYP124A1* and *CYP121A1*, indicating the collective efforts of these P450s in generating complex lipids in mycobacterial species, especially in MTBC species, as none of the other four mycobacterial category species has this type of combination in its *CYP128* gene clusters. Future studies should include cloning of the entire *CYP128A1* gene cluster (containing *CYP121A1* and *CYP124A1*) and analyzing the effect of the cluster molecule in *M. tuberculosis* pathogenesis.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1422-0067/21/14/4816/s1, Figure S1: 3VRM/CYP128A1 alignement., Figure S2: 2D structures of substrate (menaquinone 9) and azole compounds used in the study, Figure S3: RMSD *vs* time during MD simulations of MK9-CYP128A1 complexes (colors are explained in the legend of each figure), Table S1: Comparative amino acid conservation acid analysis of CYP128 P450 family with top 10 ranked families, Supplementary Dataset 1: Comprehensive comparative analyses of CYP128 P450s and those associated with secondary metabolite biosynthetic gene clusters in mycobacterial species, Supplementary Dataset 2: CYP128 P450 sequences identified and annotated in mycobacterial species, Supplementary Dataset 3: A high-resolution phylogenetic tree of CYP128 P450s.

**Author Contributions:** Conceptualization, K.S.; data curation, N.S.N., Z.E.C., W.C., J.-H.Y., D.R.N., J.A.T., J.P. and K.S.; formal analysis, N.S.N., Z.E.C., W.C., J.-H.Y., D.R.N., J.A.T., J.P. and K.S.; funding acquisition, K.S.; investigation, N.S.N., Z.E.C., W.C., J.-H.Y., D.R.N., J.A.T., J.P. and K.S.; methodology, N.S.N., Z.E.C., W.C., J.-H.Y., D.R.N., J.A.T., J.P. and K.S.; project administration, K.S.; resources, K.S.; supervision, K.S.; validation, N.S.N., Z.E.C., W.C., J.-H.Y., D.R.N., J.A.T., J.P. and K.S.; visualization, N.S.N., W.C., J.P. and K.S.; writing—original draft, N.S.N., W.C., J.-H.Y., D.R.N., J.P. and K.S.; writing—review and editing, N.S.N., W.C., J.-H.Y., D.R.N., J.P. and K.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work presented in this article is part of a National Research Foundation (NRF), South Africa grant awarded to Khajamohiddin Syed (Grant No. 114159), where all the international authors involved in the study are listed as international collaborators. Master's study bursaries were provided for Nokwanda Samantha Ngcobo (year 2019) and Zinhle Edith Chiliza (year 2018) from the same grant. Zinhle Edith Chiliza also thanks the NRF, South Africa for a DST-NRF Innovation Master's Scholarship for the year 2019 (Grant No. 117182). Khajamohiddin Syed expresses sincere gratitude to the University of Zululand Research Committee for funding (Grant No. C686) and for the laboratory facilities.

**Acknowledgments:** The authors want to thank Barbara Bradley, Pretoria, South Africa for English language editing.

**Conflicts of Interest:** The authors declare no conflict of interest.
