Journal of Pharmacy And Bioallied Sciences

: 2021  |  Volume : 13  |  Issue : 4  |  Page : 387--393

Development of Hif1α pharmacogenomic mutation models to study individual variations in drug action for tumor hypoxia: An in silico approach

Vaisali Balasubramaniam, P K Krishnan Namboori 
 Computational Chemistry Group (CCG), AMMAS Research Lab, Amrita School of Engineering, Amrita Vishwa Vidyapeetham, Coimbatore, Tamil Nadu, India

Correspondence Address:
Dr. P K Krishnan Namboori
AMMAS Research Lab, Amrita School of Engineering, Amrita Vishwa Vidyapeetham, Ettimadai, Coimbatore - 641 112, Tamil Nadu


Objective: Tumor hypoxia, a predominant feature of solid tumor produces drug resistance that significantly impacts a patient's clinical outcomes. Hypoxia-inducible factor 1-alpha (HIF1α) is the major mutation involved in establishing the microenvironment. As a consequence of its involvement in pathways that enable rapid tumor growth, it creates resistance to chemotherapeutic treatments. The propensity of medications to demonstrate drug action often diverges according to the genetic composition. The aim of this study is therefore to examine the effect of population-dependent drug response variations using mutation models. Methods: Genetic variations distinctive to major super-populations were identified, and the mutated gene was acquired as a result of incorporating the variants. The mutated gene sequence was transcribed and translated to obtain the target amino acid sequence. To investigate the effects of mutations, protein models were developed using homology modeling. The target templates for the backbone structure were identified by characterization of primary and secondary protein structures. The modeled proteins were then validated for structural confirmation and flexibility. Potential models were used for interaction studies with hypoxia-specific molecules (tirapazamine, apaziquone, and ENMD) using docking analysis. To verify their stability under pre-defined dynamic conditions, the complexes were subjected to molecular dynamics simulation. Results: The current research models demonstrate with the pharmacogenomic-based mutation of HIF1α the impact of individual variants in altering the person-specific drug response under tumor hypoxic conditions. It also elucidates that the therapeutic effect is altered concerning population-dependent genetic changes in the individual. Conclusion: The study, therefore, asserts the need to set up a personalized drug design approach to enhance tumor hypoxia treatment efficacy.

How to cite this article:
Balasubramaniam V, Namboori P K. Development of Hif1α pharmacogenomic mutation models to study individual variations in drug action for tumor hypoxia: An in silico approach.J Pharm Bioall Sci 2021;13:387-393

How to cite this URL:
Balasubramaniam V, Namboori P K. Development of Hif1α pharmacogenomic mutation models to study individual variations in drug action for tumor hypoxia: An in silico approach. J Pharm Bioall Sci [serial online] 2021 [cited 2022 Jun 29 ];13:387-393
Available from:

Full Text


It is critical to maintain homeostasis between the supply and consumption of oxygen in a cell. The inadequacy of oxygen compared with the normal range required for physiological functions is known as hypoxia. Hypoxic regions occur as a modification to survive the wide-varying oxygen levels in rapidly growing tumors.[1] Hypoxia-inducible factor 1A (HIF1A) expression has been found to be more significant in establishing resistance to therapeutic interventions in solid tumors. Under a mutated condition, the HIF1A pathway is altered for cancer cell survival by promoting drug resistance with the microenvironment. HIF1A potentially supports tumor growth by modifying cellular functions, stimulating angiogenesis, promoting cellular proliferation, and increasing treatment resistance.[2]

The adaptation exhibits certain variations that could be considered as the genetic signatures. Identifying the biomarkers of the mutation would help to stratify the patients based on their proneness to different risk categories. Single-nucleotide polymorphism (SNP) has been considered to be the influential genetic marker for understanding disease susceptibility.[3] Based on the individual variants, the microenvironment varies for each individual, and it is reflected on the variation in drug response.[4]

There is evidence of increased reactive oxygen species (ROS) in the tumorous environment, and the drug delivery systems focus on the over-production of ROS. Because of the increased hypoxia, signaling is impaired, making it difficult for the drug delivery system to recognize the location. Meanwhile, the tumor hypoxic microenvironment promotes cell proliferation and tumor growth.[5] Despite the fact that many pro-drugs were developed to treat this condition, they were found to be less effective because of a lack of molecular profiling. Only if the mutation and its drug interaction variability are extensively studied using in silico studies can a molecular level interpretation be obtained.[4],[6] The Food and Drug Administration clinical data report that the drug responses of tirapazamine, apaziquone, and ENMD differ in terms of the demographic variability of the patients enrolled in the study.[7],[8],[9] Previous research in this field of study addresses other strategies of treating the issue from a therapeutic perspective. The current study, on the other hand, focuses on the molecular level understanding of the variations and their impact on drug action changes, with the novel idea of designing mutation models of demographically varying targets corresponding to the five distinct super-populations (American, East Asian, South Asian, European, and African) based upon the 1000 genome project.

 Materials and Methods

