# Theoretical Investigation of the D83V Mutation within the Myocyte-Specific Enhancer Factor-2 Beta and Its Role in Cancer

Authors: Oleksandr Yakovenko; Ryan Morin; Ganna Vashchenko; Steven J. M. Jones
Journal: Journal of Theoretical Chemistry
 (2013)

---

## Abstract

The D83V mutation in the myocyte-specific enhancer factor-2 beta (MEF2B) gene is frequently observed in lymphomas. Surprisingly, this apparent gain-of-function mutation is within a protein that is involved in the promotion of apoptosis in B cells. To investigate the oncogenic effects of this alteration and explain its predominance over other known loss-of-function mutations of MEF2B, we propose a hypothesis that this mutation influences the dynamic folding of the C-terminal loop of the N-terminal domain of MEF2B. According to our hypothesis, the mutation allows MEF2B to bind promiscuously to a wider variety of gene promoters. A large set of molecular dynamic simulations (MD) was conducted to investigate the effects of D83V mutation in silico and support the hypothesis.


---

## Body

## 1. Introduction

The myocyte-specific enhancer factor-2 (MEF2) family of transcription factors plays an essential role in myogenesis [1, 2] and the regulation of the proapoptotic gene Nur77 [3–6]. This family is composed of four members, MEF2A, -B, -C, and -D, which normally function either as homo- or heterodimers [7]. All of these factors recognize the DNA consensus sequence C/TTA(A/T)4TAG/A, which is present within many gene promoters (see [8, 9] and references therein), and consist of a MADS-box containing a N-terminal domain responsible for DNA recognition and dimerization and an as yet uncharacterized C-terminal domain [10]. The sequence differences between the N-terminal domain of MEF2 forms are mainly within the MEF2S subdomain (residues 58–95) of the N-terminal domain which is responsible for the determining of DNA binding specificity [11–18].Recently,MEF2B was found to be altered in many B-cell-derived lymphomas [19]. Morin et al. [19] and Youn and Liu [20] proposed a model of MEF2B involvement in the development of lymphomas through Ca2+ dependent signaling downstream of the B-cell receptor. In normal germinal centre B-cells MEF2B binds to the promoter of Nur77, a proapoptotic gene that plays a crucial role in TCR-mediated apoptosis [21–23]. However, MEF2B is kept transcriptionally inactive by the suppressor Cabin1 which binds to its N-terminal domain preventing the transcription of Nur77. This repression is calcium-dependent [3, 4]: when Ca2+ enters the cytoplasm in response to B-cell receptor activation, it folds the intrinsically unstructured Calmodulin, which becomes a more preferred partner of Cabin1. Calmodulin in complex with Ca2+ releases Cabin1 from its complex with MEF2B, allowing MEF2B to initiate Nur77 transcription and recruit other activating factors such as HDAC7, HDAC9, and CBP/p300 [24]. We also note that interactions discussed above are not the only way to control MEF2B activity; it is also regulated by a significant number of posttranslational modifications [25] and through direct interactions of other regulatory proteins with the uncharacterized C-terminal domain of MEF2 [10, 26].Surprisingly, the recurrent D83V mutation ofMEF2B is the most frequent alteration of this protein in lymphoma cells. In contrast, the other characterized MEF2B alternations (K4, R24, E77, Y69) in lymphoma are all likely loss-of-function events [8, 19, 20, 27] and prevent apoptosis through constitutive inactivation of the Nur77 pathway. However, a simple loss-of-function model cannot explain the specific and recurrent D83V mutation which suggests a novel neomorphic role for this mutated protein form in the regulation of genes expression. This gained functionality may be of two modalities: it either perturbs the inherent mechanism of MEF2B activity or causes the binding of some other cofactor(s) which then, indirectly, alters the function. The hypothesis about indirect altering of MEF2B functions by its binding to a cofactor is less probable though. The hypothetical new functionality of the complex formed due to D83V mutation is most likely related to the transcriptional function of MEF2B; otherwise, it is unclear how a transcriptional factor would contribute to oncogenesis apart from DNA. Moreover, the new transcriptional function of mutated MEF2B in complex with a cofactor, whatever it would be, has to create a new interaction with DNA as it follows from the domination of the recurrent D83V over all other mutations in lymphoma [19]. Otherwise it is statistically impossible for a tumor to collect by chance the destructive recurrent mutation in the one particular residue whilst there are numerous other residues (R3, K4, K5, R15, R17, K23, R24, K30, K31, etc.) whose mutations are known to diminish (by different degree) the interactions of MEF2B and DNA [8]. However, the new hypothetical interaction between DNA and a complex of MEF2B with a hypothetical cofactor would be strongly mitigated by competition with the wild selectivity of MEF2B alone, which is known to be sufficiently high for keeping MEF2B constitutively bound to the Nur77 promoter [6]. To resolve the competition and allow a new functionality to truly arise, one has to consider the decreased selectivity of the mutant protein, which is actually in favor of perturbing the inherent mechanism of MEF2B activity rendering a new cofactor binding fully optional. Moreover, this hypothesis is not consistent with available structural data—though the domain is altered by D83V mutation, it is absolutely not clear how the distant alternation can perturb the interactions (Figure 1). MEF2B binds its cofactors in the same manner to the same site within the N-terminal domain in all resolved crystal structures: of DNA-MEF2 complexes alone [10, 28] and in complex with the activators CBP/p300 [29], HDAC9 [30], and the repressor Cabin1 [31]. Particularly the only essential difference between the last two structural models is the interactions of Phe150 residue in HDAC9 peptide and Ala2175 in Cabin1 located at the same corresponding 3D positions of spatial complex structure, which results in relative shift of the peptides by approximately 1.5 Å along binding site on N-terminal domain of MEF2. On the other hand, none of resolved structures show evidence that conformational changes induced by the D83V mutation can affect transcriptional regulation. The only region of MEF2 where the published studies are not in perfect agreement is the s3 loop at the end of the N-terminal domain (residues 79–86, here and subsequently we will use the previously described MEF2 structural component numeration [32]), within which the recurrent D83V mutation resides (see Figure 1). In the study of Santelli and Richmond [28] this fragment was removed, whereas Han et al. [30, 31] reported surprisingly large conformational variation within this region by comparing two very similar crystal structures. Only Wu et al. [32] reported high quality structural information within this region, although it was obtained by the artificial packing of two asymmetrical structures per lattice, whereas the results of a NMR study of MEF2-DNA complexes suggested that this region should not have any regular secondary structure at all [10]. Though these structural differences were obtained in studies of different forms of MEF2, mutagenesis studies indicate that specific functionality resides in this interdomain region [8, 27], where it was shown that integrity of this region between the N- and C-terminal domains is required for successful transcription initiation.Figure 1
A symmetrical MEF2B dimer (shown as a multicolor cartoon) is bound to asymmetrical DNA (indicated in violet). The asymmetry of a complex is further enhanced with binding of a HDAC peptide ligand (indicated in dark blue) to the N-terminal domain of MEF2B. The D83 residue of MEF2B is shown as sticks.In this study, we used molecular dynamic simulations (MD) to theoretically investigate the structural regulation of MEF2B activity and the influence of D83V mutation upon it. MD provides an atomic-level simulation of interactions in complex multimacromolecular systems and allows a numerical evaluation of corresponding free energies. Free energies of folding and binding were evaluated with umbrella sampling [33] and the Bennett acceptance ratio (BAR) [34, 35] methods. The umbrella sampling method is rigorous but computationally intensive approach to perform numerical estimation of absolute free energy from the trajectory of a process [36–38]. The umbrella sampling method has been successfully used in numerous studies of systems of comparable size: studying ions passing through channels [39], the affinity of polyamines to DNA [40], and the probability of processes in lipid membranes (see [41] and references there in). Here we applied umbrella sampling along multiple trajectories to reliably evaluate the folding free energy of the MEF2B interdomain region. Estimations of DNA binding affinities, however, are too computationally expensive for multiple umbrella samplings so the BAR method was used to predict differences in binding free energies of DNA to wild-type (wt) or mutant (mt) dimers of MEF2B.


