Comparative molecular dynamics simulation of apo and holo forms of the P53 mutant C176F: a structural perspective

Abstract

Zinc fingers represent a highly diverse structural domain, with P53 being a notable example among zinc-dependent transcription factors. The folding patterns of proteins in the cell are heavily influenced by the concentration of zinc. The potential for zinc loss arises due to its dysregulation and the frequent occurrence of tumorigenic P53 mutations. This could lead to significant consequences such as protein misfolding and a reduction in tumor-suppressing capabilities. To gain deeper insights into the structural conformation, stability, flexibility, and compactness of the zinc-binding mutation C176F, a comprehensive comparative computational analysis was conducted on the wildtype (WT) and mutant (MT) P53 in the presence (Holo) and absence (Apo) of zinc. This analysis was performed using molecular dynamic simulation. The overall observation was that the mutation in C176F reduces the metal binding affinity and results in less stability in the Apo and Holo P53 MT protein.

KEYWORDS:

• P53
• tumour suppressor
• zinc binding residues
• DNA-binding domain
• zinc finger
• molecular dynamic simulation

1. Introduction

P53 was first identified in 1989 and categorized as a cellular (simian virus 40) SV40 large T antigen-binding [1,2]. The tumour suppressor P53 gene, which is found on chromosome 17p13.1, is frequently referred to as “the guardian of the genome” or “the cellular gatekeeper of growth and division” [3,4] due to its function in maintaining genomic stability. Other names, such as the “Death Star” [5] and “good/bad cop” [6], further emphasize its dual role in promoting cell survival or death. Its dynamic nature in tumorigenesis has also led to it being described as an “acrobat” [7]. These names all underscore the multifaceted and pivotal role of the P53 gene in cellular function and cancer biology. The P53 gene consists of 11 exons, transcribing to 2.8 kb of mRNA and translating into a 53-kDa protein. According to The Cancer Genome Atlas (TCGA), mutations in the P53 gene occur in approximately 56% of all cancer types [8,9]. This makes P53 the most frequently mutated gene in human cancers, as it is altered in over half of all cancer cases [10]. Unlike many other tumour suppressor genes that often undergo deletion or truncation in cancer cells, most P53 allele mutations are missense mutations. These result in a single amino acid substitution in the full-length mutant (MT) protein [11].

Standard classifications for P53 mutations include DNA contact, structural/stability, and zinc binding. The N-terminal, the central DNA binding domain (DBD), and the C-terminal constitute its three domains. The DBD is a 24.7 kDa protein with a protruding 1, 2, and 3 loop core and a β-sandwich core with DNA interaction and consists of a Zn2 + ion [12]. Our study focuses on the zinc-binding class, where the Zn2 + ion has been found to have a role in P53 activity by composing the motions of several structural components of the P53 protein necessary for the DNA binding region [13]. P53 interacts with zinc, a d-block transition metal, stabilizes protein conformation, helps proteins interact with one another, and sets active sites for enzyme catalysis and electron transport reactions [14]. Studies through molecular approaches revealed that zinc triggers refolding of MT P53 protein by higher reactivity to the pAb1620 antibody and reconstructing WT P53 binding to promoter p21 in response to the drug [15]. Apo and holo protein structures of P53 was visualized using BIOVIA Discovery Studio Visualizer Figure 1(A and B). Zn2 + ion coordinates the residues Cys176, His179, Cys238, and Cys242 in a pseudo-tetrahedral coordination geometry with the distance of 2.5, 2.06, 2.4, and 2.5 Å respectively that are located in loop 2 and 3 [Figure 1(B)]. Stabilizing the loop L2 (Cys176-His179) and L3 (Cys238-Cys242) in the domain helps loop 3 to position into the DNA groove for consensus sequence accession. Studies found that the loss of zinc ions impairs loop 3, loop 2, and surrounding regions, disrupts DBD specificity, and affects the ability to distinguish consensus from non-consensus DNA [16]. Mehrian-Shai et al. reported that specific cancers are inactivated by the overexpression of metallothionein and zinc restriction of WT P53 [17]. Conformational flexibility with zinc and without zinc in p53 was studied, and it was demonstrated that zinc-free DBD (apo-DBD) differs structurally from DBD and is prone to aggregation [16,18,19].

Figure 1. The WT P53 protein in apo and holo forms were visualized using BIOVIA Discovery Studio Visualizer (A) Apo form of P53 protein, (B) Holo form of P53 protein. The distance between the Zn2 + ion and the binding residues are represented.

The zinc-chelating residue C176 is next to the most prevalent p53 mutation in cancer among other zinc-binding residues, significantly reducing metal binding affinity [20]. We studied the mutation position C176, which resulted in six prevalent mutations (C176F, C176Y, C176S, C176W, C176R, and C176G), and they are reported among various cancer types. According to reports retrieved from COSMIC (Catalogue Of Somatic Mutations In Cancer) data on the P53 gene, the mutation position of C176F was found to have the highest count of 301 [21]. It was observed that C176 mutations were detected in the tumour samples obtained from the patients, where the number of samples with detected mutation across various cancers is given in Figure 2(A), out of which C176F mutation samples were more in number than other mutations. TCGA data reported that C176F was found in 13 cancer types among various mutations recorded compared to other mutations [Figure 2(B)]. Hence, our study's primary focus is to validate the impact of mutation C176F in the P53 gene. P53 C176F lies within the DNA-binding domain of the p53 protein [22]. C176F was predicted to experience a loss of function to the p53 protein due to the inability to induce CDKN1A (p21) transcription in cell culture [23].

