An Insight into the Structural Requirements and Pharmacophore Identification of Carbonic Anhydrase Inhibitors to Combat Oxidative Stress at High Altitudes: An In-Silico Approach

Carbonic anhydrases (CA) inhibitory action could be linked to the treatment of a number of ailments, including cancer, osteoporosis, glaucoma, and several neurological problems. For the development of effective CA inhibitors, a variety of heterocyclic rings have been investigated. Furthermore, at high altitudes, oxygen pressure drops, resulting in the formation of reactive oxygen and nitrogen species, and CA inhibitors having role in combating this oxidative stress. Acetazolamide contains thiadiazole ring, which has aroused researchers’ interest because of its CA inhibitory action. In the present study, we used a number of drug design tools, such as pharmacophore modeling, 3D QSAR, docking, and virtual screening on twenty-seven 1,3,4-thiadiazole derivatives that have been described as potential CA inhibitors in the literature. An atom-based 3D-QSAR analysis was carried out to determine the contribution of individual atoms to model generation, while a pharmacophore mapping investigation was carried out to find the common unique pharmacophoric properties required for biological activity. The coefficient of determination for both the training and test sets were statistically significant in the generated model. The best QSAR model was chosen based on the values of R2 (0.8757) and Q2 (0.7888). A molecular docking study was also conducted against the most potent analogue 4m, which has the highest SP docking score (−5.217) (PDB ID: 6g3v). The virtual screening revealed a number of promising compounds. The screened compound ZINC77699643 interacted with the amino acid residues, Pro201 and Thr199, in the virtual screening study (PDB ID: 6g3v). These interactions demonstrated the significance of the CA inhibitory activity of the compound. Furthermore, ADME study revealed useful information regarding compound’s drug-like properties. Therefore, the findings of the present investigation could aid in the development of more potent CA inhibitors, which could benefit the treatment of oxidative stress at high altitudes.


Introduction
At high altitudes, oxygen pressure drops, resulting in the formation of reactive oxygen and nitrogen species (RONS), which have been linked to a variety of oxidative stress then into pIC50 value. To generate the predicted pIC 50 values for 3D-QSAR analysis, the complete dataset was divided in 7:3 training and test sets. The complete process of the present study is depicted in Figure 1. Table 1. Compounds taken in QSAR study with CA inhibitory activity (IC 50 and pIC 50 value) [21].

S. No.
Compounds Structures IC 50 Value (µM) pIC 50 Value 1 4a Figure 1. The schematic workflow of 3D QSAR, pharmacophore model, docking, and virtual screening performed in the present study.  Figure 1. The schematic workflow of 3D QSAR, pharmacophore model, docking, and virtual screening performed in the present study.  Figure 1. The schematic workflow of 3D QSAR, pharmacophore model, docking, and virtual screening performed in the present study.

Preparation of Ligands
The LigPrep [23] module was used to prepare the ligand via Maestro, which was specifically configured to provide input structures for the Glide and PHASE modules. Clean up wizard can process one ligand per second at a time, effectively converting large datasets from 2D to 3D structures, as well as essential steps in pharmacophore development and docking studies, utilising unique algorithms. To discover the optimal alignment and common features for 3D QSAR model generation, the molecules were aligned based on their most common core structure, as shown in Figure 2.

Preparation of Ligands
The LigPrep [23] module was used to prepare the ligand via Maestro, which was specifically configured to provide input structures for the Glide and PHASE modules. Clean up wizard can process one ligand per second at a time, effectively converting large datasets from 2D to 3D structures, as well as essential steps in pharmacophore development and docking studies, utilising unique algorithms. To discover the optimal alignment and common features for 3D QSAR model generation, the molecules were aligned based on their most common core structure, as shown in Figure 2. Dataset of 27 ligands were taken in the present study [21]. All of the ligands' structures were drawn in ChemDraw Ultra 12.0 software [22] and saved in '.mol' format. The biological activity was reported in literature in terms of IC50 values (in µM), which were converted to pIC50 for QSAR studies, as shown in Table 1. pIC50 is represented as the negative logarithmic value of IC50, so first the IC50 value was converted into molar concentration, then into pIC50 value. To generate the predicted pIC50 values for 3D-QSAR analysis, the complete dataset was divided in 7:3 training and test sets. The complete process of the present study is depicted in Figure 1.