## 2. Methods

Molecular dynamic simulations were carried using methods similar to those previously described [36], with only several changes. To facilitate the sampling all DNA molecules were restrained to their optimal geometry throughout the simulations. Multiple reaction trajectories were sampled with the umbrella sampling method in order to get convergence in estimations of the potential of mean force (PMF). The BAR method was used for relative free energy estimations of mt and wt MEF2B dimer affinity to DNA where the computational complexity of multiple umbrella samplings would be unfeasible. All simulations were carried out with the GROMACS package version 4.5.5 [42–45]. Analysis of the pulling simulations was performed with the weighted histogram analysis method (WHAM) [46] implemented as a g_wham routine [47] from the GROMACS package; BAR analysis was carried out with g_bar routine from the GROMACS package. For the visual inspections of MD trajectories the VMD viewer [48] was used.Several different models were used in this study.(i)
From the resolved MEF2-DNA complex structures, the most accurate one determined by X-ray diffraction [32] was chosen as our reference structure. In this model (PDB accession 3KOV) the s3 fragment of interest is folded and resolved at a high definition with high quality electron density maps. Because this chosen model was a DNA-MEF2A complex, a complex reflecting the MEF2B form was obtained using SWISS-MODEL homology modeling and 3KOV as a template [49, 50]. In order to diminish the terminal effects of solvating hydrophobic DNA in a polar water environment, each chain was modified through the addition of three flanking GC pairs at both ends plus an unpaired adenine residue on the 5′ terminus. These modifications are known not to affect the binding of MEF2B to DNA [32]. The final complex model of flanked DNA bound to wt MEF2B dimer is referred to as the N model.(ii)
Models of triple complexes with resolved partners, DNA-MEF2B2-Cabin1 peptide (C model) and DNA-MEF2B2-HDAC9 peptide (D model), were created from corresponding crystal structures [31] (PDB accession 1N6J) and [30] (PDB accession 1TQE) by superposition and substitution of the MEF2-DNA fragment with those of N model. In addition, the three last C-terminal residues of Cabin1 peptide in the C model were deleted to obtain a peptide of the same length as in the D model.(iii)
It is known from sequencing studies that D83 residue can be mutated in lymphomas not only into D83V (observed in 20 out of 27 cases) but also to D83A (6 out of 27 cases) and to D83G (1 out of 27 cases) [19]. Therefore, three models of the mutant MEFB protein with D83V (V model), D83A (A model), and D83G (G model) were prepared from the N model via manual replacement of the D83 residue and while maintaining the geometry of the sidechain as allowable without collisions with other atoms, the coordinates of backbone and beta carbon atoms (except of G model) were not changed.The MEF2B homodimer itself is a symmetric complex, but due to the asymmetric nature of DNA the triple complex is asymmetrical. This asymmetry becomes more pronounced after binding of the effector molecule to N-terminal domain (Figure1). Because of this asymmetry residues within the two MEF2B monomers are not under identical conditions and are expected to behave differently. In order to distinguish them, the monomer for which the a3 helix is in proximity to the Phe150 of the HDAC peptide in the D model (and the Ala2175 of the Cabin1 peptide in the C model) is referred to as the first monomer and the other as the second. The structural fragments of the MEF2B monomer are denoted with a slash index; for example, s3 sheet of the second monomer is denoted as s3/2, and the a3 helix of the first monomer is designated as a3/1.Prior to the simulations, all models (N, C, D, V, A, and G) were placed in triclinic box of simple point charge (SPC) water [51] to which 100 mM NaCl equivalent was added including neutralizing counterions. All the six models were approximately the same size, 200 nm3, and consisted of approximately 6.8×104 atoms each. Periodic boundaries were applied in all directions. The N- and C-termini of all protein molecules were ionized. All other titratable nucleic and amino acids were assigned to their canonical state at physiological pH. Energy terms from the GROMOS96 43a1 parameters set [52] were applied to all molecular species in the system. In all calculations short-range electrostatic interactions were cut off at 0.9 nm and van der Waals at 1.4 nm. Long-range electrostatic interactions were modeled with the particle mesh Ewald (PME) algorithm [53].To adjust the localization of solvent molecules and ions around the protein, the whole system was relaxed by l-bfgs minimization [54] and 20 picoseconds (ps) of molecular dynamic simulation with restrained positions of protein and DNA heavy atoms was conducted. The following fast equilibration consisted of two stages. The first stage involved simulation for 50 ps under a constant volume (NVT) ensemble. Following NVT equilibration, 100 ps of constant pressure (NPT) equilibration was performed. The weak couplings method [55] was used to maintain pressure isotropically at 1.0 bar and the temperature constant at 310 K; protein with DNA and other atoms were coupled to separate temperature coupling baths. All subsequent productive runs were performed with the more accurate Nose-Hoover thermostat [56, 57] with a temperature coupling time constant of 0.1 ps and the Parrinello-Rahman barostat [58, 59] with a pressure coupling time constant of 1.0 ps. This combination of thermostat and barostat ensured that a true NPT ensemble was sampled.In order to measure the stability of s3/1 and s3/2 folding the corresponding melting experiments were carried out for all models. In all cases the joined fragment of residues 61–74 of both MEF2B monomers (a2 helices) was used as the anchor. Residues 76–92 of either the first or second monomer were selected as a mobile fragment. The mobile fragment was pulled away from the anchor fragment in order to unfold the s3 loop using a standard artificial spring force of 1000 kJ mol−1 nm−2 and a pull rate of 0.001 nm ps−1 (0.01 Å ps−1). The coordinates of the atoms of the anchor fragment were not constrained artificially as the fragment remains stable with nonbonded interactions. Pulling forces were applied along all axes. The overall time of the simulation was 2.5 ns. An artificial separation between the center of mass (COM) of the anchor and the mobile fragments was used as a reaction coordinate. During the pulling simulation the s3 loops melt first followed by the a3 helices. In this experiment, however, we integrated a trajectory over 2.1–3.2 nm along the reaction coordinate until the s3 loops melted completely, but the a3 helices remained almost intact. From the obtained trajectories, snapshots were taken to generate the starting configurations for the umbrella sampling with 13 windows that covers the range of reaction coordinate ξ (0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, and 1.2 nm of artificial COM separation). In all cases, 12 ns of productive MD with fixed COM separation per window were performed but the first 2 ns were skipped during analysis to allow for equilibration.The estimations of free energies of binding with equilibrium methods, one of which is umbrella sampling, are crucially dependent on applied pulling rate and strength of pulling force. In contrast, nonequilibrium methods suffer from the inaccurate distribution of low energy. To cope with these variations in estimation, a set of multiple paths was generated with pulling under slightly different conditions and independently sampled with umbrella sampling methods for the adequate distribution of high- and low-energy trajectories. The group was then averaged in a similar way as with the Jarzynski nonequilibrium method [60]. To create a source of different trajectories, different temperatures of pulling and sampling were used: the melting trajectories were generated and sampled at 300 K, 305 K, 310 K, and 315 K for each model in this study (in total 16 samplings per COM separation per model). The resulting estimations were averaged to converge upon the most probable path-independent free energy of unfolding.To measure the differences in DNA binding between mt and wt MEF2B, the free energy of alchemical transformation of 83th residueD83↔V83 was evaluated with the BAR method. As an initial position for the alchemical simulations conformation with both melted s3 loops of mt (model V) and wt (model N) MEF2B was used. The alchemical perturbation was performed in both directions: from D model into V model and vice versa. In order to preserve zero-charged state during the perturbation, the alchemical perturbations were coupled to perturbation of corresponding counterions into methane molecules. The perturbation consisted of 21 intermediate stages for each direction with λ = 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00. Each alchemical λ was sampled for 120 ns of productive MD the first 2 ns being skipped to allow equilibration.


