Computational analysis of missense variant CYP4F2*3 (V433M) in association with human CYP4F2 dysfunction: a functional and structural impact

Data collection

The FASTA sequence of the human CYP4F2 protein was obtained from UniProt (ID: P78329) [28]. The clinical information of variant CYP4F2*3 and its relation to warfarin dose was retrieved from the PharmGKB (ID: PA166169424) [29] and PubMed databases [8–10].

Structural modeling and verification

Because the three-dimensional structure of human CYP4F2 does not exist, homology modeling was performed using the I-TASSER web server for both native and mutant form [30]. Then, the quality of the modeled structure was evaluated using Procheck, ERRAT, and Verify3D in the SAVES web server (https://saves.mbi.ucla.edu/). The Procheck platform checks the stereochemical quality of the model using a Ramachandran plot [31]. A protein structure with more than 80% amino acids in the core and allowed regions indicated high quality. The ERRAT webserver statistically analyzed the non-bonded interactions between different atom types in the protein structure [32]. Verify3D webserver has characterized the compatibility of an atomic model (3D) with its amino acid sequence (1D) by analyzing a structural class according to its location and environment (alpha, beta, loop, polar, nonpolar, etc.) [33]. A protein structure with a Verify3D value greater than 80% is defined as having good compatibility with 1D and 3D protein structures and good quality. ProsA is a web server program that is frequently used to check 3D models of protein structures for possible errors. The overall quality score computed by ProSA for a particular protein structure is represented in a plot that indicates the scores (Z-score) of all experimentally characterized protein chains recently available in the Protein Data Bank (PDB) [34]. The protein structure Z-score, which is located in the experimentally characterized protein structure range, indicates a high-quality protein structure.

Protein stability prediction

Five different prediction tools, I-mutant2, MUpro, DUET, SDM, and mCSM were utilized to identify the structural impact of the CYP4F2*3 (V433M) variant on the CYP4F2 protein. The server of I-Mutant 2.0 (http://folding.biofold.org/i-mutant/i-mutant2.0.html) estimates the impact of amino acid substitution on protein stability via a support vector machine-based method. This program was trained on ProTherm, the most comprehensive dataset of experimental information on protein mutations [35]. The MUpro tool (https://mupro.proteomics.ics.uci.edu/) uses neural networks and support vector machine methods to assess any alteration in the stability of the protein with approximately 84% accuracy [36]. Three structure-based algorithms, DUET, SDM, and mCSM (http://biosig.unimelb.edu.au/duet/) are closely intertwined algorithms and their input can be supplied together in the form of a native protein structure, while separately bringing up the position of relevant SNP [37].

Protein property and evolutionary conservation analysis

GETAREA is a structure-based tool (https://curie.utmb.edu/getarea.html) used to determine the effect of the V433M variant on the solvent accessibility surface area (SASA) of the CYP4F2 protein. The percentage of side-chain surface area to random coil value per residue (% SASA) estimates whether the residue is buried in the solvent or exposed [38]. The server of Project HOPE (http://www.cmbi.ru.nl/hope/input/) by utilizing the integrating information from the UniProt database and Distributed Annotation System (DAS) servers, predicts the effect of amino acid substitution on the protein structure [39]. To visualize the impact of the V433M variant at the structural level, Chimera 1.13 software was used. This software is an extended system for the interactive visualization and evaluation of molecular structures and related data [40]. The ConSurf web server (http://consurf.tau.ac.il) was used to define the conservation sites in the protein compared to its homologs during evolution through conservation scores ranging from 1 to 9, which are classified as variable, intermediate, and conserved [41].

Functional impact prediction

To determine the functional consequences of the V433M variant on the CYP4F2 protein, five different tools including PolyPhen-2, FATHMM-MKl, PANTHER-PSEP, CADD, and PhD-SNP^g were employed. The server of PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/) determines the possible effect of amino acid substitution on protein properties, based on position-specific independent count (PSIC) scores [42]. The difference between the scores categorizes the variants as probably damaging, possibly damaging, or benign. FATHMM-MKL (https://fathmm.biocompute.org.uk/) is a web-based tool in which a prediction score (p-value) above 0.5 identifies the variant that has a harmful effect on protein function [43]. Based on the substitution position-specific evolutionary conservation (subPSEC) score, PANTHER (http://pantherdb.org/) predicts the probability of a variant generating a harmful effect on the protein. If the score of subPSEC is ≥ 0.5, the variant is deemed deleterious [44]. CADD (http://cadd.gs.washington.edu/) as a linear kernel support vector machine-based, predicts functionally important variants using a Phred-scale C score. CADD considers amino acid substitutions as the most deleterious when they have a C-score of 20 or greater (≥ C20) [45]. PhD-SNP^g (http://snps.biofold.org/phd-snpg./) is a machine learning-based server that identifies the impact of variants in both non-coding and coding regions using evolutionary data [46].

Molecular dynamic (MD) simulation analysis

Gromacs 2019 was applied to run an MD simulation with an optimized potential for liquid simulations (OPLS) force field. Simple point charge (SPC) water was added to the simulation box cubic with a 10 A° [47, 48]. The simulation systems’ charges were neutralized using Na⁺ and Cl⁻ ions. To minimize the simulation system, a steepest descent minimization integrator with a maximum force of less than 100 KJ.mol^− 1. nm^− 1 and 5000 minimization steps were used. The simulation system was then equilibrated by a constant number of particles, volume, and temperature (NVT) and a constant number of particles, pressure, and temperature (NPT) ensemble. For NVT equilibration through 100 picoseconds, the leap-frog integrator and linear constraint solver (LINCS) algorithm were used [49]. The PME (particle mesh Ewald) algorithm with a grid spacing of 0.16 nm and a radius cutoff of 1.0 nm was considered for calculating electrostatic interaction. The equilibrated system was put into an MD simulation with a time step of 2 femtoseconds (fs) for 200 nanoseconds (ns) of simulation. The output trajectories were analyzed using root-mean-square deviation (RMSD), radius of gyration (Rg), solvent access surface area (SASA), total number of intramolecular hydrogen bonds, root-mean-square fluctuation (RMSF), secondary structure by DSSP, principal component analysis (PCA) and free energy landscapes (FELs) to determine the stability and protein structure changes during simulation.

Molecular docking analysis

As VK1 is a substrate for the CYP4F2 protein and its effective interaction is essential for the proper function of CYP4F2, molecular docking calculations were performed using the HADDOCK web server to determine the effect of the CYP4F2*3 variant on the native protein structure [50]. This web server uses a flexible docking approach to define the best protein-ligand interaction mode. In addition, the PRODIGY web server was applied to compute the binding affinity of the protein-small molecule [51] and the bond-type analyses were done by Discovery Studio 4.5 [52].

Structural modeling validation and quality estimation

Owing to the lack of tertiary structure of CYP4F2, homology modeling was applied using the I-TASSER webserver (Fig. 1A). The quality of the predicted structure was analyzed using different web servers (Table 1). Checking the quality of a model by Ramachandran plot with Procheck web server illustrated that more than 80% of the protein structures were located in the core and allowed regions, indicating the high physicochemical quality of the protein structure. A Verify3D value of more than 80% showed high compatibility between the primary and tertiary structures of CYP4F2. The ERRAT rate of 91.01 is interpreted good overall quality of protein structure. The Z-score of -8.66 was located in regions of proteins that were identified experimentally, indicating the high quality of the predicted model of CYP4F2 (Fig. 1B). The position of wild-type residue (V433) and mutant-type residue (M433) of V433M variant in CYP4F2 protein was shown in Fig. 1C and D.

Protein stability prediction

I-mutant2, MUpro, DUET, SDM and mCSM tools predict protein stability changes in the CYP4F2*3 (V433M) variant based on the ΔΔG value (change in Gibbs free energy). The ΔΔG value is a metric used to calculate the impact of a single nucleotide variant on protein stability. A ΔΔG score less than zero indicated lower stability. From the results illustrated in Table 2, the CYP4F2*3 (V433M) variant exhibited reduced protein stability. Destabilization of a protein structure can change its biological function and disrupt the signaling cascades and normal pathways of the protein [53].

Protein property and evolutionary conservation analysis

GETAREA calculates the SASA of individual residues. The threshold is described as %SASA > 50 and %SASA < 20, indicating exposed and buried residues to the solvent, respectively [38]. Based on the results obtained from the GETAREA online for the V433M variant, both wild-type (V433) and mutant-type (M433) amino acids were buried or located in the interior of the CYP4F2 protein. However, the mutated residue had a lower % SASA value (% 5.9) than the wild-type V433 (% 9). This result was compatible with the deep view of the locations of the wild-type (V433) and mutant-type (M433) amino acids (Fig. 1D), which showed that M433 residue had more buried side chain than V433 residue.

The Project HOPE server investigated the effect of the V433M variant on the physiochemical characteristics, intermolecular interactions, and functional and structural properties. It was predicted that the wild-type residue (V433) would be replaced by a mutant-type residue (M433), which has a larger side chain than the V433 moiety and may not fit into the new position. For physical and structural reasons, if a substituted amino acid does not fit into a protein, it leads to structural changes that are harmful [54]. Furthermore, the Project HOPE server predicted that the V433 residue is located near a highly conserved region and in contact with residues in a domain that is important for protein activity.

Structural analysis of neighboring amino acids around the V433 residue in less than 5-angstrom regions using the Chimera software revealed that amino acids including H103, P104, T427, N430, P431, A432, W434, P435, R453, S454 and A457 were located close to the V433 residue. For the V-to-M conversion at position 433, the mutated residue with a larger side chain than the wild-type residue caused a clash score with neighboring residues, especially with the P104 residue when compared to wild-type (Fig. 2A and B).

The ConSurf web server determined that V433 residue is positioned in the variable site (conservation score of 2) in the protein, while it is located close to residues that are highly conserved in protein when compared to homologous proteins of CYP4F2 during evolution (Fig. 2C).

Functional impact prediction

Five tools have been utilized to identify the functional effect of the CYP4F2*3 (V433M) variant on the CYP4F2 protein (Table 3). PolyPhen-2, FATHMM-MKl, PANTHER-PSEP, CADD, and PhD-SNP^g tools have identified the CYP4F2*3 (V433M) as a deleterious variant.

MD simulation findings

200 ns MD simulations were conducted to analyze the deviation between native and CYP4F2*3 (V433M) variant proteins in physiological environments (Fig. 3). RMSD deviations produced during MD simulations were used to estimate the stability of the protein [55]. As presented in the RMSD plot (Fig. 3A), the native protein attains a relatively stable conformation after 20 ns. In comparison, the V433M variant was stable after approximately 70 ns of simulations suggesting that the V433M variant experienced more significant structural rearrangement before reaching a stable structure. In addition, the RMSD average values verified that the mutant structure had a slightly higher RMSD value (0.49 ± 0.05 nm) than the native structure (0.46 ± 0.04 nm), which means that the mutant form caused instability in the protein structure of CYP4F2.

Rg is an essential parameter for evaluating the dynamic adaptability of proteins and other biopolymers [56]. The Rg plot showed that the V433M mutant fluctuated more than that of the native form (Fig. 3B). In addition, the average Rg value analysis showed that the V433M mutation had a higher average Rg value (2.45 ± 0.02 nm) than the native protein (2.51 ± 0.01 nm) during the simulation time. This result indicates that the V433M mutation is less compact and more flexible than native CYP4F2.

SASA denotes the portion of the macromolecule’s surface accessible to the water solvent [57]. As shown in Fig. 3C, SASA analysis revealed that the V433M mutation had a higher average SASA value (250.74 ± 6.13 nm) than native CYP4F2 (245.17 ± 5.48 nm), suggesting that structural rearrangements in the V433M mutation resulted in increased protein expansion.

The total number of intramolecular hydrogen bonds was estimated during the 200 ns MD simulation as shown in Fig. 3D. From the analysis, it can be observed that native CYP4F2 forms a greater number of H-bonds, with an average of 371 ± 12.67, while the V433M variant exhibits fewer H-bonds, with an average of 350 ± 12.74. The lowest number of hydrogen bonds indicated an increase in the overall flexibility of CYP4F2 upon V433M variant.

RMSF based on residue displacement explains the thermal stability, local flexibility, and heterogeneity of macromolecules over MD simulation time [58]. The RMSF plot indicated that residues in the N terminal region had a major difference in fluctuation between the native CYP4F2 structure and the V433M variant after 200 ns MD simulation. In addition, the V433M variant exhibited a higher fluctuation region in the residues 244–249 of the Gˊ helix and residues 251–260 of the G helix, followed by increased flexibility in the residues 400–470 when compared to the native CYP4F2 structure (Fig. 3E). It was also noticed that the average of RMSF was increased from 0.08 ± 0.03 in native CYP4F2 to 0.11 ± 0.07 in V433M variant, suggested that the V433M mutation altered structural flexibility. These data are compatible with the structural alignment of native and mutant forms after molecular dynamics simulation, indicating that mutant form had a more open structure than the native form especially in region 400–470 which near to site of mutation V433M and showed significant structure changes (Fig. 4).

Changes in the secondary structure provide insight into the folding mechanism and conformational behavior of a protein. The DSSP plot (Fig. 5) showed that the native protein (Fig. 5A), and V433M variant (Fig. 5B) exhibited almost the same secondary structural arrangement. The data illustrated that native CYP4F2 contains an overall 60% secondary structure, including β-bridges, β-sheets, α-helices, and turns (Table 4). However, the V433M variant showed a slight decrease in the secondary structure profile compared that with of the native protein. Significant changes in the secondary structure content of the V433M variant can be seen due to bend and turn formation and loss of the A-Helix and 3-Helix without any change in β-sheets and coils.

PCA uses covariance matrices of Cα atoms to calculate the significant motions of atom pairs associated with vital biological functions. The first two principal components (PC1 and PC2) of the native and V433M variant proteins were generated by projecting the trajectories onto the respective eigenvectors. As seen in Fig. 5C, the V433M variant covered the slightly higher subspace spanned along the two eigenvectors (PC1: -3.56 nm to + 14.12 nm and PC2: -11.27 nm to + 8.76 nm) compared to the native (PC1: -2.31 nm to + 15.65 nm, and PC2: -5.2 to + 15.12 nm). It was also noticed that the trace of the diagonalized covariance matrix of the native form and V433M variant was calculated to be 64.13 and 69.99 nm, respectively. The rmsf eigenvector plot (Fig. 5D), based on atom number analysis, illustrated a higher V433M motion than the native form. In addition, the Gibbs free energy landscapes (FELs) of the first and second PCA showed that the V433M variant (Fig. 5E) had wider global energy minima (18.5 KJ/mol) than the native form (17.9 KJ/mol) (Fig. 5F), which indicates more conformational changes of the protein in the defined space. These results confirmed an increase in the overall flexibility of the V433M variant, which is in agreement with the RMSD, Rg, H-bonds, and RMSF values, as well as the SASA and DSSP analyses.

Molecular docking

The enzymatic activity of the CYP4F2 depends on its interaction with VK1. For this reason, molecular docking using the HADDOCK web server was performed to figure out the effect of the V433M variant on the interaction of the CYP4F2 protein with VK1 (Fig. 6). The binding site of VK1 on CYP4F2 protein was characterized according to a previous work of Li et al. (2018), that predicted binding site of VK1 in this protein. Based on this research, the amino acids including Trp59, Trp61, Met92 Phe124, His236, Phe327, Glu328, Val397, and Leu504 are responsible for interacting with VK1 [59]. The docking result (Table 5) showed that although the mutated form could still interact with VK1, it had a lower binding affinity (-8.5 ± 0.31 kcal/mol) than the native form (-10.0 ± 0.07 kcal/mol), which caused less activity of the mutated state compared to the native form. Bond-type analyses illustrated that in native form, VK1 formed strong forces such as hydrogen bonding, Pi-Pi stacking, and Pi-sulfur interactions with Val397, Met92 and Phe124 amino acids, whereas, in the mutant form, only Phe60 residue formed a Pi-Pi stacking with VK1. This caused VK1 to have more strong interaction with the native form than the mutant (Fig. 6).

Competing interests

The authors state that they have no interest conflict.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.