Preparation of Ligands
The LigPrep [23] module was used to prepare the ligand via Maestro, which was specifically configured to provide input structures for the Glide and PHASE modules. Clean up wizard can process one ligand per second at a time, effectively converting large datasets from 2D to 3D structures, as well as essential steps in pharmacophore development and docking studies, utilising unique algorithms. To discover the optimal alignment and common features for 3D QSAR model generation, the molecules were aligned based on their most common core structure, as shown in Figure 2.
The LigPrep [23] module was used to prepare the ligand via Maestro, which was specifically configured to provide input structures for the Glide and PHASE modules. Clean up wizard can process one ligand per second at a time, effectively converting large datasets from 2D to 3D structures, as well as essential steps in pharmacophore development and docking studies, utilising unique algorithms. To discover the optimal alignment and common features for 3D QSAR model generation, the molecules were aligned based on their most common core structure, as shown in Figure 2.

Pharmacophore Mapping
Pharmacophore mapping is the process of defining and placing unique pharmacophoric features as well as using alignment algorithms to overlay 3D conformations [24]. The PHASE module of the Schrodinger Maestro software was used to carry out the pharmacophore mapping study. There are various features present in developed pharmacophore hypothesis, such as hydrogen bond acceptors (A), hydrogen bond donor (D), aromatic ring (R), and a hydrophobic (H) group. The top ranked hypothesis consists of two unique pharmacophoric features based on structural requirements for biological activity: hydrogen bond donor (D) and ring aromaticity (R). As shown in Figure 3, each feature in this pharmacophoric map signifies a common particular portion associated with the selected compounds.