Figure 2. Various cancer types with C176 mutation samples (A) Data reported by a database of Catalogue Of Somatic Mutations In Cancer (COSMIC) and (B) TCGA cancer report.

On the other hand, molecular dynamics simulation (MDS) is a crucial tool for understanding the impact of C176 mutations on the structure of the P53 protein since it promptly offers information about the protein at the atomic level. Earlier molecular dynamics was used in various studies to examine the effects of P53 mutations [24–28]. We performed a molecular dynamic simulation of apo WT and MT C176F (absence of zinc) and compared those results with the holo WT and MT C176F (presence of zinc) simulations to investigate the mechanisms of destabilization and effects of mutation on P53 structure and dynamics. Initially, we screened C176F for pathogenicity analysis using various in silico tools to identify its deleterious effect on the structure. We concentrate on the apo and holo simulations of WT and MT structure C176F for our study. Overall, the conformational flexibility of P53 in apo and holo WT and C176F MT structures will be revealed through extensive analysis of molecular dynamic simulation performed for 500 ns.

2. Methods

2.1. Stability prediction

Several stability prediction techniques were performed to analyze the impact of the amino acid changes at the 176th position on the stability of the WT P53 protein. The ConSurf web server analyzed the evolutionary conservation [29]. The stability of the protein was predicted using the webservers, MuPro (accessed on 15 May 2023) [30], I-Mutant 2.0 (accessed on 15 May 2023) [31], PremPS (accessed on 16 May 2023) [32], Dynamut 2 (accessed on 16 May 2023) [33], and INPS-MD (Impact of non-synonymous variations on protein stability multi-dimension) (accessed on 16 May 2023) [34]. The evaluation and prediction of the protein stability were estimated, and the sequence was mutated using the CUPSAT programme (accessed on 16 May 2023) [35]. This technique uses potentials of torsion angle and structural environment-specific atoms to calculate the difference between WT and MT proteins in the free energy of unfolding. The i-Stable (accessed on 16 May 2023) used three algorithms to compute the stability of mutated P53 protein using sequential input [36]. The HOPE server (accessed on 16 May 2023) was used to generate predictions regarding the effects of the mutation on the 3D structure and related functions.

2.2. Molecule preparation

The X-ray crystal structure for the P53 protein was obtained from RCSB-PDB, with an ID: 2OCJ, a homotetramer with 2.05 Å (96–289 in length) resolution in the absence of DNA, which resulted in new tetrameric arrangements [37]. The P53 gene complex was mutated using SWISS-PDBviewer [38]. The prepared models apo WT (absence of Zn2 + ion), apo MT C176F (absence of Zn2 + ion), holo WT (presence of Zn2 + ion), and holo MT C176F (presence of Zn2 + ion) were considered for further analysis.

2.3. Molecular dynamic simulation

A comparative MDS of 500 ns was carried out for apo and holo (WT) and apo and holo (MT) C176F of P53 protein using CHARMM-GUI [39] to create four systems as apo WT, apo MT, holo WT, and holo MT. Each system was developed using the input generator's option solution builder. A rectangular TIP3 water box with a 10 Å distance between its edges was used to solvate each system. The CHARMM36 all-atom force field analyzed each of the four system's topologies, and coordinated files were generated. The MDS was carried out at 300 K physiological conditions and neutral pH. They were all neutralized by adding counter ions to each system in the simulation box. The GROMACS v2019.4 simulation software was utilized for the simulation process [40–44]. The command pdbgmx2 was used to add hydrogens to the heavy atoms and convert the protein coordinates. The SPC water molecule was replaced with Na⁺ and Cl^- counter ions to maintain salt concentration using the plugin genion from the GROMACS programme to neutralize water atoms. Emtol convergence criteria were used to confine the bond length by increasing the force, and the maximum force applied was 1000 KJ/mol, and the steepest technique was used to keep the temperature constant with 300 K for 50,000 steps. To remove irregular torsions, the neutralization was restricted to energy minimization. The electrostatic interaction has been calculated using the Particle Mesh Ewald (PME) method. Using the number of atoms, volume, and temperature of the system (NVT), the pressure and temperature of the system (NPT), and to further equilibrate the system, progressive removal of the constraint on heavy atoms was applied with those two steps. Conformational changes were made to tune the processing system to 500 ns for MDS.