Clinical data and mutation analysis

The clinical trial data have been analyzed from the medical study result repository maintained by the National Institutes of Health (NIH), which helps to know the effectiveness of the drug under study. The major mutations involved in the signaling and the genetic variants corresponding to susceptibility of the genes toward mutation have been studied from the National Centre for Biotechnology Information (NCBI). Based upon the functional annotations, SNP Nexus has been used to identify the potentially significant SNPs present in the super-population.[10]

Target modeling

In order to obtain the three-dimensional (3D) target structure that comes as an outcome of the mutated gene, protein homology modeling is adapted. The mutated gene sequence has been transcribed and translated to obtain the target amino acid sequence using the Expert Protein Analysis System (ExPASy) Translate tool.[11] For the template, proteins given as a backbone for the modeling process are chosen based upon characterization of the protein.[6] The gene-specific proteins were retrieved from the archive of the Protein Data Bank. The primary structure characterization has been performed using the integrative portal of the bioinformatics resource, the ExPASy ProtParam tool.[12] The secondary structure characterization was studied using the self-optimized prediction method (SOPMA).[13] The sub-cellular location of a protein was found using the support vector machine-based predictor Bacello (balanced sub-cellular localization predictor).[14] The target 3D structures were built using an automated structural bioinformatics server named SWISS Model.[15]

Evaluation of the mutation model

The modeled proteins were then validated with their Q Mean Z score, Errat Quality factor, Overall G factor, and Ramachandran plot using SAVES v6.0.[16] The 3D structure of the protein was analyzed for potential errors using the Protein Structure Analysis tool, ProSA.[17] Monte Carlo dynamics was adapted to study the flexibility and fluctuation of the protein structure when placed in a dynamic environment using CABS-flex 2.0.[18]

Identification of active sites

In order to facilitate the interaction, the protein surface had been characterized and the functional sites have been found using Computed Atlas of Surface Topography of proteins (CASTp).[19]

Interaction analysis with molecular docking

The docking study has been carried out to analyze the interaction between the ligands and designed macromolecules using the Autodock 4.2.6 suite.[20] The electrostatic potential of the protein was provided by applying the Kollman charges, and they were assigned AD4 type. The ligand was equilibrated for electronegativity using the Gasteiger charges. The grid along the x-, y-, and z-axes was defined as the rigid frame of active sites, and docking was carried out using the Lamarckian Genetic Algorithm. The best binding confirmation was visualized to see its interaction sites and bonds using BIOVIA Discovery Studio Visualizer.

Molecular dynamics simulation of the complex

The molecular dynamics (MD) simulation of the complex was carried out using Nanoscale Molecular Dynamics (NAMD) with Visual Molecular Dynamics (VMD) v1.9.3.[21] For high performance accuracy, the Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field has been applied using the CHARMM-GUI simulation system. The topology and parameter files of the ligand with its coordinates have been obtained using the CHARMM General Force Field (CGenFF) application. To the merged protein structure files (PSFs) of ligands and proteins, the solvation box is added by setting the minimum and maximum (0–5) box padding values and solvated. With the set minmax values, the box cell vector and the periodic boundary conditions are generated. The MD simulation is then carried out for the set parameters under a 310 K temperature, which matches the human body temperature, for a time period of 100 ns. The schematic representation of the methodology is given in [Figure 1].{Figure 1}

 Results and Discussion

Protein characterization