## 3. Results

We tested our hypothesis that the D83V mutation occurring within the s3 loops of MEF2B can inherently alter the function of the MEF2B dimer without implicating any other interacting proteins. The transcriptional activity of the MEF2B dimer requires coordination between the C- and N-terminal domains that theoretically could be achieved through the constraining of the C-terminal domains at the proper positions in space, which in turn is facilitated by the folding of interdomain s3 loop only after the DNA binding sequence is recognized by the N-terminal domains. Otherwise, in solution, the s3 loop is unfolded disrupting any interdomain communication due to entropy effects. This is supported by the observation that the s3 loop only is folded in structures which are bound to DNA. It can be proposed that such dynamic folding in response to sequence-specific DNA recognition (and presumably further enforced by cofactor binding) prevents transcriptional activation at inappropriate DNA sequences. Likewise, neofunctionality, which is introduced by the D83V mutation, can also be proposed to interfere with DNA site-specific interactions. In order to support the hypothesis, the free energies of the s3/1 and s3/2 sheets versus their unstructured forms were calculated for C, D, N, V, A, and G models and the free energy of the MEF2B binding to DNA was computationally evaluated for the mutant (V) and wild-type (N) models.


### 3.1. D83V, D83G, and D83A Mutations Are All Predicted to Destabilize the Folding of S3 Loops

