Molecular dynamics simulations reveal the effect of mutations in the RING domains of BRCA1-BARD1 complex and its relevance to the prognosis of breast cancer

Abstract The N-terminal RING-RING domain of BRCA1-BARD1 is an E3 ubiquitin ligase complex that plays a critical role in tumor suppression through DNA double stranded repair mechanism. Mutations in the BRCA1-BARD1 heterodimer RING domains were found to have an association with breast and ovarian cancer by a way of hampering the E3 ubiquitin ligase activity. Herein, the molecular mechanism of interaction, conformational change due to the specific mutations on the BRCA1-BARD1 complex at atomic level has been examined by employing molecular modeling techniques. Sixteen mutations have been selected for the study. Molecular dynamics simulation results reveal that the mutant complexes have more local perturbation with a high residual fluctuation in the zinc binding sites and central helix. A few of the BRCA1 (V11A, I21V, I42V, R71G, I31M and L51W) mutants have been experimentally identified that do not impair E3 ligase activity, display an enhanced number of H-bonds and non-bonded contacts at the interacting interface as revealed by MD simulation. The mutation of BRCA1 (C61G, C64Y, C39Y and C24R) and BARD1 (C53W, C71Y and C83R) zinc binding residues displayed a smaller number of significant H-bonds, other interactions and also loss of some of the hotspot residues. Additionally, most of the mutant complexes display relatively lower electrostatic energy, H-bonding and total stabilizing energy as compared to wild-type. The current study attempts to unravel the role of BRCA1-BARD1 mutations and delineates the structural and conformational dynamics in the progression of breast cancer. Communicated by Ramaswamy H. Sarma

The BRCA1-BARD1 heterodimer complex has numerous known cellular targets, however the molecular mechanism of recognition of DNA DSB and ubiquitination specificity for substrates remain unknown.This lack of mechanistic and atomic level understanding has limited the knowledge of the molecular functions of this heterodimer and how mutations in both BRCA1 and BARD1 relate to the breast as well as ovarian cancer susceptibility.Thus, studying the pathogenicity of cancer-predisposing mutations identified in different domains of BARD1 and BRCA1 is of outstanding importance (Choudhary et al., 2018, Choudhary et al., 2016;Kumar et al., 2022, Barua et al., 2022).Several researchers have attempted to comprehend the atomic level molecular mechanism of PPI through MD simulation studies in order to pinpoint the hotspot residues for drug designing and the influence of mutation on protein's structure and function (Sarma et al., 2022;Sharma & Sastry, 2015;Badrinarayan & Sastry, 2011, 2014;Srivastava et al., 2011;Kumar et al., 2021;Sarma & Sastry, 2022;Sarvagalla et al., 2016Sarvagalla et al., , 2021;;Jha et al., 2019;Priyadarsinee et al., 2022).
There is a great interest in the study of BRCA1-BARD1 BRCT and RING domain and they have highlighted how some of the selected mutations play an important role in the ubiquitination mediated DNA repair mechanism and also its relevance to the prognosis of breast cancer (Choudhary et al., 2016;Choudhary et al., 2018;Kumar et al., 2022).Barua et al. 2022 have screened missense variants on the BRCT regions of both BRCA1 and BARD1 from cBioPortal database, and the pathogenicity of unknown variants was evaluated and a total of thirteen mutations were identified.About six variants of the BARD1 BRCT, along with the two wild type BRCA1 & BARD (BRCT domain) proteins were identified and subjected to MD simulations.In a very recent study, Kumar et al., 2022, attempted a massive screening of 2118 mutation of BRCA1, among which, 222 mutations were predicted to be deleterious by various in silico tools.Four mutants of BRCA1 RING domain V11G, M18K, L22S and T97R were selected for structural analysis as these four mutants fall under conserved protein-protein interaction (PPI) region with destabilizing effect and pathogenicity.In addition to these, a number of in silico mutational studies are also reported for other diseases viz., laminopathy, cancer and tuberculosis to evaluate the effect of mutations on protein's dynamics and function (Rajendran et al., 2012;Rajendran & Sethumadhavan, 2014;Rajendran, 2016;Al-Subaie & Kamaraj, 2021;Kamaraj et al., 2021).
The foregoing discussion points to a great contemporary interest in the area, both from a computational and experimental point of view.However, an exhaustive study of the mutations directed towards single mutations of the important ring domain is missing.The current study illustrates the influence of a single point mutation on BRCA1 and BARD1 RING domain to evaluate the conformational change on BRCA1-BARD1 dimer and its relation to breast cancer progression.This is the first study wherein sixteen BRCA1-BARD1 RING domain mutations were modeled to assess how these mutations modulate the structure and function of BRCA1-BARD1 heterodimer in the dynamic system by employing a range of computational tools.
A total of eleven point mutations of BRCA1 and five for BARD1 RING domain have been considered having an association with breast cancer by interfering with the E3 ligase activity and the mutants of BRCA1 that act akin to that of wildtype were modeled for the study.Thus, a total of 17 complexes, including the wild type and sixteen mutated BRCA1-BARD1 complexes were subjected to 200 ns Molecular Dynamics (MD) simulations each.Consequently, binding energies were estimated by employing MM-PBSA (Srivastava & Sastry, 2012) and MD trajectory average conformers were used to analyze each protein-protein interaction profiling.The study provides the characterization of the structural and conformational change induced by mutation and unravel the role of mutation in BRCA1-BARD1 RING domain at atomic level in the progression of breast cancer.