Pharmacophore Mapping
Pharmacophore mapping is the process of defining and placing unique pharmac phoric features as well as using alignment algorithms to overlay 3D conformations [24 The PHASE module of the Schrodinger Maestro software was used to carry out the pha macophore mapping study. There are various features present in developed pharmac phore hypothesis, such as hydrogen bond acceptors (A), hydrogen bond donor (D), ar matic ring (R), and a hydrophobic (H) group. The top ranked hypothesis consists of tw unique pharmacophoric features based on structural requirements for biological activit hydrogen bond donor (D) and ring aromaticity (R). As shown in Figure 3, each feature this pharmacophoric map signifies a common particular portion associated with the s lected compounds.

Pharmacophore Hypothesis Generation
Using the PHASE module's "Developing a pharmacophore model" [25,26], 20 h potheses were created to explain how active molecules bind to receptors with a box siz of 1 Å and a minimum inter site distance of 2 Å. A pharmacophore site is a characterist feature of a conformation that has been mapped to a specific location. Common pharm cophoric features were identified from a set of variants-feature types that identify a pu tative pharmacophore. Using a scoring function as a survival score, the common pharm cophore hypotheses were investigated. The site, vector, volume, selectivity scores, an number of matches were calculated for each of the generated hypotheses. The PHAS

Pharmacophore Hypothesis Generation
Using the PHASE module's "Developing a pharmacophore model" [25,26], 20 hypotheses were created to explain how active molecules bind to receptors with a box size of 1 Å and a minimum inter site distance of 2 Å. A pharmacophore site is a characteristic feature of a conformation that has been mapped to a specific location. Common pharmacophoric features were identified from a set of variants-feature types that identify a putative pharmacophore. Using a scoring function as a survival score, the common pharmacophore hypotheses were investigated. The site, vector, volume, selectivity scores, and number of matches were calculated for each of the generated hypotheses. The PHASE module has six in-built pharmacophore features, including hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), ring aromaticity (R), positively ionizable (P), and negatively ionizable (N) groups (Supplementary File; Table S1).

An Atom Based 3D-QSAR
Schrodinger Maestro software was used to build an atom-based 3D-QSAR model using a collection of aligned ligands in order to predict activities for additional compounds [27]. The training set consisted of 70% of the dataset's compounds, while the test set consisted of 30%. The compounds were clustered using a PLS factor of 4. The Atom type fraction segment demonstrated the fraction owing to each atom type in the QSAR model for each number of PLS factors included in the model. Furthermore, actual activity versus predicted activity for compounds in the training as well as the test sets were plotted to form a scatter plot.

Virtual Screening
The virtual screening study was carried out using the ZINC data base, and the AAHRR 1 hypothesis was employed to screen ZINC compounds using the Lipinski rule of 5 (Supplementary File; Table S2) [28]. Using different filters, a total of 3568 molecules were screened. The GLIDE module of Schrodinger was used to screen these compounds using different docking methodologies. The Swiss Target Prediction tool was used to predict the target of these selected compounds. This tool predicts all the available targets of the molecule. Different targets were identified in the current investigation, with CA being the most important one for the particular molecule.

Molecular Docking
The docking studies were carried out using the Schrodinger Maestro software's Glide module [29]. The score function in the software was used to rank and group distinct possible adduct structures generated by molecular docking [30] (Supplementary File; Table S3). Based on the binding properties of the ligand and target, it predicts the three-dimensional structure of any complex. The method of predicting ligand conformation and orientation (or posing) within a specific binding site is referred to as docking. In Maestro wizard, the protein structure was pre-processed using the "protein preparation wizard". Hydrogen atoms and certain essential bonds were introduced to the missing site of the protein molecule by automatically generating states and refinement steps phases of the module. After the optimization process, receptor grid generation was processed and docking scores were analysed with different docked ligand conformations [31,32].

ADME Properties Prediction
The drug-like activity of the ligand molecule was predicted using the Schrodinger software's QikProp module [33]. It predicts both physicochemical descriptors and pharmacokinetic properties with the objective of increasing the success rate of compounds for further development. The ADME study suggests that drug-like properties of ligand molecules were validated using Lipinski's rule of five [34,35]. The rule states that the molecular weight <500, Log P o/w <5, hydrogen bond donor ≤5, and hydrogen bond acceptor ≤10. The compounds that meet these criteria are classified as drug-like compounds (Supplementary File; Tables S4 and S5).

Pharmacophore Mapping: Selection of the Best Pharmacophore Hypothesis
All the selected compounds from the database were used to generate pharmacophoric hypothesis or the minimum features of the compound required to bind with receptor. Among 20 pharmacophore hypotheses, we screened best pharmacophoric hypotheses on the basis of various scores described in Table 2. Three hydrogen bond donors (D) and two ring aromaticity (R) were identified as common pharmacophore features. DDDRR 1 and DDDRR 2 were found to be the best among the 20 hypotheses generated (Supplementary File; Table S1) by the PHASE module, based on a scoring function of each parameter value listed in Table 2 [36]. The "survival" scoring (S) function was responsible for determination of features from described models and to rank all hypotheses. The scoring algorithm involves selectivity, number of ligands matched, relative conformational energy, and activity. However, the models should also differentiate between inactive and active molecules. If an inactive molecule scores well, the hypothesis could be invalid as it does not discriminate between active and inactive molecules. For this reason, adjusted survival score S_I was calculated by subtracting the inactive score from the survival score.

Selection of Atom Based QSAR Model
The QSAR findings for the hypothesis are shown in the statistical table obtained for the training and test set molecules. Several statistical parameters were employed to examine the robustness of the QSAR model, including SD, R 2 , F, P, RMSE, Q 2 , and Pearson-R [37]. High R 2 (higher than 0.6), Q 2 (greater than 0.5), Pearson-R (greater than 0.5), and F values characterize a good QSAR model. All the four models generated by the module are shown in Table 3. Based on the above information, the fourth model was chosen as the best QSAR model due to higher Q 2 and R 2 values as 0.7888 and 0.8757, respectively.
To demonstrate the uniform distribution of training set molecules over the straight line passing through the origin (0, 0), a scattered plot was drawn between experimental activity and predicted activity of training and test set compounds (Figure 4).

Evaluation of Contour Map
The effects of different substituents on biological activity were determined by the study of contour maps, which are shown in Figure 5A-E. It also examines the variation of substituents and their biological activity. The contour map is represented by color coding, including different substituents on the core moiety. This color coding has a specific indication of the substituent and is useful in the expansion of novel compounds as CA inhibitors. The blue color showed enhanced activity, whereas the red color indicated decrease in activity. The electron-withdrawing group substitution on the phenyl ring showed a decrease in activity. Hydrogen bond donor group substitution on phenyl ring with para substitution showed an increase in activity. Hydrophobic substitution at the ortho position displayed good activity. Negative ionic group showed good activity at the para position and deceased activity in the meta position. The positive ionic group substitution in the phenyl ring showed decrease in activity.

Evaluation of Contour Map
The effects of different substituents on biological activity were determined by the study of contour maps, which are shown in Figure 5A-E. It also examines the variation of substituents and their biological activity. The contour map is represented by color coding, including different substituents on the core moiety. This color coding has a specific indication of the substituent and is useful in the expansion of novel compounds as CA inhibitors. The blue color showed enhanced activity, whereas the red color indicated decrease in activity. The electron-withdrawing group substitution on the phenyl ring showed a decrease in activity. Hydrogen bond donor group substitution on phenyl ring with para substitution showed an increase in activity. Hydrophobic substitution at the ortho position displayed good activity. Negative ionic group showed good activity at the para position and deceased activity in the meta position. The positive ionic group substitution in the phenyl ring showed decrease in activity.

Molecular Docking Analysis
The Glide module was utilized to conduct molecular docking between the potent derivatives and the target protein [38,39]. The potent analogues 4m, 4o, 4s, 4p, and Table 4. The purple arrows indicate hydrogen bond interactions, whereas the green arrows indicate π-π-stacking interactions, as seen in Figure 6 for compound 4m, Figure 7A for compound 4o, and Figure 7B for compound 4s [40,41].

Molecular Docking Analysis
The Glide module was utilized to conduct molecular docking between the potent derivatives and the target protein [38,39]. The potent analogues 4m, 4o, 4s, 4p, and Table 4. The purple arrows indicate hydrogen bond interactions, whereas the green arrows indicate π-π-stacking interactions, as seen in Figure 6 for compound 4m, Figure 7A for compound 4o, and Figure 7B for compound 4s [40,41].

Virtual Screening
The molecular docking methodologies such as HTVS, SP, and XP were applied to screen the compounds from the ZINC database. Each step screened 20% of the highestscoring compounds with high docking scores. The top compounds, ZINC77699643, ZINC89275054, ZINC77671412, and ZINC70762033, were selected with xp docking scores of −6.178, −5.743, −5.561, and −5.535, respectively, out of 333 compounds screened using the SP docking approach. These compounds were chosen as the final ZINC compounds for further investigation and assessed using MMGBSA to calculate binding interaction energy. The compound ZINC77699643 interacted with Pro201 and Thr199 amino acid residues, as shown in Figure 8A. These interactions demonstrated the significance of CA activity.
The compound ZINC89275054 showed binding interactions with amino acid residues His94, Gln92, and Thr199 in the same cavity, as shown by a crystal ligand ( Figure 8B). The compound ZINC77671412 demonstrated binding interactions with several amino acids at the receptor's binding region, including His200 and Gln92, which are critical for activity ( Figure 9A). The compound ZINC70762033 showed binding interaction with different amino acids such as His200, His94, Gln92, and Thr199 at the binding site of receptor ( Figure 9B).

ADME Properties Prediction
The ADME properties of all the compounds from the dataset were predicted through the QikProp module of Schrodinger; the results are given in Table 5. The software predicted the compounds' physicochemical properties, lipophilicity, drug-like behaviour, water solubility, permeability through BBB, pharmacokinetics, and synthetic accessibility [42,43]. All compounds have a molecular weight ranging from 259.34 to 399.24, number of hydrogen bond acceptors ranging from 4.5 to 7.0, and number of hydrogen bond donors ranging from 3 to 4. All the predicted ADME properties satisfied the Lipinski's rule of five for all the compounds.

Optimization of Novel Ligands
Optimization of ligands revealed that the substitution of the hydrogen bond donor group on the phenyl ring at the para position, hydrophobic group at the ortho position, and negative ionic group at the para position showed increase in activity. However, substitution of the negative ionic group at the meta position of the phenyl ring, and positive ionic group substitution on the phenyl ring displayed a decrease in activity ( Figure 10). substitution of the negative ionic group at the meta position of the phenyl ring, and pos tive ionic group substitution on the phenyl ring displayed a decrease in activity (Figu 10).

Conclusions
The computational design of 1,3,4-thiadiazole derivatives as potential CA inhibito was discussed in this study. The PHASE module (Schrodinger) was used to investiga the pharmacophore model on 27 compounds taken from a dataset. The compound 4m [( (2-(dimethylamino)benzylidene)-N-(5-thioxo-4,5-dihydro-1,3,4-thiadiazol-2-yl)hydrazine-1-carboxamide)] showed the best SP docking score of −5.217. The predicted ADM properties revealed that all the compounds followed the Lipinski's rule of five, indicatin their drug-like potential. Optimization of ligands can help to generate novel thiadiazo derivatives by using the response of a contour map. The generated contour map showe the substitution of the electron-withdrawing group on the phenyl ring responsible for d crease in CA activity. However, the hydrophobic group substitution at the ortho positio and hydrogen bond donor substitution at the para position of the phenyl ring showed a increase in CA activity. These results may aid in the development of novel thiadiazo compounds as CA inhibitors, which might have a major influence on the treatment oxidative stress induced by low oxygen pressure at high altitudes.
Supplementary Materials: The following supporting information can be downloade at:www.mdpi.com/xxx/s1, Table S1: The pharmacophore study carried out for the dataset; Table S The HTVS docking study carried out for the ZINC compounds; Table S3: The docking study carrie out for the dataset; Table S4: The QIKPROP study carried out the ZINC database; Table S5: Th ADME study carried out for the database.

Conclusions
The computational design of 1,3,4-thiadiazole derivatives as potential CA inhibitors was discussed in this study. The PHASE module (Schrodinger) was used to investigate the pharmacophore model on 27 compounds taken from a dataset. The compound 4m [(2-(2-(dimethylamino)benzylidene)-N-(5-thioxo-4,5-dihydro-1,3,4-thiadiazol-2-yl)hydrazine-1carboxamide)] showed the best SP docking score of −5.217. The predicted ADME properties revealed that all the compounds followed the Lipinski's rule of five, indicating their druglike potential. Optimization of ligands can help to generate novel thiadiazole derivatives by using the response of a contour map. The generated contour map showed the substitution of the electron-withdrawing group on the phenyl ring responsible for decrease in CA activity. However, the hydrophobic group substitution at the ortho position, and hydrogen bond donor substitution at the para position of the phenyl ring showed an increase in CA activity. These results may aid in the development of novel thiadiazole compounds as CA inhibitors, which might have a major influence on the treatment of oxidative stress induced by low oxygen pressure at high altitudes.
Supplementary Materials: The following supporting information can be downloaded at:https:// www.mdpi.com/article/10.3390/cimb44030068/s1, Table S1: The pharmacophore study carried out for the dataset; Table S2: The HTVS docking study carried out for the ZINC compounds; Table S3: The docking study carried out for the dataset; Table S4: The QIKPROP study carried out the ZINC database; Table S5: The ADME study carried out for the database.