The stability of the s3/1 and s3/2 loops was measured as described in the methods section—by averaging of the PMFs calculated for a set of different unfolding trajectories. The four different melting trajectories were generated through the pulling of the s3/1 or s3/2 fragment of MEF2B away from the joined a2 fragments of both monomers at temperatures 300 K, 305 K, 310 K, and 315 K. Such an approach provided a set of different energetic paths of unfolding as described by the force at reaction coordinate plots (Figure2). The behavior of the loops is however similar in all cases: because of the applied external force in each trajectory, the s3 loop melts first followed by stretching of the a3 helix (an example of the process is shown in Figure 3) which is explaining why most paths have two peaks and a valley between them on the force at reaction coordinate plots.Force strength along the reaction coordinates under different temperatures: 300 K (blue), 305 K (cyan), 310 K (green), and 315 K (red).
(a)
Unfolding s3/1 loop
(b)
Unfolding of s3/2 loopTrajectory of s3 loops melting. (a) Melting of s3/1 loop of V model at 310 K and reaction coordinateξ = 2.1, 2.6, 3.4, and 4.5 nm. (b) Melting of s3/2 loop of N-model at 310 K and reaction coordinate ξ = 2.1, 2.4, 3.4, and 4.5 nm.
(a)
(b)In order to study the stability of the s3 loops further in the different models, each of the generated unfolding trajectories was than sampled with the umbrella sampling method again at four slightly different temperatures 305 K, 310 K, 315 K, and 320 K, which resulted in 16 different individual PMFs per model per loop (see Figure S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2013/313419). The upper limit of PMF integration was set to the “valley” on the force at reaction coordinate diagram, the region ξ = 3.3 nm of artificial separation—at the point where s3/1 loops are completely melted but the a3/1 helixes remain almost intact in all models (the lower limit was set at ξ = 2.1 nm—the initial distance from the source pdb file).Although each individual trajectory shows a strong variation in the value of PMF, their average (Figure4) is appropriate for ranking and consistently different between models. Because the melting of the s3 loop and the stretching of the a3 helix are two events that cannot be clearly separated during the simulations, the PMF is still growing along the reaction coordinate at the end of the integration. To prove that the estimations are still valid, a relative PMF growing diagram for the models was built (Figure 5): even if the exact upper limit of the integration would be shifted, the difference between averaged PMFs would remain constant along the “valley” region (3.0–3.4 nm of artificial COM separation).PMFs along the reaction coordinate of s3 loops unfolding for D model (blue), C model (cyan), N model (green), V model (red), A model (orange), and G model (yellow).
(a)
PFM of s3/1 loops melting
(b)
PMF of s3/2 loops meltingAveraged PMFs of s3 loops melting for D model (blue), C model (cyan), N model (green), V model (red), A model (orange), and G model (yellow) and their difference from N model.
(a)
(b)
(c)
(d)The results of these experiments show that the activating of HDAC protein binding to theMEF2B2-DNA complex (model D) dramatically (by more than 20 KJ mol-1) increases stability of the active s3/2 fold, an energy change which is not observed with the binding of the Cabin1 suppressor (model C). In accordance with our hypothesis, the dynamic folding MEF2B protein is destabilized by the D83 mutations. The s3/1 fold is significantly destabilized by approximately 10 KJ mol-1 in the D83V mutant (model V). It is also interesting that the mutants exhibit different trends: the rare D83G mutation (model G) shows a similar yet weaker pattern of destabilization as D83V, whereas, the D83A change surprisingly has converse effects—stabilizing the active folding of s3/1 but destabilizing that of s3/2 by approximately 10 KJ mol-1.