For trajectory analysis, the GROMACS utilities “gmx rms”, “gmx rmsf”, “gmx gyrate”, and “gmx sasa” which stands for root-mean-square deviation (backbone), root mean square fluctuation, the radius of gyration, and solvent accessible surface area respectively were used. Dynamic cross-correlation matrices were used to determine the secondary structure of the protein. The kernel Density Estimation (KDE) was done by comparing Rg and SASA for Protein compactness and surface area. The Gibbs Free Energy landscape can be used to understand the protein’s stability, folding, and function [45]. Proteins-associated movements are extracted by covariance analysis, often called principal component analysis (PCA) or essential dynamics (ED), to understand the motions that are crucial to a protein's function [46]. PCA is an extensively utilized method for explaining the observed motion changes in the protein across the simulation system. A linear mapping’s eigenvalue measures the distortion caused by the transformation, while eigenvectors reveal the distortion pattern. The patterns of secondary structure were observed by DSSP analysis [47]. The resulting data were evaluated using the software called PyMol v2.5 [48], VMD v1.9.3 [49], and Chimera v1.16 [50] to visualize the figures. The XmGrace tool was used to plot the graph.

3. Results

3.1. Stability prediction of C176F using in silico tools

Consurf was used to determine the degree of conservation for a given amino acid. Cysteine at position 176 is shown as 9, marked as a highly conserved region [Figure S1(A)]. For stability detection of C176F, i-Mutant 2.0, MUpro, Prem-PS, INPS-MD, and iStable scored −1.61, −1.04, 0.98, −2.76, and 0.60 Kcal/mol with decreased stability. Dynamut 2 and CUPSAT with a score of −0.99 and −1.5 Kcal/mol exhibited destabilizing effects (Table 1). C176F scored 0.99 by the HOPE server, which makes the prediction score based on dbNSFP, and the MetaRNN score ranges from 0.0 to 1.0. The score is directly proportional to the pathogenicity, which implies that with the increase in the score, the pathogenicity increases.

Table 1. Prediction of stability change by DDG for the variant C176F.

Tools	DDG kcal/mol	Stability
i-Mutant2.0 SEQ	−1.61	Decrease
MUpro	−1.04	Decrease
PremPS	0.98	Decrease
INPS-MD	−2.76	Decrease
CUPSAT	−1.5	Destabilizing
DynaMut 2	−0.99	Destabilizing
iStable	0.604838	Decrease

3.2. Molecular dynamic simulation

Mutations in the P53 gene can cause conformational changes, resulting in no or poor protein synthesis [22,51]. To comprehend the structural and functional characteristics of the reported deleterious mutation, we used MDS to compare the WT and MT (C176F) proteins in or absence of zinc. GROMACS was executed on four systems for 500 ns MDS. Trajectory files of the apo WT and MT, and holo WT and MT structures were prepared and subjected to comparative studies such as the Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), Radius of Gyration (Rg), Solvent-Accessible Surface Area (SASA). Defining the secondary structure of the protein (DSSP) was performed to analyze the secondary structure pattern. KDE analysis, DCCM, FEL, and PCA after the MDS were used for free energy calculation.

3.3. Trajectory analysis

3.3.1. Structural stability

To determine the stability of the P53 protein, we examined the RMSD of the backbone for WT and MT in apo and holo forms. The acquired modifications by setting 500 ns as a function of time notable changes in the patterns of RMSD average values were analyzed (Table 2). We performed the triplicates for the RMSD of the backbone and Cα atoms for all the structures as Run 1, 2, and 3 [Figure S2(A–D)]. During apo simulations, the scores obtained by WT and MT were up to ∼0.35 nm and ∼0.30 nm at 500 ns, respectively. Higher deviation patterns were observed throughout the simulation process in apo MT between 250 and 350 ns compared to apo WT, which showed reduced deviations (Figure 3(A)). Holo simulations expressed WT at ∼0.3 nm at 500 ns (Figure 3(B)), and MT C176F showed deviations at ∼0.45 nm at 250–300 ns, while more deviations were found between 100–400 ns and ∼0.27 nm at 500 ns as shown in Figure 3(B) as unstable due its higher deviation compared to WT throughout the simulation. Therefore, the MT forms in apo and holo indicate the cause of instability during 500 ns. Given that proteins require a suitable and stable structure to function, from the acquired data, we conclude that Apo and Holo MT might alter the structure of the P53 protein.

Figure 3. The Backbone Root Mean Square Deviation (RMSD) and value for WT and MT C176F of P53 protein in Apo and Holo simulations at 500 ns (A) Apo WT (Grey), Apo MT (red), (B) Holo WT (Grey), Holo MT (Blue). The Root Mean Square Fluctuation (RMSF) of Cα atoms for WT and MTC176F of P53 protein in Apo and Holo simulations at 500 ns (C) Apo WT (Grey), Apo MT (red), (D) Holo WT (Grey), Holo MT (Blue).

Table 2. The average value of time of structural properties for RMSD, Cα-RMSF, Rg, and SASA are calculated for WT and MT of Apo and Holo-structure of P53 protein.

	Apo		Holo
	WT	MT	WT	MT
Backbone RMSD (nm)	0.25	0.28	0.26	0.28
Cα-RMSF (nm)	0.43755	0.6329	0.4926	0.50905
Rg (nm)	1.665488957	1.673102703	1.66618413	1.680861024
SASA (nm^2)	108.003	108.859	108.889	109.114
Trace of the covariance matrix diagonalized (nm^2)	5.60441	6.59661	5.35048	7.08271

3.3.2. Structural flexibility