Retrieval of PPI complex, mutational screening and preparation of mutant complexes
Sixteen mutant variants of BRCA1-BARD1 heterodimer complexes, identified from literature, along with the wild type complex were subjected to further analysis and computations as described below.A domain identification tool, InterPro (Blum et al., 2021), was employed to discern the domains and active sites of a protein (see Figure S2).
InterPro performs a comparative analysis on the analogous protein families.The BRCA1-BARD1 complex was obtained from PDB (PDB ID: 1JM7) (Brzovic et al., 2001).In order to perform the MD simulation, UCSF chimera was used for the protein complex preparation by removing ions and solvents (Pettersen et al., 2004).The cleaned wild type (WT) structure was mutated by each one of the mutants identified from previous studies (Hashizume et al., 2001;Brzovic et al., 2001).A total of eleven mutants were chosen for BRCA1 such as V11A, I21V, I31M, C24R, C39Y, I42V, L51W, C61G, C64Y, K65R and R71G (Brzovic et al., 2001;Ruffner et al., 2001), and five mutants selected for BARD1 includes C53W, C71Y, C83R, V85M, (Stewart et al., 2018) and R99E.Among these R99E is a designed mutation found to obstruct Ub ligase activity vigorously (Densham et al., 2016).The mutant complexes were modeled by making a selected point mutation in the WT protein structure using PyMol (DeLano, 2002).The selected sixteen point mutations on BRCA1-BARD1 complex are depicted in Figure 1A.The detail structural information of BRCA1-BARD1 heterodimer is represented in Figure 1B.All the 16 mutants and WT complexes were subjected to MD simulation for a production run of 200 ns each using Gromacs5.0.4 package (Abraham et al., 2015).

Molecular dynamics simulations and binding energy calculation
All MD simulations were carried out utilizing Gromacs 5.0.4 package (Abraham et al., 2015) with CHARMM36 force field (Best et al., 2012) and simple point charges (SPC) water model was used for the water simulations (Fuhrmans et al., 2010).To eliminate the boundary effect, periodic boundary conditions (PBC) were imposed on the system.The systems were energy minimized with a steepest descent method with 5000 steps.Prior to MD simulations, all the complex systems were equilibrated for 500 ps under a constant number of atoms, volume/pressure and temperature (300 K), of the system called NVT and NPT using a modified Berendsen thermostat (Berendsen et al., 1984(Berendsen et al., , 1995)).Long-range electrostatics was computed using PME (Darden et al., 1993) and the SETTLE algorithm for the solvent molecules (Miyamoto & Kollman, 1992), while all the bonds were restricted using the LINCS technique (Hess et al., 1997).Production run was carried out for 200 ns and MD trajectories analysis was done using xmgrace tool (Turner, 2005).Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) was employed to calculate the binding free energy (BFE).The energetic values were calculated and depicted below: The following equation was used to evaluate the free energies: The term DE MM in the equations above refers to variations in molecular mechanical energy in the gas phase.Internal energy, commonly referred to as DE bonded , as well as the van der Waals and electrostatic contributions are both included in DE nonbonded .The solvation energy (both polar and nonpolar) was calculated.G polar and G non-polar together forms DG solvation .The solvent accessible surface area is typically used to estimate the surface area (SASA), i.e.G non-polar , where c is the surface tension of 2.27 kJmol À 1 nm À 2 and b is a constant having a value of 3.85 kJmol À 1 .By using the final 20 ns (i.e. 180 ns to 200 ns) of a 200 ns trajectory, MMPBSA was performed (Jamir et al., 2022).The MmPbSaStat.py and MmPbSaDecomp.pyprogram binding free energies and residue-specific contributions were calculated for all the PPI complex systems (Kumari et al., 2014).