### 3.2. The D83V Mutation Is Predicted to Improve Nonspecific Binding of MEF2B to DNA

Because the destabilization of MEF2B folding is actually a loss-of-function event, a change in the spectrum of DNA sequences recognized could explain apparent neofunctional nature of the mutations (otherwise its high recurrence in tumours remains unexplained). Therefore, the suggested mechanism is not the mitigation of the protein affinity to DNA (which would be trivial loss-of-function event) but to strengthen nonsequence-specific affinity to polyanionic frame of DNA. To test this, an alchemical mutation of MEF2 in complex with a DNA double helix was performed and the affinity difference between mutant and wild-type MEF2B binding was evaluated. Though the accurate direct estimation of free energy of DNA binding to mutant (ΔGmt) and to wild-type (ΔGwt) form of MEF2B is a very difficult computational problem, the difference in their affinity ΔΔG=ΔGwt-ΔGmt can be efficiently evaluated via thermodynamic cycle by BAR method, which is believed to be an accurate approach (see BAR method references).Unfortunately the system appears to be too complex for sampling, as it follows from the huge differences between estimations of direct and inverse mutations of D83 into the V83 residue. It is likely a consequence of the mutation conferring significant conformational shifts in the protein (in the models the DNA molecules are restrained to facilitate the sampling). Due to co-occurrence of the conformation shift, the distribution of molecular states during the binding has to be sampled multiple times under different values of the parameter of conformation shift. So the volume of configuration space to scan is a product of two reaction coordinates one of which is associated with molecular binding (distance between MEF2B and DNA) and another is associated with conformation shift itself (the distance between folded core of N-terminal domain MEF2B and its inter-domain fragment). These two forms a grid of parameters to sample the distribution of configurations on it, which is overwhelming as for today computation facilities. Despite using a long productive run for data collection—120 ns of productive MD perλ and 5.04 μs of sampling totally—even one of the most accurate procedures failed to yield a convergence of estimations and was not able to sample the problem well. Considering absence of even trends of convergence in the BAR estimations with increasing of simulation time, authors conclude that there is probably no procedure today which is sophisticated enough to be able to solve such complex system.Nevertheless both transformations (Figure6) suggest a significant increase in the stability of DNA complexes with a mutant form of MEF2B, especially forV→D. Because the method does still sample significant zones of the configuration space and as two independent samplings resulted in the same qualitative outcome, it is likely that the D83V mutant binds DNA more strongly than the wild-type protein. However, the increased affinity is achieved in a sequence independent manner by the removing of the coulomb repulsion between D83 and polyanionic DNA’s skeleton, so the overall effect of the D83V mutation is an increase in affinity to DNA coupled with a loss of sequence-specific recognition by the mutant dimers.Figure 6
The thermodynamic cycle for evaluation of ΔΔG (KJ·mol-1) of alchemical transmutation D83↔V83.