To clarify how CYS 176 and the mutation to PHE 176 affect the dynamic behaviour in apo and holo forms, Cα-RMSF was calculated. In apo WT showed fluctuations of 0.45 nm (Figure 3(C)), while apo MT C176F showed a slightly higher fluctuation of 0.46 nm (Figure 3(C)). The holo simulation structure’s RMSF calculation of WT showed fluctuations at 0.3 nm (Figure 3(D) and Table 2), while holo MT C176F fluctuations were at 0.55 nm with higher fluctuation than holo WT (Figure 3(D)). The apo MT C176F showed higher average Cα-RMSF values than other structures (Table 2). These results reveal that apo MT C176F is more flexible. A fluctuation value of less than 2Å is acceptable for a small protein. The L1 and L2 loop region of the DBD domain shown in Figure 3(C and D) had the majority of fluctuation in all cases of WT and MT of both the apo and holo forms, and residues from 160 to 251 experienced considerable fluctuation, between where our target binding site 176 of Zn2 + ion is located. The results imply that the mutation C176F influenced the protein’s overall flexibility rather than the residual level. In our case, the C176F affects the flexibility of the P53 protein in both the apo and holo form. The RMSF results and the RMSD readings are in accord.

3.3.3. Structure compactness

The mass-weight root mean square of the distance between the group of atoms is calculated using the Rg of particular atoms from the common center of mass. Rg provides details about the protein's structure. Rg is a crucial factor in determining proteins’ dynamic stability and compactness. The association between a protein's flexibility and compactness is widely established. The amount of the protein surface exposed to solvents is also influenced by how compact the protein appears. Rg for apo WT and MT was 1.66 and 1.67 nm, respectively; apo MT was slightly less compact than apo WT (Figure 4(A)). The Rg of protein was computed for holo WT and MT to understand the influence of Zn2 + ion on protein compactness, and the scores obtained were 1.67 and 1.71 nm, respectively. The results implied that the holo MT possessed less compactness comparatively (Figure 4(B)). The fluctuation score of Rg is given in Table 2. Rg analysis suggests that the holo MT is more flexible in conformation, agreeing with RMSD and RMSF values.

Figure 4. Radius of gyration (Rg), for Apo and Holo P53 protein simulations for WT and MT at 500 ns (A) Apo simulations for WT and MT C176F, (B) Holo simulation for WT and MT C176F, and Solvent Accessible Surface Area (SASA) during 500 ns representing Apo and Holo simulation versus time, (C) SASA of Apo WT and MT C176F, and (D) SASA of Holo WT and MT C176F.

3.3.4. Analysis of surface area

To interpret how zinc ions affect the surface of the P53 protein, we computed solvent-accessible surface area (SASA). SASA evaluates how much an amino acid involves interaction geometrically (the solvent and the core of the protein). The number of liberated zinc-binding residues affects the tertiary structure and its SASA level. This relationship between amino acid and how it is released under various environmental conditions in nature is inverse. Higher SASA values are seen in dispersed protein structures, while lower SASA values are found in compact protein structures. A rise or reduction in the SASA values indicates a structural change in the protein. According to the examination of SASA simulations, apo WT varied at 108 nm² and MT at 109 nm², and holo WT and MT varied with time at 110 and 117 nm², respectively, as shown in Figure 4(C and D). WT displayed a lower average total SASA of 108.003 and 108.889 nm² (Figure 4(C)). In contrast, MT structure displayed 108.859 and 109.114 nm² (Figure 4(D)) of apo and holo-structure, respectively. MT structures showed higher SASA values compared to WT (Table 2). Therefore, SASA values also suggest following RMSD, RMSF, and Rg values.

3.4. Free energy calculation

The Kdeplot, also known as a Kernel Distribution Estimation Plot, is a graphical representation of the probability density function of non-parametric or distribution-free data variables. The data are plotted against multiple or bivariate variables using the Seaborn Kdeplot function. The area under the depicted curve serves as a representation of the probability distribution of the data values. We created a Kdeplot with additional capabilities using the Python Seaborn module. The colour palette was given using “cmap” parameters. Protein flexibility and compactness were associated. The compactness of apo and holo P53 protein also influences the surface area exposed to the solvent. SASA was calculated to examine the effects of MT metal-binding residues on protein surface area, and Rg was computed to examine the influence of metal-binding residues on protein compactness. The KDE analysis of SASA and Rg were compared and plotted, as shown in Figure 5. The Rg and the SASA value was at 1.66 nm and 107 Ų in apo WT (Figure 5(A)), with the highest 1.67 nm Rg value, and 111 Ų SASA value in apo MT (Figure 5(B)). In Holo simulation, the Rg and SASA value of holo WT (Figure 5(C)) was 1.67 nm and 110 Å^2, respectively. Holo MT's Rg value was 1.68 nm and 113 Ų SASA (Figure 5(D)). Therefore, it was observed that holo MT loses compactness with a higher surface area. In apo MT form, the compactness is the same but increases with surface area.

Figure 5. The Kernel Density Estimation (KDE) correlation plot of Rg and SASA. (A) Apo WT, (B) Apo MT, (C) Holo WT, and (D) Holo MT. The probability distribution is shown by concentric circles that are differentiated by the intensity.