The drug action-specific HIF1A gene-associated proteins, namely, 1H2K, 1L3E, 1LM8, 1LQB, 2ILM, 3HQR, 4AJY, 4H6J, 5JWP, 5L9B, and 5LA9, were taken for the characterization. Proteins with a higher molecular weight (Mw) must be chosen because smaller protein templates are incapable of obtaining a good backbone structure and forming a stable target. The isoelectric point (pI) denotes both the cationic nature and anionic nature of a protein; hence, a protein with pI toward an acidic pH (<7) would be ideal for studying cancer environments. With a half-life of 30 hours or more, the protein structure is a good target for the ligand to bind to for a longer period of time and provide drug action. An instability index ≤40 is found to be ideal for a stable protein, so the proteins 3HQR, 4H6J, 5L9B, and 5LA9 have evolved toward stability. Proteins with long aliphatic chain groups have a higher index ranging from 66.5 to 84.33, indicating greater thermal stability; with the exception of the 1L3E protein, all were found to be thermally stable. The GRAVI value indicates the measure of the hydrophobic nature of the protein. The negative values of the chosen proteins denote that they have hydrophilic nature. The protein solubility plays a major role in choosing the protein as a template; hence, a moderately soluble protein has to be chosen [Table 1]. From the secondary structure characterization, the alpha helix structure denotes the percentage of protein structure that is nascent or non-functional in nature. The beta turn elucidates the portion of protein that is evolved and has got functional side chains. The random coil percentage is the proportion of the protein that is evolved and functional in nature for the interaction to happen [Table 2]. Therefore, the protein 4h6j, which is evolved and has more functional sites and satisfies all the characteristic needs, has been chosen as the template.{Table 1}{Table 2}

Evaluation of the protein models

The designed protein models have been evaluated using the Q Mean Z score, Errat Quality factor, Overall G factor, and Ramachandran plot and through MD simulation [Table 3]. The protein models belonging to the five different populations have been obtained and assessed. The Q Mean Z score helps to know the degree of nativeness of the model with the naturally available proteins. An Z score of –4 or below is considered to be an unfit model, and an Z score above –4 shows structural similarity with the proven protein structures. The Errat Quality factor depicts the overall quality of the model, and the optimum range of it is above 50. The Overall G factor describes about the steriochemical property of the model; a value of –0.5 or above is said to be optimum to satisfy the geometrical substantiality of the model. The Ramachandran plot helps to know how much proportion of the structure falls to the less steric hindrance region. The protein stability is further validated using the MD simulation, and a root-mean-square deviation (RMSD) value below 3.5 Å is found to be the ideal difference between the initial and final structures. As a result, the models were assessed, and all five population models were found to satisfy the evaluation criterion, making them potential candidates for the interaction analysis.{Table 3}

Interaction study and MD simulation of the complex

The molecular docking results reveal the best conformational interaction pose of the protein model and the ligand molecules. The best docking complex is the one with minimum free binding energy obtained when the ligands bind to the protein in the hotspot. The free binding energy is obtained as a sum of the interaction energy and steric hindrance. The drug action is observed for a well-bound complex depicted with a lower overall energy in kilocalories per mole. From the interaction study of ENMD [Figure 2], it is evident that the drug shows more drug action specific to the American population with a binding energy of –7.26 kcal/mol. It suits the least for the European population with a binding free energy of –6.25 kcal/mol. The apaziquone interaction [Figure 3] shows different drug actions for different populations, and it is also comparatively found to be more suitable for the American population with a free binding energy of –5.91 kcal/mol. It shows less drug action for the European population with a free binding energy of –5.40 kcal/mol. The tirapazamine interaction [Figure 4] shows different drug actions for the same set of populations as observed in the clinical trial study and yet again was found to be most suitable for the American population with a free binding energy of –5.65 kcal/mol [Table 4]. Tirapazamine shows less drug action for the East Asian population with a free binding energy of –5.39 kcal/mol. Comparing the results, the variations involved in producing the population-specific mutation affect the disease setup, which is expressed in the change in drug action.[22]{Figure 2}{Figure 3}{Figure 4}{Table 4}

The target-ligand complexes were subjected to MD simulation to observe their stability in a dynamic environment. Their time-evolved behavior has been studied by analyzing the RMSD from the initial and final docked confirmations under the set parametric conditions. A deviation of 4 Å or below is ideal for the complex to be stable in a human environment and sustain the biological dynamics. Thus, all the complexes are found to fall within the optimum range, proving them to be potential complexes to elucidate variation in drug action with respect to population changes [Table 4].


Several diagnostic and therapeutic measures are suggested toward treatment of tumor hypoxia. Because of the influence of the population-specific genetic variants, the variation in drug action and side effects are still observed. The proposed work investigates the alteration in drug response observed in the clinical trial analysis using pharmacogenomic-based mutation models. The molecular docking followed by MD analysis helps to visualize the person-specific drug response altering with respect to the population class. The results clearly show that with altering individual variations, the mutated target is unique and that the same drug shows altering drug action. Hence, following a general drug for the condition for all populations would not be suitable as an efficient treatment procedure. The study recommends that a customized approach to drug design is required to improve therapeutic efficacy and safety. To provide successful therapeutic care with the right drug at the right time, a precise molecular understanding is required and the associated personalized drug design technique must be followed.