## 3.1. D83V, D83G, and D83A Mutations Are All Predicted to Destabilize the Folding of S3 Loops

The stability of the s3/1 and s3/2 loops was measured as described in the methods section—by averaging of the PMFs calculated for a set of different unfolding trajectories. The four different melting trajectories were generated through the pulling of the s3/1 or s3/2 fragment of MEF2B away from the joined a2 fragments of both monomers at temperatures 300 K, 305 K, 310 K, and 315 K. Such an approach provided a set of different energetic paths of unfolding as described by the force at reaction coordinate plots (Figure2). The behavior of the loops is however similar in all cases: because of the applied external force in each trajectory, the s3 loop melts first followed by stretching of the a3 helix (an example of the process is shown in Figure 3) which is explaining why most paths have two peaks and a valley between them on the force at reaction coordinate plots.Force strength along the reaction coordinates under different temperatures: 300 K (blue), 305 K (cyan), 310 K (green), and 315 K (red).
(a)
Unfolding s3/1 loop
(b)
Unfolding of s3/2 loopTrajectory of s3 loops melting. (a) Melting of s3/1 loop of V model at 310 K and reaction coordinateξ = 2.1, 2.6, 3.4, and 4.5 nm. (b) Melting of s3/2 loop of N-model at 310 K and reaction coordinate ξ = 2.1, 2.4, 3.4, and 4.5 nm.
(a)
(b)In order to study the stability of the s3 loops further in the different models, each of the generated unfolding trajectories was than sampled with the umbrella sampling method again at four slightly different temperatures 305 K, 310 K, 315 K, and 320 K, which resulted in 16 different individual PMFs per model per loop (see Figure S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2013/313419). The upper limit of PMF integration was set to the “valley” on the force at reaction coordinate diagram, the region ξ = 3.3 nm of artificial separation—at the point where s3/1 loops are completely melted but the a3/1 helixes remain almost intact in all models (the lower limit was set at ξ = 2.1 nm—the initial distance from the source pdb file).Although each individual trajectory shows a strong variation in the value of PMF, their average (Figure4) is appropriate for ranking and consistently different between models. Because the melting of the s3 loop and the stretching of the a3 helix are two events that cannot be clearly separated during the simulations, the PMF is still growing along the reaction coordinate at the end of the integration. To prove that the estimations are still valid, a relative PMF growing diagram for the models was built (Figure 5): even if the exact upper limit of the integration would be shifted, the difference between averaged PMFs would remain constant along the “valley” region (3.0–3.4 nm of artificial COM separation).PMFs along the reaction coordinate of s3 loops unfolding for D model (blue), C model (cyan), N model (green), V model (red), A model (orange), and G model (yellow).
(a)
PFM of s3/1 loops melting
(b)
PMF of s3/2 loops meltingAveraged PMFs of s3 loops melting for D model (blue), C model (cyan), N model (green), V model (red), A model (orange), and G model (yellow) and their difference from N model.
(a)
(b)
(c)
(d)The results of these experiments show that the activating of HDAC protein binding to theMEF2B2-DNA complex (model D) dramatically (by more than 20 KJ mol-1) increases stability of the active s3/2 fold, an energy change which is not observed with the binding of the Cabin1 suppressor (model C). In accordance with our hypothesis, the dynamic folding MEF2B protein is destabilized by the D83 mutations. The s3/1 fold is significantly destabilized by approximately 10 KJ mol-1 in the D83V mutant (model V). It is also interesting that the mutants exhibit different trends: the rare D83G mutation (model G) shows a similar yet weaker pattern of destabilization as D83V, whereas, the D83A change surprisingly has converse effects—stabilizing the active folding of s3/1 but destabilizing that of s3/2 by approximately 10 KJ mol-1.