3.5. Dynamic cross-correlation analysis

Comparative dynamic cross-correlation matrix (DCCM) analysis was carried out for apo and holo of WT and apo and holo of MT to calculate the contrast in the structure and dynamics (Figure 6). DCCM plots evaluated correlated movements between residues from various regions, observed throughout 500 ns trajectories of each structure complex. The correlation range in the DCCM plot was from +1 and −1, representing variation in the colour intensity; reddish orange-colored contours correspond to positively correlated movements, while yellow-colored regions represent zero correlation or uncorrelated movements, and dark blue-colored contours indicate anti-correlated movements. The characteristic features in DCCM plots are shown in squared regions labelled 1, 2, and 3 to illustrate C176F mutation on the structure and dynamics of the apo and holo complex. In apo WT regions, 1 and 2 correspond to positively correlated motions (Figure 6(A)), whereas in apo MT, region 1 has zero to negatively correlated motions, and region 2 with positive to zero correlation motions (Figure 6(B)). In holo WT, regions 1 and 2 had mixed correlations (Figure 6(C)); holo MT contained positive to zero correlations (Figure 6(D)). Region 4 and 5 represent loops 1 and 2, respectively, positively correlated.

Figure 6. Dynamic Cross-correlation matrix (DCCM) of backbone atoms. (A) Apo WT, (B) Apo MT, (C) Holo WT, and (D) Holo MT. DCCM maps were obtained from 500 ns trajectories that contain contours of different colours. Correlated motion and anti-correlated motions are represented by reddish-orange and dark blue contours, respectively. The square shape represents residue motion and the circle and ellipse shapes resemble loop 1 and loop 2 regions respectively.

3.6. Gibbs free energy landscape

The projection of their first (PC1) and second (PC2) eigenvector of each WT and MT structure of the apo and holo complex was used to examine the Gibbs energy landscape using gmx covar, gmx analog, and gmx sham (Figure 7). These free energies indicate the stable conformation state of the molecules. The folding pattern of a protein is crucial for its proper function. Red-colored regions represent the lowest energy or stable conformation states, whereas blue-colored regions correspond to higher energy with unstable conformation states.

Figure 7. Gibbs free energy landscape of PC1 & PC2 (initial two principal components) of projection of simulated trajectories. (A) Apo WT, (B) Apo MT, (C) Holo WT, and (D) Holo MT obtained during 500 ns MD simulation.

The Apo WT shows minimum red contour with more blue regions leading to higher energy and unstable states within the energy landscape (Figure 7(A)) than the Apo MT separated by more red areas with lower energy and stable conformation (Figure 7(B)). The Holo WT structure represents a minimum red region (Figure 7(C)) compared to Holo MT which showed more red areas representing lower energy and stable conformation (Figure 7(D)). FEL examined the fluctuation in the two structures for all Cα atoms of apo WT and MT was 0–16.9 and 0–17, while Holo WT and MT ranged from 0 to 16.7 and 0 to 16.6 respectively. The Apo WT complex was found to have higher energy with unstable conformation from the four complexes.

3.7. Principal component analysis

PCA was then utilized to forecast the variation in dynamic behaviour between WT and MT P53 protein molecules. The principal component of a covariance matrix was created using its eigenvectors. Changes in the trajectory files generate a projection (eigenvector changes) by projecting acquired trajectories onto the two major components (PC1 and PC2). A matrix was built using the eigenvectors of the atoms in the molecule. This matrix illustrated how the protein's atoms shifted. Eigenvectors define the crucial subspace with the highest associated eigenvalues in which most protein dynamics occur. The apo WT covered the space in PC 1 and the apo MT with less space (Figure 8(A)), whereas the holo WT structure occupied more space, while the holo MT covered more significant regions (Figure 8(B)). The trace of the diagonalized covariance matrix of Cα atomic variations in position was used to quantify the protein's flexibility. The trace of covariance matrix value for apo WT was 5.60441 nm², apo MT was 6.59661 nm², and holo WT was 5.35048 nm², and holo MT was 7.08271 nm² consisting among all, the holo MT protein seemed to possess higher values from observation (Table 2). Hence, from the data obtained, the MT protein structure showed a higher flexibility at 300 K than the WT structure. In summary, the study's validity is enhanced by the concordance between the PCA and RMSD, RMSF, Rg, and SASA results. In addition, we used an ED analysis to support our MDS findings and pinpoint that the WT and MT motions were associated with one another along the trajectory through simulations. Since the total of the eigenvalues represents the system's overall motility, we displayed the eigenvalues for the first 40 modes of motion for WT and MTs during simulations against the corresponding eigenvector index in Figure 8(C). Most of the protein's internal motion is confined to a narrow dimension in the critical subspace, as shown by the fact that only a few eigenvectors in the apo and holo simulations had significant eigenvalues for both WT and MT. Figure 8(C) of the eigenvalue spectrum shows that the system's primary fluctuations were concentrated within the first forty eigenvectors.