Financial support and sponsorship


Conflicts of interest

There are no conflicts of interest.


1Al Tameemi W, Dale T, Al-Jumaily R, Forsyth N. Hypoxia-modified cancer cell metabolism. Front Cell Dev Biol 2019;7. doi: 10.3389/fcell. 2019.00004.
2Zhang Z, Niu N, Gao X, Han F, Chen Z, Li S, Li J. A new drug carrier with oxygen generation function for modulating tumor hypoxia microenvironment in cancer chemotherapy. Colloids Surfaces B Biointerfaces 2019;173:335-45.
3Anand CL, Vyshnavi AH, Francis G, Karthikeyan S, Kumar PS, Namboori PK. Population wise variation of breast and ovarian cancer-A pharmacogenomic approach. Mater Today Proc 2018;5:16106-10.
4Spiegelberg L, Houben R, Niemans R, de Ruysscher D, Yaromina A, Theys J, et al. Hypoxia-activated prodrugs and (lack of) clinical progress: The need for hypoxia-based biomarker patient selection in phase III clinical trials. Clin Transl Radiat Oncol 2019;15:62-9.
5Lv Y, Zhao S, Han J, Zheng L, Yang Z, Zhao L. Hypoxia-inducible factor-1α induces multidrug resistance protein in colon cancer. OncoTargets Ther 2015;8:1941.
6Siniprasad P, Nair B, Balasubramaniam V, Sadanandan P, Namboori PK, Nath LR. Evaluation of kaempferol as AKT dependent mTOR regulator via targeting FKBP-12 in hepatocellular carcinoma: An in silico approach. Lett Drug Des Discov 2020;17:1401-8.
7A Study of ENMD-2076 in Ovarian Clear Cell Cancers - 2021. Available from: [Last accessed on 2022 Jan 31].
8Tirapazamine Combined With Chemo and RT in Limited-Stage Small Cell Lung Cancer - 2021. Available: [Last accessed on 2022 Jan 31].
9Single Dose Intravesical Apaziquone Postoperative in Patients Undergoing TURBT for Noninvasive Bladder Cancer. 2021. Available from: [Last accessed on 2022 Jan 31].
10Vaisali B, Parvathy CR, Hima AM, Krishnan PK. Tumor hypoxia Diagnosis using deep CNN learning strategy a theranostic pharmacogenomic approach. Int J Prognostics Health Management 2019; 10.
11ExPASy - Translate tool. 2021. Available from:
12“ExPASy - ProtParam tool”,, 2021. Available from: [Last accessed on 2022 Jan 31].
13I. UCBL, “[email protected]: SOPMA secondary structure prediction. 2021. Available from: [Last accessed on 2022 Jan 31].
14Pierleoni A, Martelli PL, Fariselli P, Casadio R. BaCelLo: A balanced subcellular localization predictor. Bioinformatics 2006;22:e408-16.
15Sangeetha M, Saranya TS, Sathianarayanan S, AM HV, Namboori PK. Design and development of potential flavonoid moiety for Pbp2a inhibition for Mrsa therapy-A computational technique. Biomed Pharmacol J 2020;13:687-92.
16“SAVESv6.0-Structure Validation Server”,, 2022. Available from: [Last accessed on 2022 Jan 31].
17Wiederstein M, Sippl M. ProSA-web: Interactive web service for the recognition of errors in three-dimensional structures of proteins. Nucleic Acids Res 2007;35:W407-10.
18Kuriata A, Gierut AM, Oleniecki T, Ciemny MP, Kolinski A, Kurcinski M, et al. CABS-flex 2.0: A web server for fast simulations of flexibility of protein structures. Nucleic Acids Res 2018;46:W338-43.
19Binkowski TA, Naghibzadeh S, Liang J. CASTp: Computed atlas of surface topography of proteins. Nucleic Acids Res 2003;31:3352-5.
20Forli S, Huey R, Pique M, Sanner M, Goodsell D, Olson A. Computational protein–ligand docking and virtual drug screening with the AutoDock suite. Nat Protoc 2016;11:905-19.
21Available from: 2016. [Last accessed on 2021 May 04].
22Hicks N, Giffen SR, Culviner PH, Chao MC, Dulberger CL, Liu Q, et al. Mutations in dnaA and a cryptic interaction site increase drug resistance in mycobacterium tuberculosis. PLOS Pathog 2020;16:e1009063.