## 3.2. The D83V Mutation Is Predicted to Improve Nonspecific Binding of MEF2B to DNA

Because the destabilization of MEF2B folding is actually a loss-of-function event, a change in the spectrum of DNA sequences recognized could explain apparent neofunctional nature of the mutations (otherwise its high recurrence in tumours remains unexplained). Therefore, the suggested mechanism is not the mitigation of the protein affinity to DNA (which would be trivial loss-of-function event) but to strengthen nonsequence-specific affinity to polyanionic frame of DNA. To test this, an alchemical mutation of MEF2 in complex with a DNA double helix was performed and the affinity difference between mutant and wild-type MEF2B binding was evaluated. Though the accurate direct estimation of free energy of DNA binding to mutant (ΔGmt) and to wild-type (ΔGwt) form of MEF2B is a very difficult computational problem, the difference in their affinity ΔΔG=ΔGwt-ΔGmt can be efficiently evaluated via thermodynamic cycle by BAR method, which is believed to be an accurate approach (see BAR method references).Unfortunately the system appears to be too complex for sampling, as it follows from the huge differences between estimations of direct and inverse mutations of D83 into the V83 residue. It is likely a consequence of the mutation conferring significant conformational shifts in the protein (in the models the DNA molecules are restrained to facilitate the sampling). Due to co-occurrence of the conformation shift, the distribution of molecular states during the binding has to be sampled multiple times under different values of the parameter of conformation shift. So the volume of configuration space to scan is a product of two reaction coordinates one of which is associated with molecular binding (distance between MEF2B and DNA) and another is associated with conformation shift itself (the distance between folded core of N-terminal domain MEF2B and its inter-domain fragment). These two forms a grid of parameters to sample the distribution of configurations on it, which is overwhelming as for today computation facilities. Despite using a long productive run for data collection—120 ns of productive MD perλ and 5.04 μs of sampling totally—even one of the most accurate procedures failed to yield a convergence of estimations and was not able to sample the problem well. Considering absence of even trends of convergence in the BAR estimations with increasing of simulation time, authors conclude that there is probably no procedure today which is sophisticated enough to be able to solve such complex system.Nevertheless both transformations (Figure6) suggest a significant increase in the stability of DNA complexes with a mutant form of MEF2B, especially forV→D. Because the method does still sample significant zones of the configuration space and as two independent samplings resulted in the same qualitative outcome, it is likely that the D83V mutant binds DNA more strongly than the wild-type protein. However, the increased affinity is achieved in a sequence independent manner by the removing of the coulomb repulsion between D83 and polyanionic DNA’s skeleton, so the overall effect of the D83V mutation is an increase in affinity to DNA coupled with a loss of sequence-specific recognition by the mutant dimers.Figure 6
The thermodynamic cycle for evaluation of ΔΔG (KJ·mol-1) of alchemical transmutation D83↔V83.


## 4. Discussion