Figure 8. Principal Component Analysis of the first two eigenvectors (A) Apo WT (black) and MT structure (red), (B) Holo WT (mint green) and MT (orange) structures and (C) eigenvalues for the first 40 modes of motion.

3.8. Secondary structure analysis

We utilized the DSSP (Define Secondary Structure of Protein) tool to evaluate the differences in secondary structure variations between WT and MT. This tool assigns secondary structure labels to the protein residue using H-bonding patterns and several geometrical properties. The do_dssp function, a GROMACS built-in feature, was used to examine how the secondary structure developed over time. In addition, the surface area of the secondary structures at the start of the simulation and at specific time steps where conformational changes occurred at a higher range was analyzed between the WT and MT forms. The patterns of α-helices, β-sheets, coils, bends, and turns are evaluated in the average percentage of structure 54%, 53%, 54%, and 50% in the apo WT, apo MT, holo WT, and holo MT, respectively (Table 3). The Four secondary structure features of apo and holo-structures of C176 show slight changes between MT and WT structures at different nanoseconds, namely 0 ns, 250 ns, and 500 ns (Figure S3). It was observed that during the simulation, the apo MT structure was coiled [Figure S3(B)], the Holo WT structure into a beta-sheet [Figure S3(C)], and the holo MT forming beta bridge (Figure S3(D)). An investigation of time-dependent secondary structure variations with distance was calculated using gmx distance by GROMACS, calculated for the centre of mass between the Loop 1 (PHE113:N-THR123:O), Loop 2 (LYS164:N-CYS176:O, CYS182:N-LEU194:O) and Loop 3 (MET237:N-PRO250:O) distances between apo WT and MT protein with the C176 binding site found to be 2.8 nm and 2.8–3 nm, respectively [Figure S4(A)]. In Holo WT and MT structures, distances were calculated for the C176 binding site and Zn2 + ion that were 2.7 and 2.4 nm, respectively [Figure S4(B)].

Table 3. The average percentage of secondary structure values from trajectory analysis of 500 ns for WT and MT of Apo and Holo P53 protein.

	Type	Structure (%)	α-Helix (%)	β-Sheet (%)	Coils (%)	Bends (%)	Turns (%)
Apo	WT	54	6	32	31	15	13
	MT	53	5	34	29	15	14
Holo	WT	54	5	33	29	15	15
	MT	50	5	33	33	15	11

4. Discussion

The P53 protein, which carries a heterozygous C176F mutation, retains its function as a transcription factor. The ConSurf tool predicts that the Cysteine at position 176 is highly conserved. The HOPE servers suggest that the C176F mutation could damage the P53 protein. When Cysteine in P53 is mutated to Phenylalanine, it is not properly positioned due to its larger size and its location in a region with known splice variants. This disrupts the interaction with the Zn2 + ion [Figure S1(B)]. The C176F mutation has been linked to sporadic cancer. The substitution involves replacing a hydrophobic amino acid at position 176 with a polar side chain cysteine, which is then replaced with the non-polar Phenylalanine. This can disrupt bonding or cause a structural change in the protein [52]. The wildtype (WT) C176 forms a hydrogen bond with H179. However, due to the size difference of the mutant type (MT) residue C176F, there is no hydrogen bond as in the WT residue [Figure S1(C)]. Therefore, the mutation of the C176 residue to Phenylalanine in P53 eliminates the coordination with zinc, which is crucial for maintaining the structural integrity of P53. This could significantly impact the structural stability of the P53 protein [53]. To understand the effects of this mutation on both the apo and holo forms, we examined the amino acids within a 10 Å radius at the 0 and 500 ns time points in the trajectory files (Table S1). In the apo form, both the WT and the MT shared 20 amino acid residues at 0 ns and 22 residues at 500 ns. For the holo form, both the WT and MT had 21 common amino acid residues at both 0 and 500 ns time points.

We performed MDS to draw in silico validations on the apo and holo form structures to investigate destabilization mechanisms and the effects of mutation on the P53 structure. Our reasoning through various computational investigations was supported by the fact that under usual physiological conditions, a significant fraction of P53 may exist in the zinc-free state or apo P53. P53 and apo P53 were thermodynamically stable and kinetically accessible [16]. To understand how the variant C176F changes the balance of P53 conformation, we excluded the interference of external organic molecules, such as DNA and other proteins, to avoid a negative impact on simulation speed and increased system complexity. Zinc was believed to support protein structural integrity without interacting with DNA [12]. Our study of the RMSD and RMSF analysis found that apo and holo MT had altered structural instability; thus, MT P53 might significantly influence susceptibility towards cancer. The KDE analysis for RG and SASA was done to find the association of protein flexibility, compactness, and surface area exposed to solvent, which was also influenced by the compactness of the MT P53 protein. The compactness in the apo MT structure was the same but with increased surface area, while the holo MT structure lost its compactness and surface area. Covariance plots of the apo complex and DBD exhibited similar patterns for secondary structure analysis [13]. Our results from the covariance matrix for the apo MT structure showed two regions with zero to negative and positive to zero correlation, respectively. At the same time, holo MT during DCCM had positive to zero correlations, FEL analysis showed lower energy and stability as holo WT, and PCA covered a larger space. A study on secondary elements by Duan and Nilsson on P53 DBD predicted that L2 was highly flexible in apo DBD [13]. We conducted a study on the flexibility of the loops in both the apo and holo forms of the WT and MT proteins, specifically at the 0 and 500 ns time points, to identify any structural alterations in the loop region (Figure 9). The angle calculations were executed using the 2021 version of Maestro from Schrodinger, which allowed us to observe any structural modifications in the P53 protein. The central amino acid residues ALA 119, GLY 244, and GLY 187 of loops L1, L2, and L3 respectively, were kept constant during this process.