Effect of the mutant variants on protein stability and pathogenicity prediction
I-Mutant 2.0 server was employed to observe the effects of a mutation on the structural stability of all complex systems.By using the protein structure or sequence as an input, the server applies Support Vector Machine (SVM) to forecast how the structural stability of the protein would change in response to a single point mutation (Capriotti et al., 2005).The free energy change (DDG) value is calculated, negative DDG and positive DDG indicates, decrease structural stability and increase structural stability of a protein.The pathogenicity of the mutant variants was predicted using PredictSNP (https://loschmidt.chemi.muni.cz/predictsnp1/)(Bendl et al., 2014).PredictSNP is a user-friendly web interface that allows easy access to eight common pathogenicity prediction tools.

Protein-protein interaction profiling
The average PPI complexes were analyzed for interface profiling using three servers' viz., PDBsum, PIMA and PPcheck (Laskowski et al., 2018;Mathew & Sowdhamini, 2016;Sukhwal & Sowdhamini, 2015).PDBsum provides an imagebased structural depiction such as secondary structure information of the input macromolecules and bound small molecules with a protein, protein-protein interaction and protein-DNA interaction at the atomic level.PIMA and PPcheck server utilize the inter chain interactions of the protein complex.These servers take into account all possible chain combinations in the PDB structure that has been uploaded and identify the real interactions based on a number of characteristics, including electrostatic energy, hydrogen bonds, interface residues, short contacts, total stabilizing energy, van der Waals pairs and their energy, salt bridges, etc.Both PIMA and PPcheck server quantifies the strength of a protein-protein interface at 10 Ả (PIMA) and 8 Ả (PPcheck), respectively.The residual network analysis was carried out for initial WT, average WT and the mutant complexes using NAPS (http://bioinf.iiit.ac.in/NAPS/) (Chakrabarty & Parekh, 2016) and RING (Residue Interaction Network Generator) server (Clementel et al., 2022).Both the servers were used to identify residue-residue network information of the given PPI complexes.The protein-protein complex option is chosen in NAPS server and Ca network type and unweighted edge weight options were selected.Within the lower and upper thresholds (by default, the upper threshold is 7 and the lower threshold is 0), an edge reflects the separation between two C atoms.For the long-range interaction networks, default parameters and a minimum residue separation of one were taken into consideration.The difference in the formation of total number of edges, diameter, radius, average degree, average shortest path and clustering coefficient have been measured for WT and mutant BRCA1-BARD1 complexes.
RING server was utilized to extract the important part of the network i.e. the interface of the complex, to identify the interface hub residues across PPI interface (Jebamani et al., 2023).Other than H-bond and salt bridges this server gives information on pi-pi stack, pi-cation and electrostatic interaction.For the selected PPI complexes, residue-residue networks obtained from this server in .jsonfile format was imported in Cytoscape3.9.1 and analyzed the residue interaction by utilizing the Cytohubba module by calculating the top 50 nodes based on high degree.The nodes having highest interaction with neighboring residues considered to have highest degree and referred as hub residue (Shannon et al., 2003).

Hotspot residue identification
Only a small proportion of amino acid residues at the PPI interface are mainly responsible for stabilizing energy and have significant effect on the binding-free-energy, those residues are considered to be a hotspot.(Keskin et al., 2005;Moreira et al., 2007;Caffrey et al., 2004;Lockless & Ranganathan, 1999).The conversion of those residues to alanine results in a substantial rise in the binding free energy (DDG) from 1.5 to 2.0 kcal mol À 1 .The hotspot identification was carried out across the PPI interface using publicly available servers DrugScore PPI (Kr€ uger & Gohlke, 2010), Robetta (Kim et al., 2004;Kortemme et al., 2004) and KFC2 (Darnell et al., 2008).For the wild type and the sixteen mutant complexes the average structure was extracted from last 50 ns (150-200 ns) MD trajectory and subjected to the hotspot residue prediction servers.

Principal component analysis (PCA) and free energy landscape (FEL)
To extract the functional movements for proteins, PCA is frequently utilized.Proteins' large-scale, associated collective motions are extracted using PCA to pinpoint the motions that affect the complex's primary function.Covariance matrices of the protein-protein complex simulations were constructed using backbone atom trajectory.Gromacs utilities g_covar and g_anaeig were used for analyzing the trajectories for PCA by following the standard protocol (Amadei et al., 1993).It is advantageous to characterize the mechanism of protein folding via FEL analysis (Altis et al., 2008).The 'g_sham' module was used to produce the Gibbs free energy landscapes from the two main components (PC1 and PC2) in order to capture the conformer with the lowest energy.
Protein folding dynamics can be quantitatively described for a protein structure using FEL.

Molecular dynamics (MD) simulations and MMPBSA analysis
The MD simulations results for the WT BRCA1-BARD1 and selected sixteen mutant complexes were obtained for 200 ns MD trajectories.MD trajectories were analyzed to understand the stability and fluctuations of the complexes on the impact of specific mutations.In case of the WT, the RMSD values displayed an increase between 0 to 50 ns (3-5 nm) and converged to a constant value of around 0.45 nm (Table S1).In contrast, all variants displayed a higher converging point than the WT (Figure 2A).Among all the sixteen mutants, V11A, C24R, I31M, C39Y, I42V, R71G of BRCA1 and R99E of BARD1 mutant complex structures displayed a continuous perturbation and convergence point is not clear.While remaining mutants of BRCA1 (I21V, L51W, C61G, K65R) and BARD1 (C53W, C71Y, C83R, V85M) displayed a higher converging point with more fluctuations (6 to 8 nm) than the WT (see Figure S3 and S4).The RMSD of only the RING domains of BRCA1 and BARD1 is showed in Figure S5 and S6.The WT RING domains RMSD reached a convergence with less fluctuation with the RMSD value of 0.4 nm.While in case of few mutants the RING domains attain the convergence with high RMSD value with minor local fluctuations (V11A, C24R, I31M, L51W and R71G in BRCA1 & C53W, C71Y, C83R and V85M in BARD1).MD trajectories of the proteinprotein complexes were further evaluated by Rg, SASA and the existence of a total number of intermolecular hydrogen bonds.
Radius of gyration (Rg) designates the compactness of the system (Figure 2B & Figure 3B).In case of the WT, the Rg value displayed fluctuation in the initial simulation period and then showed a decline after 100 ns production run with 1.95 Rg value.The Rg value of the mutant heterodimer complex showed similar pattern in the mutants except for V11A, C24R, I31M and I42V of BRCA1 and C71Y, V85M, R99E of BARD1 which showed an increase in the Rg value after 125 ns (Figure S7, S8).These observations signify that the mutation in these protein complexes results in a change in the compactness.The lower Rg value indicates compactness, whereas a higher Rg value indicates a lesser level of compactness, which indicates that the PPI complex is prone to have significant fluctuations.
The RMSF value of chain A (BRCA1) and chain B (BARD1) of the complex was also analyzed to understand the impact of the mutation on the residual flexibility.The average RMSF value of WT BRCA1 is 0.26 nm and WT BARD1 is 0.27 nm.The comparative analysis indicates that the mutant variants have shown relatively more local fluctuation as compared to the WT (Figure 2C & 2D; Figure 3C & 3D).The point mutations at BRCA1 display more local fluctuation in BRCA1 residues range from 21 to 82 which contain two Zinc binding site and one central helix (see Figure 2C; Figure S9 & S10), while BRCA1 mutated complexes do not display a huge fluctuation on BARD1 residues ranges 50-99 (contains Zn þ binding site and central helix) in the heterodimer complex (see Figure 2D).Among the BRCA1 mutant complex, V11A, C24R, I31M, L51W displayed relatively more fluctuation at second Zn þ binding region and central helix of BRCA1, as compared to WT (see Figure 4).The point mutations of BARD1 (Chain B) displayed more local fluctuation in the BARD1 residues (res 50-99) as compared to BRCA1 (res 21-82) in the heterodimer complex (see Figure 3D, Figure 4, Figure S11 & S12).
Furthermore, measure of the expansion of PPI complexes' volume was analyzed by estimating the SASA (solvent accessible surface area) from the MD trajectories and average SASA value of the WT was observed to be 130.20 nm.The average SASA values of the mutant variants did not show The total binding free energy (BFE) for the WT and mutated BRCA1-BARD1 complex, was estimated using the MMPBSA method.The WT binding free energy was observed to be À 284.238 kJ/mol.Although the total binding energy in WT is slightly less as compared to many other mutants, however the van der Waals energy contributes more (-676.414kJ/mol) in WT as compared to mutants (except I31M and K65R) which specify the close proximity of atoms in both the proteins during the interaction and complexation process, thus maintaining the structural veracity and constancy of the WT dimmer during the MD simulation.The electrostatics energy contribution (-813.541kJ/mol) was also observed to be higher in WT than few of the mutants (except I31M, L51W, K65R, C53W and C71Y) (see Table 1).However, high van der Waals and electrostatic energy contributions (see Table 1) can be seen in the I31M and K65R mutant complexes with high binding free energy than the WT and these two mutant variants are considered to be non-cancerous.

Domain analysis of BRCA1-BARD1
InterPro is a database which analyzes the function of a protein from the sequences and classifies them into protein families and predicts the presence of functional domains and active sites.Two domains were predicted in BRCA1 protein which includes the BRCA1 associated protein domain (res 7-110) and the zinc finger ring domain (res 24-65).The server also predicted three domains for the BARD1 protein, viz., BARD1 zinc finger domain (res 16-80), an ankyrin repeat domain (res 16-118) and a ring domain (res 1-97) (Figure S2).

Protein stability and pathogenicity prediction of mutant variant
I-Mutant predicts the stability of the protein upon a single point mutation.It was observed that the stability of the point mutation decreased in all sixteen mutant variants (see Table 2).Further, all the mutant variants were subjected to nine tools (MAPP, nsSNPAnalyzer, PANTHER, PredictSNP, PhD-SNP, Polyphen-1, Polyphen-2, SIFT and SNAP) to predict their pathogenicity.The variants which are predicted to be deleterious by three or more prediction tools have been considered as true deleterious as shown in Table 3 with the distribution of the prediction along with accuracy percentage.Based on the prediction tool, C24R, C39Y, C61G, R71G, L51W, C53W, V11A of BRCA1 mutants and C71Y, C83R, C53R, R99E and V85M of BARD1 are predicted as deleterious mutant variants by three or more servers.Among them, BRCA1 mutants C24R, C39Y, C61G, C64Y have been predicted to be deleterious by eight servers, these mutants are also associated with breast cancer progression (Brzovic et al., 2003;Ruffner et al., 2001).R71G, L51W predicted as deleterious by six servers and V11A predicted as deleterious by five servers and these three BRCA1 mutants do not impair E3 ligase activity in the experimental assay.In case of BARD1, C53W, C71Y, C83R mutant variants predicted to be deleterious by seven servers and these three mutants are found to be associated with breast cancer, V85M and R99E predicted to be deleterious by least number of servers that is three and five servers respectively.While V85M is a less studied mutant variant having association in cancer progression and R99E is a designed mutant which has been experimentally shown to be inhibiting the E3 ligase activity and involved in breast cancer progression.

Interaction profile analysis
A comparative profile analysis of the WT and selected mutants was done for the initial and the last 50 ns average structure using the PDBsum server (see Table 3, Table S1 and Table S2) which represents the interaction statistics showing the interacting residues, interface area, H-bond, salt bridges and non-bonded interactions.In case of the WT, a small reduction in the number of interacting residues after MD simulation was observed.The interface area did not show any significant change, but the number of H-bonds between the two proteins doubled from 3 to 7 after the simulation in case of WT.The number of salt bridge remained same during the simulation.In the initial stage of the complex, salt bridges were formed between ARG7 (A) and ASP117 (B), ASP96 (A) and HIS36 (B), while in the post MD structures the bridge was formed between ASP40 (A) and LYS96 (B), GLU85 (A) and ARG43 (B).The interacting residues of the two chains before and after MD simulations for the wildtype have been represented in Figure 5. with the encircled hotspot residues.
In addition to PDBsum, PIMA and PPcheck servers were utilized to analyze the interaction profile wherein the cutoff distance across interface considered as 8 Ᾰ and 10 Ᾰ respectively.From these two servers we can see that the number of interface residues was high in the final MD conformer of the WT; however, a reduction was observed in the interface residues in most of the mutants except BRCA1 mutant K65R and I31M, where an increase in the interface residues was observed in the final MD conformer.The Hbond energies, electrostatic energies, VDW energies and total stabilizing energies between BRCA1 and BARD1 were also analyzed (see Table S2).The H-bond and total stabilizing energies are found to be relatively higher in the initial (pre-MD) and final (post-MD) WT and in few mutants (BRCA1 V11A, I31M, I42V, K65R; BARD1 V85M) which does not impair ubiquitination ligase activity except in BRCA1 C64Y.The mutations which obstruct ubiquitination activity have shown lesser H-bond and total stabilizing energy of the complex.
We could see high electrostatic energies in case of WT post-MD structure as compared to all the mutants (see Table S2).Therefore, a greater number of residues involved in the electrostatic interactions between BRCA1 and BARD1 WT complex.The BRCA1 (chain A) ALys10 formed one electrostatic interaction with BARD1 BLys110 (chain B); A38Lys formed two elec.Int. with B94Lys and B96Lys; A40Asp formed three elec.Int. with B102Asp, 96Lys, 99Arg; A78Arg formed two ele.Int. with B63Leu and 64Gly; A7Arg has one electrostatic int.with B119Glu; A85Glu with B43Arg, A88Lys with B43Arg and A96Asp with B36His and B38Arg.The residue involved in interface electrostatic interaction is obtained from PIMA server.
The post-MD analysis revealed that the interface area did not show any significant change except for a slight decrease in variant C24R, but weighty changes have been observed in the profile of the average structure of the mutant variants.The number of H-bonds between the two proteins increased up to three-fold from 3 to 10 in the mutant variants after the simulation.The number of H-bonds was highest (i.e.10) in variant R99E and I42V while the other mutants maintained the H-bonds at four or more.The number of salt bridge remained the same for the WT, but a reduction was observed for the mutant variants during the simulation  with some mutants displaying complete dissociation of the salt bridges (I21V, C24R, I31M, I42V, L51W and R71G).Significant differences in the non-bonded contacts were also observed in the mutant variants of the post-MD profile.The greatest difference in the number of non-bonded contacts was observed in variant BRCA1 C61G where 67% of the contacts dissociated and this mutation is a wellknown breast cancer associated mutation (Brzovic et al., 1998).The other variants also showed similar pattern where the number of non-bonded contact reduced by about 60% (see Table S2 & S3).

Hotspot and hub residue prediction
To evaluate the importance of key interface residues across BRCA1-BARD1 complex, hotspot residues were identified using three servers-DrugScore PPI , Robetta and KFC2.Hotspot residues can provide important information about the protein function as they contribute to most of the binding free energy towards the formation of PPI complex.The prediction was performed for the WT and mutant average structure extracted from last 50 ns trajectories.The resultant values were also compared with the per-residue binding score obtained from MMPBSA where higher value of per-residue energy contribution was observed in the hotspot residues (Table S3 and S4).In case of the pre-MD initial WT, a total of thirteen residues were predicted as hotspot by two and more servers.The predicted residues are ARG-7-A, VAL-11-A, VAL-14-A, MET-18-A, ILE-21-A, PHE-79-A, LEU-82-A, ILE-89-A, TRP-34-B, LEU-47-B, MET-104-B, LEU-107-B and LEU-114-B (Table S3, whereas, the WT 50 ns average structure presented a total of fifteen hotspot residues predicted by two or more servers.The fifteen predicted hotspot residue includes: ARG-7-A, VAL-11-A, VAL-14-A, MET-18-A, ASP-40-A, ARG-78-A, LEU-82-A, ILE-89-A, TRP-34-B, LEU-44-B, LEU-47-B, TRP-91-B, MET-104-B, LEU-107-B and LEU-114-B (Table S4).Eleven common hotspots residues were observed after mapping the initial WT and average WT results.The common residues are ARG-7-A, VAL-11-A, VAL-14-A, MET-18-A, PHE-79-A, LEU-82-A, ILE-89-A, LEU-47-B, MET-104-B, LEU-107-B and LEU-114-B.Ten out of these eleven residues are also in the top twenty-five per-residue binding score obtained from MMPBSA analysis (see Figure S17).The analysis of the predicted hotspots revealed that the composition of amino acid residue leucine was the highest in all three servers.Leu constituted 34% of the total hotspots predicted by KFC, 20% in DrugScore PPI and 29% in Robetta for the post-MD complex.Leucine is involved mainly in the protein synthesis by promoting the energy metabolism.Therefore, this amino acid residue can serve as a key interface residue whose mutation may affect the binding of the two proteins.
The number of interface hotspot residues predicted for the post MD mutant complex variants, were comparatively lesser than the post-MD WT where one common hotspot residue (TRP34B in BARD1) was observed in fifteen mutant variants, except I21V (see Table 5).The BARD1 TRP34 residue formed H-bond with BRCA1 LEU6 in WT as well as with thirteen mutant complexes (except BRCA1 V11A and C39Y, I21V BRCA1-BARD1 Decrease À 1.49 2.
BRCA1-BARD1 C53W BARD1) of MD simulated structures and this could be considered as one of the important hotspot residues involved in the formation of complex that may maintain the structural integrity of heterodimer complex.The detailed Hbond and salt-bridge information are shown in Table S3-S21.The loss of hotspot residues might obstruct the structural integrity and which in turn destabilize the BRCA1-BARD1 heterodimer thus failing to carry out their normal activities like ubiquitination and DNA repair.
Network analysis was carried out to identify interface hub residue by utilizing RING sever followed by the Cytoscape tool.The residues Trp34B, Leu114B of BARD1 and BRCA1 Ile89A, Glu85A were found to be the hub residues based on its maximum number of interaction (degree) with neighboring residues.Among these Trp34B, Leu114B and Ile89A (see Table 6 and 7) have also been predicted as interface hotspot residue in the WT and in most of the mutants.As an example, the hub residue result of C24R BARD1 complex is shown in Figure 6 and for WT and other mutant complexes see Figure S18 to Figure S34.
Further residue-residue network analysis from NAPS server shows that the initial WT (pre-MD) complex has a greater number of edges (interaction pathways) as compared to post-MD WT and mutant complexes.The PPI complex having a greater number of edges represents more interactions, high diameter and less compactness.Some of the post MD mutant average structure of BRCA1 (V11A, I21V, I31M, C39Y, C64Y, L51W, R71Y) and BARD1 (C83R, V85M) have shown a greater number of edges as compared to post-MD WT, thus these mutants display high interactions, more fluctuation and less compactness in the PPI complexes.The number of edges, diameter, radius, average degree, average shortest path and cluster coefficient are presented in Table 8.

Principal component analysis (PCA) and free energy landscape (FEL)
The PCA of the MD simulation trajectories was performed using GROMACS.PCA display the combined motion in MD trajectories, through the eigenvectors of the mass-weighted covariance matrix (C) of the atomic positional fluctuations.
The PCA results of the protein-protein complex and the conformational differences were visualized for each trajectory data-set by plotting in accordance with PC1 and PC2 (see Figure S35 & S36).PCA helps in the better understanding of the protein function by analyzing the comprehensive movement of a protein from slow collective movement to fast local movement.The eigenvector delineates the protein movement atoms whereas the eigenvector values give the overall contribution of protein atoms group movement.The highest eigenvector value represents the most dominant collective motions.PCA analysis was performed for 200 ns MD trajectories of WT and sixteen mutant variants of BRCA1-BARD1 protein complex.The projection of the WT and mutant variant shows significant difference.The representation of clusters for WT, is localized between À 6 to 6 nm at Y axis and À 5 to 8 nm at X axis, while the mutant variants were localized between more than À 6 to 5 nm.Most of the Table 3.The pathogenicity prediction obtained using PredictSNP prediction tools.The mutant variant predicted as deleterious by three or more prediction tools has been evaluated as detrimental mutants.Variant  Table 5. Hotspot residues predicted for the mutant variants and wildtype from the initial structure and post-MD average structure (from last 50 ns trajectory) using three servers i.e.KFC2, Robetta and DrugScore PPI .The highlighted row indicates the residue common among all complexes.
mutant variant showed highest deviation from the WT, represents a significant difference in the phase space coverage.
The PC analysis displayed that, the projection spanned larger ranges in mutants than the WT system, which shows the changes in the conformation of the complex by the mutations (Figure S35 & S36).The backbone atoms were considered for generating gibbs free energy landscape for BRCA1 and BARD1 mutant complex systems (see Figures 7 and 8), red color indicating the lowest energy conformer.In case of wild type PCA plot, clear collective motions can be observed displaying three local basins containing local minima structures.However, mutant variant, indicating that, mutation affected the overall conformational stability of the system, with distorted collective motions in the phase space coverage.

Conformation and structural changes of the mutant complexes
The conformation and structural changes were visualized for all the average structures obtained from MD trajectory (last 50 ns) in UCSF chimera (see Figure 4).In case of BRCA1 mutants (V11A, I21V, I42V, R71G, I31M, L51W, K65R) that does not interfere with Ub ligase activity, we could not see huge fluctuation in the central helix of BRCA1 which act as the interface for UBE2D3 binding.In case of designed BARD1 mutant R99E which is known to obstruction of Ub ligase activity abruptly (Densham et al., 2016), interestingly we observed a loss of central helix in this variant at the end of the simulation (see Figure 4D, in gray).The loss of central helix could be the possible reason of impairing E3 ligase activity and hamper ubiquitylation of histone H2A of nucleosome core particle (NCP), because the central helix and Zn 2þ binding region of both the BRCA1-BARD1 along with UBE2D3 tether to the NCP for H2A ubiquitylation (Witus et al., 2021;Hu et al., 2021).Therefore, structural changes in these regions may affect the functional activity.Both the BRCA1 Zn 2þ coordinate binding residue mutations (C24R, C39Y, C61G, C64Y) and central helix of BRCA1-BARD1 heterodimer complexes have shown some fluctuations and conformational changes in BRCA1 variant complexes (see Figure 4 & Figure S40) as compared to wild type (see Figure S37).In these mutations, the central helix of BRCA1 and BARD1 moved much away from the initial position, shifted to a different direction and went inward in some cases which may hamper the binding of UBE2D3 through BRCA1 and binding of BRCA1-BARD1 dimer to nucleosome (H2A) thus impairing ubiquitination ligase activity.
Interestingly in case of BRCA1 C64Y which is a common inherited mutant, one of the Zn 2þ coordinate binding residues i.e.C39 position change the conformation from coil to a small b strand (see Figure S40D), which again could be the possible reason for impairing the UBE2D3 binding and interfere with Ub ligase.The DSSP graph displays the secondary structure transition along the 200 ns MD simulation (see Figure S38 & S39).The plots display the number of residues involved in the formation of a particular secondary structure

Discussion
MD analysis reveals that the mutant complexes display overall local perturbation and converged at high RMSD.The RMSF analysis confer structural flexibility of mutated complex mainly at both the Zn 2þ binding site and central helix of BRCA1 and BARD1 as compared to wild type, which may possibly influence the functional activity of the heterodimer.This trend was observed in all mutants with variation in the range of fluctuations (Figure S9-S12).In an experimental study, the nearby residues of the two zinc binding sites of BRCA1 RING domain interact with E2 (UBE2D3) (see Figure S1) allowing the heterodimer to coordinate with substrate and transfer E2 simultaneously thus leading to the E3 ubiquitin ligase activity and this process is facilitate by RING domain (Dai et al., 2021;Hu et al., 2021).The current study emphasizes the high structural fluctuation in the E2 binding site of mutated BRCA1 may impair the binding of E2 Ub, thereby inhibiting the ubiquitin E3 ligase activity leading to the breast cancer progression.The prediction from I-Mutant server also indicates that the mutation decreases the stability of the protein.
In cancer-associated mutations of BRCA1 residues at Zn 2þ binding site (C39Y, C61G, C64Y), were observed to eliminate E3 ligase activity for all known substrates by disturbing the structural integrity of the RING that probably destabilize BRCA1 and BARD1 in cells (Hashizume et al., 2001;Brzovic et al., 2001).While in another breast cancer mutation study (Stewart et al., 2018), BARD1 Zn 2þ binding residues (C53W, C71Y and C83R) mutation interfere with nucleosomes binding and hamper the ubiquitination of histone H2A.However, these mutation does not affect the ability to form heterodimer complex with BRCA1 (Stewart et al., 2018) and we assume that the complex remained together probably due to some of the identified interface hotspot residues in this study (Trp34, Leu114 in BARD1 and BRCA1 Ile89) which maintained the high degree (interaction) during simulation for WT and in most of the mutant complexes.And, interestingly, some of the Zn 2þ coordinate binding mutations of both BRCA1 (C24R, C39Y, C61G) and BARD1 (C53W and C83R) complexes display a significant reduction in the number of H-bond, salt bridges and non-bonded contact during simulations (see Table 4 and Figure S14 and S15) which might obstruct the structural integrity by destabilizing the heterodimer.Ruffner et al. (2001) reported that some of the BRCA1 mutation viz., V11A, I21V, I31M, I42V and R71G does not disturb ligase activity and act as the wild type, while combination of L51W and K65R (Stewart et al., 2017) mutant complex rescues the cancer associated variants (C61G and C64Y) (Brzovic et al., 2003;Ruffner et al., 2001).In the current study H-bonds and non-bonded contact between the BRCA1 and BARD1 increased at the end of the simulation for L51W and K65R mutants that also display overall stable convergence in RMSD as that of WT.Ruffner et al. (2001) observed, mutations of BRCA1 Zn 2þ coordinate binding residues (C39Y, C61G, C64Y) obstructs the ligase activity and display as a ligase dead mutant of BRCA1.Another study (Brzovic et al., 2003) showed C61G mutation disrupts Zn 2þ binding necessary for stability of the RING structure subsequently abolishes interaction with E2 conjugating enzymes as well as with BARD1.It was observed that C61G mutant complex displays a significant reduction in the number of H-bonds and nonbonded contacts.NAPS server analysis also displayed a lesser number of edges in the C61G complex signifying a reduction in the number of interaction (see Table 4 and Figure S15).The loss of significant H-bonds and non-bonded contacts in C61G complex may have affected the complex in a way that interfered the ability of BRCA1 to bind to UBE2D3 enzymes, which in turn interfered with the ubiquitin ligase activity.However interestingly, C61G complex structure display a stable RMSD and relatively high compactness (Rg) with less residual fluctuation observed from RMSF.To the best of our knowledge the current study is the first MD simulation study carried out on several mutations of the RING domain of BRCA1-BARD1 dimer.There are a few similar mutational molecular modeling studies focused on the BRCT domain of BRCA1 and BARD1 and analyzed the differences in the flexibility of mutated complexes and wild-type proteins.(Barua et al., 2022;Choudhary et al., 2016;Choudhary et al.,2018).

Conclusions
This is the first attempt where a comprehensive analysis on the effect of sixteen mutation was carried out on BRCA1-BARD1 RING dimer to assess the conformational change and interface profiling.This computational study provides insights into each mutation-induced structural and conformational changes and unveils a promising way for the BRCA1-BARD1 RING domain mutations' role in the breast cancer progression and susceptibility.The mutant complexes have more local perturbation with a high residual fluctuation as compared to WT according to MD simulation results.High structural fluctuation in the E2 binding site of BRCA1 in the mutant variants have been observed in this study, this might be a possible reason for interfering the binding of E2 Ub to BRCA1 and inhibiting the ubiquitin ligase activity leading to breast cancer.The BRCA1 mutations (V11A, I42V, R71G, I31M) which do not interfere E3 ligase activity, have shown more H-bond formation and non-bonded contacts at BRCA1-BARD1 interacting interface during the simulation.While the mutations of BRCA1 (C61G, C64Y, C39Y, C24R and BARD1 C53W, C83R) which hamper the ubiquitination process and are associated with breast cancer, display a reduction in the number of significant H-bond as compared to WT at the end of the simulation (with a few exceptions BARD1 C71Y and R99E).
According to some studies, although a few mutations hamper the E3 ligase activity, the BRCA1-BARD1 dimer remains intact.The reason might be the interfaces hotspot and hub residues identified in this study, may play a significant role in holding the complex together.

Figure 1 .
Figure 1.Crystal structure of the BRCA1-BARD1 complex (PDB ID: 1JM7) (A) BRCA1 (chain A) is represented purple color and BARD1 (chain B) is represented in green color respectively.(B) Representation of BRCA1-BARD1 complex binds with E2 ubiquitin through BRCA1 RING interface (PDB ID: 7JZV), the blue highlighted residues are involved in Zn 2þ coordinate binding site I, red highlighted residues are Zn 2þ coordinate binding site II.

Figure 2 .
Figure 2. Molecular dynamics simulations analysis for 200 ns trajectory for the wild type (WT) BRCA1-BARD1 complex and mutant variants of BRCA1 (A) RMSD plot (B) Rg (C) RMSF of chain A (D) RMSF of chain B (E) SASA (F)Hydrogen bond (H-bond) modulations, along the 200 ns MD trajectory.

Figure 3 .
Figure 3. Molecular dynamics simulations analysis for 200 ns trajectory for the wild type (WT) BRCA1-BARD1 complex and mutant variants of BARD1 (A) RMSD plot (B) Rg (C) RMSF of chain A (D) RMSF of chain B (E) SASA (F) Hydrogen bond (H-bond) modulations, along the 200 ns MD trajectory.

Figure 5 .
Figure 5.The interacting residues of BRCA1 (Chain A) and BARD1 (Chain B) (WT) of human breast cancer obtained from PDBsum server.(A) The interacting residues before the MD simulation (B) The interacting residues after the MD simulation using the average 50 ns trajectory.The encircled residues represent the predicted hotspot from three servers (KFC, Robetta and DrugScore PPI and validated from per-residue energy decomposition analysis).

Figure 6 .
Figure 6.Interface hub residue analysis depicted for C24R BRCA1 mutation: (a) Last 50 ns average complex structure was submitted to RING server.Selected the edge type check box 'inter-chain' to visualize the interchain residue-residue interaction network between BRCA1 (Chain A) and BARD1 (Chain B).Downloaded the network in '. json' file format.(b) Imported the network.jsonfile of each complex in Cytoscape and visualized the top 50 nodes based on highest degree (interaction) using the Cytohubba module.The nodes having highest degree ranked one and depicted in dark red in color.The dotted rectangle box displays the interface hub residue having a greater number of interaction or high degree.The color range, dark red to yellow signify highest to lowest degree of a residue.(c) Excel sheet containing list of nodes with its corresponding score presented in ascending order (highest to lowest).The star mark signifies the interface residue having high degree.
Predictor of human Deleterious Single Nucleotide Polymorphisms, MAPP: Multivariate Analysis of Protein Polymorphism, nsSNPAnalyzer: Nonsynonymous single nucleotide polymorphisms Analyzer, SIFT: Sorting Intolerant From Tolerant, SNAP: Screening for Non-Acceptable Polymorphisms, PANTHER: Protein ANalysisTHrough Evolutionary Relationship.
and transition during the simulation period.Overall, the trend of conformational changes and fluctuation is almost similar among the mutants that impair ubiquitin ligase activity and similar pattern was observed among the mutant complexes that does not impair ubiquitin ligase activity.

Table 1 .
MMPBSA binding free energy calculation of wild type and mutant BRCA1-BARD1 complexes.Energies are in kJ/mol.

Table 4 .
The interface profile information obtained from PDBsum server for the comparison of the wild type and mutated BRCA1-BARD1 complexes from the average complex structure extracted from the last 50 ns trajectory of 200 ns MD simulation.

Table 6 .
Degree of Hub residue in wild type and BRCA1 mutant complexes analyzed for top 50 nodes ranked by degree.

Table 7 .
Degree of Hub residue in wild type and BARD1 mutant complexes analyzed for to 50 nodes ranked by degree.

Table 8 .
Residue-residue network analysis using NAPS server.The final structure is the average structure extracted from last 20 ns (80-100 ns) MD trajectories.C a atom of an amino acid residue is considered as node and an edge is drawn if the C a -C a distance between a pair of residues is within a threshold distance, Rc (�7 Ả).Ca network type and unweighted edge weight, default parameters for the long-range interaction networks and minimum residue separation of 1 were considered.The final structure designates the average structure extracted from last 20 ns (80-100 ns) from MD trajectory.Spatial distance or interaction energy between amino acid residues.Radius and Diameter: Information regarding the shape of the protein.Average shortest path: Average of the degree of its immediate neighbors.Cluster Coefficient: Ratio of number of connected neighbors of a node to total number of connections possible between the neighbors.