Our hypothesis on the conformational regulation of MEF2B explains the mechanism of the experimentally observed communication between the C- and N-terminal domains. These interactions are necessary for transcriptional activation, as shown through mutagenesis studies [8], but the actual mechanism has not been previously discussed. We propose here that the folding of the s3 loops freezes the a3 helix in the major grooves of the DNA and thus permanently brings both C-terminal domains of the dimer into close proximity with each other, increasing the likelihood of the interaction with other factors, which are required for transcriptional activation [10, 26]. This mechanism can be seen as the loss of transitional entropy of the two noninteracting C-terminal domains due to the folding event in N-terminal domains. Although C-terminal domains were not modeled here, the efficiency of their restraining can be roughly estimated numerically from the known structures of MEF2 complexes. The C-terminal domains are located at approximately 3.2 nm away from the N-terminal domains dimer in the case of folded s3 loops and at least 4.3 nm away in case of unfolded s3 loops. Suggesting a homogeneous distribution of the C-terminal domains in a segment of a sphere away from the N-terminal domains dimer, the effective concentration C of the C-terminal domains is the inverse to the volume V of the segment which is proportional to the radius r of the sphere: C=1/V~1/r3, where r≥ 3.2 nm if the s3 loops of MEF2B are folded and r ≥ 4.3 nm otherwise. Considering a transcription factor X, which binds to two noninteracting C-terminal domains simultaneously, the equilibrium constant K of triple complex formation depends on the second power of effective concentration of C-terminal domains of MEF2B as K=[C]2[X]/[C2X]. Thus the probability of the triple complex formation is at least 6 times as large in the case of folded s3 loops of the N-terminal domain due to the transitional entropy loss.The proposed mechanism of the entropic interdomain communications naturally resolves apparent contradictions which were previously observed in the experimental studies of the MEF2B spatial structure. The equilibrium between the folded and disordered states of s3 loops may be affected by experimental conditions which results in the different outcomes of the experimental studies. Indeed, there are studies where the s3 loops of MEF2B are shown either as disordered [10] or with a specific structure [30–32].Another argument in favor of the proposed mechanism is the location of the recurrent D83V (A,G) mutations in the s3 loops of MEF2B. If the speculated dynamic folding is real, the recurrent mutations should affect it, and this is clearly seen in pulling MD experiments. Particularly all three mutants (D83G, D83A, and especially the most frequent D83V) show significant destabilization of loop folding (by 5–10KJ·mol-1) in contrast to the activating HDAC peptide which significantly reinforces the folding (by approximately 20 KJ·mol-1). It is important to note that the mutation is observed heterozygously in all of 27 studied cases [19]. From our modeling, a heterodimer with a mutated s3/1 loop would display a perturbation of function similar to the double mutant MEF2B dimer, allowing the mutant to act as a dominant negative in its ability to reduce the normal function of the remaining wild-type copy.However, the abolishment of the interdomain communication alone demonstrates an effective loss-of-function nature of V83 and cannot explain the higher frequency compared to other loss-of-function mutations of MEF2B encountered in lymphoma [19]. The explanation comes as the altered recognition of DNA sequence, granted by the mutation in MEF2B. The results of our thermodynamic MD calculations of the incompletely folded MEF2B dimer allow us to speculate that the mutation improves the affinity of the unfolded MEF2B dimer to DNA (by at least 24 KJ·mol-1). However, the affinity is gained in a sequence independent manner. Taken together, from the results of the loop melting experiments we conclude that D83V mutant is less efficient in activating its transcriptional program by binding DNA with increased affinity but with decreased specificity towards its normal binding sites. In other words, the mutation expands the suppression function of Cabin1 to the wide range of new genes via delocalization of MEF2B across a genome, where it may work synergistically with known repressive cancer drivers in lymphoma, such as loss-of-function insertion/deletion mutations in the activatory epigenetic mark writer MLL2 or gain-of-function mutations in repressive epigenetic marks writer EZH2 [19].Our analysis suggests that the model of MEF2B-dependent apoptosis has to be revised in the case of D83V mutated proteins. In this case, MEF2B mutant (hetero-) dimers are not located specifically at the proapoptoticNur77 gene but instead also occupy various AT-rich sequences (and thus promoters) which is typical for MADS-box proteins [15, 61]. This recruits the transcriptional suppression mediated by Cabin1 to these atypical genes. In response to B-cell receptor signaling, Cabin1 is removed via a Ca2+-Calmodulin complex from all these promoters activating transcription across a wide range of unspecific (perhaps even prosurvival) AT-rich patterns instead of a strong transcription of Nur77 gene. As a result, the apoptotic signal is weakened and the ensuing transcriptional noise rises possibly also increasing the presence of prosurvival transcripts.


## 5. Conclusions

Here we propose a new mechanism of conformational regulation for the MEF2B transcriptional enhancer. The hypothesis suggests a dynamic (un)folding of s3 loops of MEF2B in response to recognition of specific DNA patterns by N-terminal domain. The mechanism was theoretically investigated with the aid of MD pulling simulations. The proposed mechanism resolves the contradictions in previous experimental results about the structure of MEF2B, particularly the folding of the s3 fragment, which, depending on the experimental conditions, was reported either as folded or as unstructured.According to our modeling, observed mutations can interfere with the dynamic folding and contribute to the misrecognition of theNur77 promoter via the reinforcement of sequence-independent affinity of MEF2B to DNA and explain how these mutations can act dominantly. It possibly couples Cabin1 suppression to a variety of new genes assisting other known epigenetic repressors acting aberrantly in lymphomas. Minimally these mutations can titrate MEF2B mediated transcription events away from their proapoptotic function across the genome and raise transcriptional noise to interfere with apoptotic signals.Our result is relevant for the study of diffuse large B-cell lymphoma and follicular lymphoma where the MEF2B D83 amino acid residue is mutated in approximately 5% of cases [19].



---

*Source: 313419-2013-12-24.xml*