Figure 9. The flexibility of the loops was analyzed by angle calculation using the software Maestro v 2021. (A) Apo WT, (B) Apo MT, (C) Holo WT, and (D) Holo MT Colour representation as follows: Loop 1 (dark orange in colour), Loop 2 (dark green in colour), and Loop 3 (dark pink in colour) during 0 ns (Grey in colour), and Loop 1 (light orange in colour), Loop 2 (light green in colour), and Loop 3 (light pink in colour) during 500 ns (Blue in colour) of MD simulation.

In the apo WT state at 0 ns, the distances between ALA 119 (L1), GLY 244 (L2) and GLY 187 (L3) are 28.74, 23.53, and 33.16 Å respectively. These distances change to 35.95, 29.97, and 29.22 Å respectively at 500 ns. The angle between L1, L2, and L3 decreases from 78.0° at 0 ns to 51.7° at 500 ns (Figure 9(A) and Table 4). In the apo MT state, the corresponding distances at 0 ns are 30.59, 22.38, and 32.44 Å, and they change to 28.17, 27.24, and 25.27 Å at 500 ns. The angle decreases from 78.9° at 0 ns to 51.4° at 500 ns (Figure 9(B)). In the holo WT state, the distances at 0 ns are 28.27, 22.61, and 32.64 Å, and they change to 31.32, 20.31, and 35.59 Å at 500 ns. The angle increases from 79.0° at 0 ns to 84.3° at 500 ns (Figure 9(C)). In the holo MT state, the distances at 0 ns are 28.62, 23.99, and 33.03 Å, and they change to 27.62, 15.86, and 30.88 Å at 500 ns. The angle increases from 77.2° at 0 ns to 83.6° at 500 ns (Figure 9(D)). It is observed that the apo MT state at both 0 and 500 ns and the holo MT state at 500 ns significantly impact the stability of the L2 region, with a smaller angle. The angles for the apo structures (both WT and MT) are similar at 0 and 500 ns. In the Holo structure, the MT has a smaller angle than the WT form at both 0 and 500 ns. Thus, comparing the WT and MT of both apo and holo forms, the apo forms have smaller angles at 500 ns, and the holo MT has a smaller angle at 0 ns than at 500 ns. The change in angle of degree at 0 and 500 ns of apo WT, apo MT, holo WT, and holo MT are 26.3°, 27.5°, 5.3°, and 6.4° respectively. Thus, the values range from apo MT > apo WT > holo MT > holo WT with apo MT is higher and holo WT being less flexible. In summary, these observations underscore the significance of the presence or absence of a ZN and the impact of mutations on the stability and structural dynamics of the L2 region in the P53 protein. The bond angle of three loops likely refers to changes in the protein’s conformation and structural flexibility, and these observations suggest that both the ZN binding and the mutation have significant effects on these structural aspects, potentially influencing the protein’s function and behaviour.

Table 4. The calculation of angles and distances for L1, L2, and L3 was conducted using the 2021 version of Maestro from Schrodinger.

Protein Type	Apo WT		Apo MT		Holo WT		Holo MT
Simulation Period	0 ns	500 ns	0 ns	500 ns	0 ns	500 ns	0 ns	500 ns
Angle	78.0°	51.7°	78.9°	51.4°	79.0°	84.3°	77.2°	83.6°
L1 (ALA 119) to L2 (GLY 244)	28.74 Å	35.95 Å	30.59 Å	28.17 Å	28.27 Å	31.32 Å	28.62 Å	27.62 Å
L2 (GLY 244) to L3 (GLY 187)	23.53 Å	29.97 Å	22.38 Å	27.24 Å	22.61 Å	20.31 Å	23.99 Å	15.86 Å
L3 (GLY 187) to L1 (ALA 119)	33.16 Å	29.22 Å	32.44 Å	25.27 Å	32.64 Å	35.59 Å	33.03 Å	30.88 Å

Overall, the RMSD average values of apo WT have been observed to be less deviated than apo MT, holo WT and holo MT, but it is found that all the structures were at the range of 2 Å leading to stable MD trajectories. Small deviation reflects the stable nature of the protein, as zinc-free P53 is studied to be folded and stable by Butler et al in both apo and holo form [16]. In RMSF the values indicate fluctuation range from apo MT > holo MT > holo WT > apo WT as ∼0.63 nm, ∼0.50 nm, ∼0.49 nm and ∼0.43 nm respectively. From which it is not below 0.2 nm indicating the presence of loops leading to higher fluctuations (C238, C242-L2 and C176, H179-L3). Since loss of zinc disturbs the structure of L3 causing it to lose its ability to differentiate consensus from non-consensus DNA. It is also observed that apo MT and holo MT forms have higher values than apo WT and holo WT due to the size difference of C176 residue to mutant Phenylalanine with no hydrogen bond formation eliminating coordination to zn2 + ion. In vivo, Studies also found the P53 C176F mutation to be hypomorphic and destabilizing [54]. In compactness, a higher value reflects less compact, thus apo and holo MT have higher deviations of ∼1.67 nm and ∼1.68 nm respectively than apo and holo WT forms with ∼1.66 nm that leads to disturbance in hydrogen bond formation when mutated, determining all four structures are closely packed during MDS. From the angle distance calculations, it is found that holo WT and MT structures are less flexible than apo WT and MT, which is due to the coordination of Zn2 + ions in holo forms that are  intact than apo forms which are zinc free have higher fluctuation and flexibility.

According to computational studies, a P53 tetramer made up of two WT and two C176F MT molecules carried structural stability and DNA-binding interactions by forming more hydrophobic and cation-Π interactions to make up for the absence of zinc coordination [55]. Insights behind P53 dysfunction and drug development have been acquired from earlier attempts to describe DBD folding, although these studies are ultimately lacking because of insufficient data on the zinc-binding affinities in folded and unfolded states [56,57]. Several studies have suggested a core-core interaction between DBDs when bound to DNA and that the interface region was proposed to be around H1 in the L2 loop [58,59]. Additionally, a loss of Zn2 + can lead to misfolding, especially when combined with tumour-causing mutations [60,61]. It has been discovered that the patient-specific MT P53 C176F still exhibits some relevant WT p53 apoptotic activity [54].

MDS has been used in earlier research to explore the effect of P53 mutations by examining the zinc-binding residues, which helps in delivering information on the protein at the atomic level by evaluating the effect of mutation on protein structure [24–28]. In another study, MD simulations were done for P53 monomer protein on comparison with N-terminal, C-terminal, and DBD or core domain, which expressed more deviation and fluctuation in N and C terminal region compared to DBD region through RMSD and RMSF, respectively. However, at the same time, H1 through Zn2 + ions coordinated cysteine residues had conformational changes in the DBD region, including the loop-sheet-helix region, hypothesizing that helix plays a pertinent role in coordination to fulfil P53 biological functions [62].

A group did a study on P53-targeted therapies based on zinc ion delivery. They defined zinc metallochaperones (ZMC1) to reactivate MT P53 C176S by transporting Zn2 + ions from external sources into cells across the plasma membrane. It was found that it decreased zinc binding affinity, but the stability was stable [63]. Whole-exome sequencing of colorectal cancer tumour tissues from Taiwan patients revealed 50% of P53 mutations were identified in 16 samples, also, 15 variants were detected including C176F [64]. Deep sequencing revealed a P53 mutation frequency of 60%, similar to that of the cosmic database, and 48% detected by Odon and co-studied P53 sensitivity to cetuximab treatment [65]. A study using patient-derived xenografts that harboured P53 C176F mutation that was expressed in Embryonal rhabdomyosarcoma. In other cancers, such as esophageal adenocarcinoma P53, C176F was reported as a frequent mutation in Chinese patients by targeted sequencing via NGS platforms [66]. This method has worked for some mutations in mice models and cell cultures of cancer [20,67–69]. A population study was done on a group of patient samples from TCGA and other public databases with a high-grade ovarian severe cancer of P53 mutations for survival analysis, which identified better overall survival analysis in C176 mutations [70].

P53 folding and function have specific control of zinc availability and binding. Zinc deficiency is doubly deleterious [71], and various tumorigenic mutations destabilize DBD. A lack of Zinc can also lead to misfolding, especially when combined with tumour-causing mutations. Therefore, insights into structural conformation through computational approaches for apo and holo WT and MT zinc binding residue C176F through MDS and various trajectory analyses exposed that MT apo and holo structure causes effect by conformational changes into WT P53 protein leading to more suitable interaction study, thereby contributing to the rational development and design of new anticancer drugs.

5. Conclusion

The present study gives a detailed insight into the genotype-phenotype interaction of zinc-binding residues Cys176F of the P53, causing various cancers. As a result, stability prediction analysis confirmed that C176F of the P53 is destabilizing and crucial to cancer development. In both the presence and absence of Zn2 + ion, the apo and holo MT structure exhibited higher SASA values, and greater space occupied in PCA plots demonstrated their deleterious impact on the DBD of P53. The binding affinity was lower for C176F residue that disrupts structural stability. Overall, the current computational approach will provide a complete understanding of the structural insights and conformational changes that result in the P53 destabilizing processes in cancer. This knowledge may be helpful for the development of P53 MT-based therapeutic interventions against cancer types.

Author contributions

HCR and GPDC were involved in the study design and data collection. HCR and GPDC were involved in acquiring, analysing, interpreting the results, and drafting the manuscript. GPDC supervised the entire study. All authors edited and approved the submitted version of the article.

Acknowledgment

The authors would like to thank the Vellore Institute of Technology, Vellore, India, for providing the necessary research facilities and encouragement to carry out this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

No funding agency was involved in the present study.