Probing the Structural Dynamics of the Unbound MAX Protein: Insights from Well-Tempered Metadynamics

TL;DR

This study reveals a stable structure of the MAX protein, which could help in developing drugs to target cancer-related MYC/MAX interactions.

Contribution

The study identifies a stable conformation of unbound MAX protein using metadynamics, revealing potential druggable sites.

Findings

MAX protein can adopt a stable, thermodynamically favorable conformation.

Metastable structural motifs in MAX may serve as potential drug targets.

The findings provide a structural framework for future drug design targeting the MYC/MAX axis.

Abstract

One of most relevant
challenges in tumorigenesis is the
association
of MYC and MAX proteins, whose related cancers remain undrugged. In
particular, the disordered regions shown by those oncogenes make their
structural characterization and the development of new drugs a truly
hard task. To address these challenges, we employed the enhanced-sampling
well-tempered metadynamics method to systematically explore the conformational
space of the unbound MAX protein. Our results revealed, for the first
time, a well-defined and thermodynamically favorable conformation
of monomeric MAX. This is a remarkable finding, as it demonstrates
that regions of MAX previously considered persistently disordered
are capable of adopting stable, folded structures under specific conditions.
Moreover, our findings also suggest that the metastable structural
motifs observed in this work may harbor druggable sites, particularly
relevant for strategies aiming to target MAX directly or to disrupt
its interaction with MYC, thereby modulating oncogenic signaling pathways.
The present study establishes a new structural framework for understanding
the dynamics of MAX and provides a foundation for future structure-based
drug design targeting the MYC/MAX axis. Finally, our work offers a
strategic blueprint for investigating similarly challenging drug targets.

Figures

(A) Crystal structure of a MYC-MAX dimer recognizing DNA obtained in ref. (B) RMSD profiles over the full simulation time of triplicate MD simulations initiated with different velocity seeds, spanning a total of 4 μs. The frames corresponding to representative structures are highlighted. (C) Cluster density profile derived from PCA clustering. (D) Seven distinct clusters identified from PCA results, with the representative structure of each cluster indicated by a red cross. (E) The transition pathways of MAX are illustrated as follows: from cluster 2 to cluster 1 in the first MD run (rep1, shown in subfigure B), from cluster 0 through cluster 6 to cluster 4 in the second independent MD run (rep2), and from cluster 5 to cluster 3 in the third MD run (rep3). The key structural transition is highlighted by the progressive decrease in the angle between the two helices.

(A) Representative structures from the seven clusters aligned based on residues 222-278 of Helix 2 of MAX, clearly illustrating the progression of conformational transitions. Colors correspond to those shown in Figure E. (B) Illustration of two collective variables (d1 and d2), Cα atoms of Ile218 and Ile250 are shown in van der Waals (VdW) representation.

Two-dimensional free energy surface F(d1, d2) (kJ/mol) of unbound MAX in solution obtained from WT-MetaD simulations. Regions with free energy exceeding 90 kJ/mol are shown in red. “c0” to “c6” denote the coordinates of the representative structures from clusters 0 to 6 projected in this plot (see Figure ). Clusters c1 and c2 (from the rep1 MD simulation) are shown in black, c0, c4, and c6 (from rep2 MD simulation) are shown in brown, and c3 and c5 (from rep3 MD simulation) are shown in purple.

Partial free energy surface F(d1, d2) (kJ/mol) is presented here, highlighting newly discovered stable states within basins “A”-“L”. The corresponding MFEP connecting selected basins are also depicted. Representative MAX conformers from each basin are shown, highlighting the helical domain (I264-A280) near the C-terminus in all cases. Structures of MAX for basin C and I are not shown due to their short lifetime.

1D integrated free energy profiles for two CVs at 3500 ns. Basin regions marked with same labels (A-L) as in Figure . In order to directly compute the height of free energy barriers, absolute minima of each free energy are set equal to zero.

(A) MFEPs between A to CS (red), A to G (purple), and G to CS (orange) on the free energy surface F(d1, d2) (kJ/mol). (B) Aligned structures of state G and the CS. Angles between vectors along the helical domain are shown here. (C) Electrostatic potential of MAX in states A, D, and G in units of k B T/e (k B T/e ≈ 26.7 mV). The potential drug-binding pocket of each state is marked with a pink star, and positive side chains of these pockets are highlighted here. (D) Charge distribution of small-molecule c-MYC inhibitor 10058-F4.

Most favorable structures of MAX-inhibitor complex predicted by HADDOCK2.4. The inhibitor 10058-F4 is displayed in van der Waals representation in all cases. Figures were generated using VMD.

Tables

1. WT-MetaD Simulation Parameters Used in This Work

parameter	value
Gaussian width of d1 (nm)	0.2
Gaussian width of d2 (nm)	0.2
Gaussian hill (kJ/mol)	1.2
deposition stride (ps)	1
biasfactor	8
simulation time (ns)	3500

2. Coordinates of the Representative Basins (Minima) Corresponding to the Figure with the Helical Domain in Each Structure Highlighted Here

basin	d1 (nm)	d2 (nm)	free energy (kJ/mol)	helical domain
A	0.45	1.61	0.00	E248-R255 and I264-Q276
B	0.45	2.18	4.77	E248-R255 and I264-A280
C	0.45	1.05	6.68
D	1.09	1.97	11.37	A246-K256 and I264-R279
E	2.08	1.26	17.30	I242-H260 and I264-R279
F	2.63	2.87	20.77	K245-K256 and I264-A280
G	3.21	0.77	25.69	K203-L225, A240-H260 and I264-L281
H	5.05	1.96	30.55	K245-R255 and I264-L281
I	4.62	0.69	32.23
J	6.06	0.65	33.23	T247-R255 and I264-G282
K	3.69	4.22	38.40	A246-K256 and I264-R279
L	0.52	4.72	41.21	K236-R239, L243-H260 and I264-R279
CS	8.98	1.96	74.91	K203-S228 and R239-R279

3. Free Energy Barriers (ΔF) between Basins (kJ/mol) along the Corresponding MFEP and Their Transition Half-Lives (ps), as Predicted by the Nudged Elastic Band Method and Eyring-Polanyi Equation, Respectively (See Section “Methods”)

state transition	ΔF (kJ/mol)	transition half-life (ps)
A to B	12.30	15.39
B to A	6.86	1.73
A to C	12.06	13.94
C to A	2.65	0.32
A to D	17.81	139.71
D to A	6.83	1.71
B to D	11.15	5.79
D to B	3.93	0.54
D to E	9.86	5.78
E to D	3.93	0.53
E to F	7.03	1.86
F to E	3.56	0.46
E to G	15.04	46.08
G to E	6.65	1.59
D to F	13.07	20.89
F to D	3.66	0.48
F to H	14.26	33.74
H to F	4.48	0.67
G to I	8.79	3.76
I to G	2.26	0.27

Full text

Introduction

The V-myc avian myelocytomatosis viral oncogene homologue (MYC) and the MYC-Associated factor X (MAX) proteins are key transcription factors involved in regulating gene expression related to cell growth, proliferation, metabolism and apoptosis. MYC functions as a potent oncogene and requires dimerization with MAX to bind DNA at E-box sequences and activate transcription of target genes. MAX, in contrast, also forms dimers with other proteins like the MAX dimerization protein 1 (MXD1), which act as transcriptional repressors, thereby balancing MYC’s activity. The MYC-MAX network’s central role in tumorigenesis has become a significant target for anticancer therapies, although direct inhibition of MYC and/or MAX remains challenging due to its intrinsically disordered structure and widespread cellular functions.? Whereas crystal structures of MAX-MYC dimers have been determined with precision, both including ?,? and excluding DNA,? the determination of suitable drugging pockets is still quite elusive.

MAX encodes a ubiquitously expressed transcription factor that serves as an essential dimerization partner for MYC family proteins, facilitating their ability to regulate gene expression. ?,? MAX plays a pivotal role in modulating cellular processes such as proliferation, differentiation and apoptosis by forming heterodimers with MYC or other members of the MYC-MAX-MXD network. ?−? ? These interactions enable the binding of the resulting complexes to E-box sequences in DNA, where they activate or repress target genes depending on the dimeric partner involved. ?,? Similarly to MYC, MAX possesses a basic helix-loop-helix leucine zipper (bHLHLZ) domain at its C-terminus, which is crucial for both dimerization and DNA binding. ?,? However, in contrast to MYC, MAX lacks an intrinsic transactivation domain, acting instead as a scaffold that confers transcriptional activity or repression through its binding partners. ?,? The evolutionary conservation of MAX underscores its central role in balancing oncogenic and tumor-suppressive signals within the cell. ?,?

In many cancers, dysregulation of the MYC/MAX/MXD transcriptional network leads to an imbalance favoring MYC-MAX heterodimers, which drive uncontrolled cell proliferation and repress differentiation. Although MAX itself is not typically overexpressed or mutated, its obligatory role as a dimerization partner makes it a critical enabler of MYC-driven oncogenesis. ?,? MAX is expressed ubiquitously and forms stable bHLHLZ-mediated dimers with MYC, enabling binding to E-box DNA sequences and transcriptional activation of target genes involved in cell growth, metabolism and survival. ?,?

Given the structural challenges in directly targeting MYC-an intrinsically disordered protein lacking well-defined binding pockets-MAX has emerged as an attractive alternative target for therapeutic intervention. ?,? Disrupting the MYC-MAX interaction could effectively suppress MYC oncogenic activity without requiring direct inhibition of MYC itself. Several strategies have been proposed to this end, including the use of dominant-negative MYC mutants such as OmoMYC? (OMO-103, currently in clinical trials), which sequesters MAX and blocks MYC-dependent transcription.? Several small molecules, such as 10058-F4,? KJ-Pyr-9,? and MYCMI-6,? have been developed to interfere with the MYC-MAX interaction directly. Also, peptides that mimic MAX or stabilize alternative dimerization partners (such as MXD proteins) could shift the transcriptional balance toward growth suppression. ?,?

In order to identify viable drug-binding regions on MAX, it is essential to characterize its structural configurations. In this study, we leverage all-atom molecular dynamics (MD) simulations with high temporal and spatial resolution to model the conformational behavior of MAX and its potential interactions with MYC. We have chosen MD simulations due to its suitability to explore and characterize the conformational space of highly dynamic proteins. ?,? By analyzing the dynamic structural ensemble of MAX, we aim to uncover transient or stable pockets that could serve as potential drug-binding sites. Subsequently, we used the Well-Tempered Metadynamics (WT-MetaD)? simulations to map the free energy landscape of MAX in solution, which will allow us to identify low-energy conformations with sufficient structural rigidity to accommodate small molecules, paving the way for the development of novel inhibitors targeting the MYC-MAX axis in cancer therapy.

Methods

System Preparation

and MD Simulations

The amino acid sequence of the MAX protein is available in UniProt (ID: P61244). There are several tridimensional structures of the MAX protein available in the Protein Data Bank, always in the form of a complex including another protein (MYC). We have selected the structure with PDB code 1NKP obtained in a detailed study of MAX recognition functions,? but we should note that this MAX configuration is not different from that reported in other crystal structures. Hence, the MAX protein coordinates extracted from that PDB entry were employed as starting point for our simulations. Our simulation system consists of the MAX protein solvated in an explicit water box containing 0.15 M KCl to mimic the cell ionic strength. The system contains a total of 159,804 atoms. The system coordinates were generated using the CHARMM-GUI solution builder, ?−? ? and the CHARMM36m force field? was adopted. All-atom MD simulations were conducted using AMBER24 ?,? software. The system was minimized for 10000 steps followed by three 125 ps equilibrium runs while gradually reducing the harmonic constraints (NVT ensemble). Production runs were performed with the NPT ensemble, whereas Langevin dynamics with the friction coefficient 1.0 ps^–1^ was used for temperature regulation (310 K). The Monte Carlo barostat was used for the pressure regulation (1 bar). All bonds involving hydrogens were set to fixed lengths, allowing fluctuations of bond distances and angles for the remaining atoms. The time step was set to be 2 fs, and the frames were saved every 100 ps for analysis. The particle mesh Ewald method was automatically used by AMBER24 to calculate the electrostatic interactions, and the van der Waals interactions were calculated using a cutoff of 1.2 nm. Periodic boundary conditions were taken in three directions of space.

Well-Tempered Metadynamics

The enhanced sampling method WT-MetaD applies a time-dependent biasing potential along a set of collective variables (CV) by adding a Gaussian additional bias potential to the total potential in order to overcome barriers larger than k B T, with k B being Boltzmann’s constant and T the temperature. In the method, sampling is performed on a selected number of degrees of freedom (s _ i _), i.e., a chosen set of CV. For each degree of freedom, the biased potential V(s _ i _,t) is a dynamical function constructed as the sum of Gaussian functions:

$$[eqn]$$

where k is an integer, τ is the Gaussian deposition stride, W(kτ) is the height of the Gaussian and σ_ i _ is the width of the Gaussian for the ith CV. Full details are provided in refs ?,?,? . In WT-MetaD, the biased potential can force the system to move around all possible states inside a particular range of the subspace of the selected CV. In the present work, we have taken a specific approach similar to previous works ?−? ? ? where the CV were defined from distances between groups of particular amino acids. After analyzing the protein’s dynamics from the classical MD runs, we have selected two collective variables: d1 and d2, where d1 refers to the distance of center of mass of group1 (backbone atoms of residues 202-204) and group2 (backbone atoms of residues 282-284) and d2 refers to the distance between C_α_ atoms of Ile218 and Ile250 to explore the conformational space of MAX, see Section “Results and Discussion” for a visualization of these CVs over the protein structures. The WT-MetaD simulations were performed using the PLUMED plugin (version 2.9.2) implemented in Gromacs MD package (version 2023.5). ?−? ? As initial configuration for the WT-MetaD simulations we considered the same relaxed crystal structure (CS) employed in the unbiased MD simulations. The main parameters of the WT-MetaD simulation are reported in Table. The V-rescale thermostat was considered for accurate temperature control and proper thermodynamic sampling at 310.15 K. Pressure was controlled to be 1 bar by a Parrinello-Rahman piston with a damping coefficient of 5 ps^–1^. Periodic boundary conditions in the three directions of space were also considered. More details, including analysis of the convergence of the metadynamics simulations (Figures S3–S7) can be found in the “Supporting Information” (SI).

1: WT-MetaD Simulation Parameters Used in This Work

Data Analysis and Visualization

The backbone atoms (N, C_α_, C, and O) of MAX were used to calculate the root mean squared deviation (RMSD) profile through RMSD Trajectory Tool of Visual Molecular Dynamics (VMD) version 1.9.3.? To explore the conformational landscape of the MAX protein during the molecular dynamics simulations, we carried out cluster analysis using the atomic coordinates of the backbone atoms, applying Principal Component Analysis (PCA)­s.? The analysis was performed using scikit-learn? and MDAnalysis Python libraries.? Pairwise distance matrices were constructed by extracting the positions of the backbone atoms throughout the entire trajectory using MDAnalysis. These structural descriptors were then employed as input data for PCA analysis using scikit-learn. Subsequently, the KMeans clustering algorithm? was applied to the PCA-transformed data to partition the conformations into distinct clusters. Each cluster corresponds to a group of structurally similar protein conformations, potentially representing distinct metastable states sampled during the simulation. Then the representative structures from each cluster which are closest to the cluster centroid were identified for further structural characterization and visualization. This clustering approach provides a simplified yet informative view of the protein’s conformational ensemble, facilitating the interpretation of structural transitions and dynamic behavior. For each representative cluster structure, we report the angle between the two helices, calculated considering as the direction of each helix the one defined by the α carbon atoms of residues 210 and 221 of helix 1 and residues 249 and 267 of helix 2.

The R-package metadynminer ?,? was employed to analyze the WT-MetaD outputs. This package provides a set of tools for the efficient extraction, visualization, and interpretation of free energy surfaces (FES) computed from WT-MetaD trajectories. In this study, metadynminer was specifically employed to reconstruct the FES using the fast Bias Sum algorithm. Convergence of the free energy calculations was analyzed in detail in the SI. The neb function of the same package was employed to compute the minimum free energy paths (MFEP) by using the Nudged Elastic Band method.? MFEP describe the most probable transition routes between these states on the FES. The free energy profiles, free energy barriers and the transition states identified along the transitions were reported in detail in the main text and SI. The neb function of metadynminer also estimates the half-life of these states? using the classical Eyring-Polanyi equation of chemical kinetics,? with a transmission coefficient assumed to be 1.

In addition to the full free energy map and the free energy profile along the MFEP, we have also computed the 1D free energy profiles F(d 1) and F(d 2) (where d 1 and d 2 are the CVs) by thermally averaging all contributions from d 2 and d 1 respectively in the full free energy profile F(d 1, d 2). More explicitly, F(d 1) is calculated using ?,?

$$[eqn]$$

where β = 1/(k B T), k B is the Boltzmann constant and T is the absolute temperature. F(d 2) is obtained by simply exchanging the roles between d 1 and d 2 in eq. The analysis tool VMD? and Chimera? were used for visualization and analysis purposes.

Results and Discussion

Structural Transition Pattern of MAX

The 3D structure of the MAX protein is shown in FigureA, as obtained experimentally in a MYC-MAX dimer recognizing a common DNA target.? The structure shows how the two helices of MAX are used to interact simultaneously with MYC and the DNA sequence in a way that has been intuitively compared with a specialized hand, perfectly designed to perform two tasks: shake hands with other proteins and grip onto DNA. Given this structure, a natural question is whether the MAX protein adopts a stable conformation when not dimerized with MYC, and more broadly, what conformational states MAX can assume based on different arrangements of its two helices. This question should be relevant both for understanding how MAX performs its biological function but also should be relevant in the design of drugs interacting with MAX. To this end, we performed first exploratory MD simulations (described in this subsection) which were employed to identify relevant collective variables to be employed in subsequent extensive free energy calculations performed in the remaining subsections.

(A) Crystal structure of a MYC-MAX dimer recognizing DNA obtained in ref. (B) RMSD profiles over the full simulation time of triplicate MD simulations initiated with different velocity seeds, spanning a total of 4 μs. The frames corresponding to representative structures are highlighted. (C) Cluster density profile derived from PCA clustering. (D) Seven distinct clusters identified from PCA results, with the representative structure of each cluster indicated by a red cross. (E) The transition pathways of MAX are illustrated as follows: from cluster 2 to cluster 1 in the first MD run (rep1, shown in subfigure B), from cluster 0 through cluster 6 to cluster 4 in the second independent MD run (rep2), and from cluster 5 to cluster 3 in the third MD run (rep3). The key structural transition is highlighted by the progressive decrease in the angle between the two helices.

We have performed triplicate MD simulations of the MAX protein solvated in salty water (150 mM of KCl mimicking physiological salt concentration), each initiated with distinct velocity seeds and spanning a total of 4 μs. All simulations used the crystal structure of the MAX protein, extracted from the MYC-MAX-DNA complex (FigureA) as the starting structure (see Section “Methods” for full details). During the MD trajectories the time evolution of quantities such as the RMSD (FigureB) indicate the possibility of substantial changes in the protein conformation, even at the time scales explored by the unbiased MD simulations. In order to identify the essential features of the protein configurations explored during the MD trajectories, we have performed PCA analysis. We found that a 64.1% of the variance of the data can be explained by the first principal component. Adding a second component, the PC1 and PC2 variables explain 74.4% of the total variance (PC1 = 64.1%, PC2 = 10.3%). Adding more components does not significantly improve the variance explained, unless one adds a large number of PCA components (see Figure S1 of SI). Therefore, we have decided to restrict our analysis to two PCA components. Once the two PCA components were obtained, we projected the trajectory into the PCA components and clustered the data (see Section “Methods”) and obtained representative structures for each cluster (FigureC-E).

We evaluated clustering with varying numbers of clusters in order to identify the optimal choice. The results obtained with two different number of clusters (7 clusters and 8 clusters) are reported in Figures and S2, respectively. All these clusters correspond to a significant proportion of the trajectory configurations (FigureC) and the structures characterizing these clusters are clearly different, as shown by the values of the PCA components (FigureE), their different RMSD (FigureB) and their different angle between the two helices that define the MAX protein structure (FigureD). We found that 7 clusters were enough to represent the variability of the protein structures found in the MD simulation. Considering a larger number of clusters does not provide significantly different protein conformations, as illustrated with the results obtained for 8 clusters shown in the SI (Figure S2). The results from the 8-cluster analysis were largely consistent with those obtained using 7 clusters, with the additional cluster further subdividing the existing groups without revealing fundamentally new structures. This consistency indicates that clustering into 7 groups captures the essential data structure and is sufficient for interpretation (Figure).

FigureC shows that cluster2 exhibits the highest population density, followed by cluster3, cluster5, cluster0, cluster1, cluster6, and cluster4 (note that the labeling of the clusters is arbitrary, not based on any physical quantity or time sequence). The time evolution between these clusters during the simulation was identified by tracking the simulation frames corresponding to each representative structure along the unbiased MD simulation timeline. Representative structures from the seven distinct clusters identified through PCA, along with the corresponding transition pathway, are shown in FigureE. Across all three MD replicas, we observed a consistent structural transition in which the angle between the two helices in MAX progressively decreases along the trajectory. This reduction in interhelical angle reflects a gradual compaction of the MAX domain, suggesting a shift toward a more closed conformational state. These results suggest that the observed transition is a reproducible and intrinsic feature of MAX’s conformational landscape and not an artifact of a particular simulation run.

Selection of Collective

Variables in WT-MetaD Simulations

After performing unbiased MD simulations, we will perform a more detailed analysis of the possible conformation space of MAX by employing WT-MetaD simulations, as described in the Section “Methods”. Before performing the simulations, the WT-MetaD technique requires the selection of the appropriate CV employed in the exploration of the conformation space. As it was mentioned above, we will consider as CV distances between particular amino acids, as it was assumed in several previous works. ?−? ? ? The selection of these distances to be used as CV is made by analyzing the transitions between the seven different clusters obtained in the unbiased MD simulations (FigureE). The geometrical changes undergone by the protein during these transitions can be clearly seen by aligning the representative structures of these clusters, as shown in FigureA.

(A) Representative structures from the seven clusters aligned based on residues 222-278 of Helix 2 of MAX, clearly illustrating the progression of conformational transitions. Colors correspond to those shown in Figure E. (B) Illustration of two collective variables (d1 and d2), Cα atoms of Ile218 and Ile250 are shown in van der Waals (VdW) representation.

The conformational transitions are characterized by the rotation of helix 1 (residues 202-228) toward helix 2 (residues 239-281). This rotation ultimately leads to a significantly reduced angle between helix 1 and helix 2, as illustrated in FigureA and a reduced distance between the N- and C-termini (indicated in FigureB as d1), along with a marked decrease in the distance between the C_α_ atoms of residues Ile218 and Ile250 (indicated in FigureB as d2). As a result of this conformational rearrangement, the overall distance between the N- and C-termini of the 202-284 residue segment in unbound MAX is substantially reduced, leading to the formation of a more folded structure, represented by cluster4.

These observations were derived from the initial set of unbiased MD simulations conducted in the first phase of this study. Based on these structural changes, we selected d1 and d2 as the collective variables for subsequent enhanced sampling via WT-MetaD simulations.

Free Energy Landscape of MAX

After selecting d1 and d2 as CV, we performed WT-MetaD simulations. The results are summarized in Figure. We ensured that the free energy calculation is properly converged. As shown in Figure, the free energy landscape of MAX reveals pronounced structural diversity within the low free energy region, populated by several metastable states. This indicates a rich ensemble of conformational states readily accessible in solution. The different states were found during the previous unbiased MD simulations. The results demonstrate that over the 4 μs of an unbiased MD simulation, the system evolves from a high energy, crystal-like conformation of MAX (cluster5) toward a significantly lower-energy structure (cluster4). The free energy difference between the most stable state found in unbiased MD (cluster4) and the crystal-like state is ∼60 kJ/mol. Also, this figure shows that the unbiased MD simulations were not reaching the lowest free energy conformation of the protein, in spite of the long simulation time. These results emphasize that, even for a relatively small protein motif in this study (83 amino acids), the crystal structure can differ significantly from the minimum free energy configuration. They also underscore the importance of using enhanced sampling methods, such as WT-MetaD, to explore the structural variability of proteins. This enhanced sampling strategy enables us to explore the broader conformational landscape of MAX and ultimately identify previously inaccessible, stable drug-binding pockets in underexplored regions of the protein.

Two-dimensional free energy surface F(d1, d2) (kJ/mol) of unbound MAX in solution obtained from WT-MetaD simulations. Regions with free energy exceeding 90 kJ/mol are shown in red. “c0” to “c6” denote the coordinates of the representative structures from clusters 0 to 6 projected in this plot (see Figure ). Clusters c1 and c2 (from the rep1 MD simulation) are shown in black, c0, c4, and c6 (from rep2 MD simulation) are shown in brown, and c3 and c5 (from rep3 MD simulation) are shown in purple.

Looking in more detail at the free energy landscape, it is possible to observe several primary basins (Figure). The details of these molecular structures (labeled A through L in Figure) are also provided in Table.

Partial free energy surface F(d1, d2) (kJ/mol) is presented here, highlighting newly discovered stable states within basins “A”-“L”. The corresponding MFEP connecting selected basins are also depicted. Representative MAX conformers from each basin are shown, highlighting the helical domain (I264-A280) near the C-terminus in all cases. Structures of MAX for basin C and I are not shown due to their short lifetime.

2: Coordinates of the Representative Basins (Minima) Corresponding to the Figure with the Helical Domain in Each Structure Highlighted Here

Basins A, B, and C, located in the lower-left region of the free energy surface (d1 ∼ 0.45 nm, d2 ≤ 2.2 nm), represent the most thermodynamically stable conformations, as indicated by their deep free energy wells, see Table. These basins are tightly clustered, suggesting a compact ensemble of low-energy structures, associated with the native or folded state of MAX (see the corresponding structures of different basins in Figure). Basins D, E, F, and G occupy moderately low-energy regions (11-26 kJ/mol) and may represent alternative metastable states that are structurally related to the native ensemble but differ in loop orientations or side-chain packing. In contrast, basins H-L, that are located in higher-energy regions (between 31-42 kJ/mol), refer to less stable and partially unfolded conformations. This free energy landscape reveals broader and more diffuse features in this region with greater structural flexibility. We observed a substantial free energy difference of 74.91 kJ/mol between the most stable basin (A) and the experimentally resolved crystal structure (CS). This substantial difference between the most stable configuration of MAX and its crystal structure has important implications for the development of drugs targeting MAX. This significant free energy gap underscores the difficulty of targeting MAX’s CS conformation with small molecules due to its low state population. The stable conformations identified in this study may offer promising targets for the development of effective inhibitors against MAX-associated diseases, much better than the crystal structure of MAX.

Identification of Transition Paths and Transition

States

The WT-MetaD technique reveals not only the free energy basins (the most stable configurations) but also the free energy barriers that the protein needs to overcome to shift between stable states. Furthermore, the maxima along each path give indications of the approximated location of transition states (TS) between stable states. In Table we show the free energy barriers between the different basins of interest seen in Figure. The transition half-lives between these basins, estimated using the Eyring-Polanyi equation (see Section “Methods”), are also provided. The corresponding free energy profiles along the MFEP for each case are shown in Figure S8 of SI, including the estimated locations of local TS between stable basins, and the raw data can be found in our Github repository.? The energy barriers between low-energy basins and higher-energy regions (e.g., A to D, E to G, and F to H) suggest that large-scale rearrangements are thermodynamically unfavorable and likely require crossing high activation barriers. The half-life values along the MFEP between basins are consistent with the number of conformers sampled for each basin in the WT-MetaD trajectories. A shorter half-life indicates a faster kinetic transition between states for MAX, leading to lower accumulation in individual state populations (i.e., fewer structurally relevant conformers). Notably, the structures of MAX corresponding to basin C and I, which exhibit very short transition half-lives (0.32 and 0.27 ps, respectively), were not captured in the WT-MetaD output, where frames were saved every 2 ps, indicating that they are transient states.

3: Free Energy Barriers (ΔF) between Basins (kJ/mol) along the Corresponding MFEP and Their Transition Half-Lives (ps), as Predicted by the Nudged Elastic Band Method and Eyring-Polanyi Equation, Respectively (See Section “Methods”)

In addition to the complex picture provided by the MFEP, it is possible to obtain a more simplified but useful picture by considering one-dimensional (1D) profiles where free energy depends of only one CV, after the other CV has been integrated out over their Boltzmann factor. Previous works? have shown the relevance of these 1D free energy profiles since, under appropriate conditions, they can be directly related to experimentally measured free energies ΔG. In Figure we show the results for the 1D free energy profiles F(d 1) and F(d 2) obtained by thermally averaging all contributions from d 2 and d 1 respectively (see the Section “Methods” for further details).

1D integrated free energy profiles for two CVs at 3500 ns. Basin regions marked with same labels (A-L) as in Figure . In order to directly compute the height of free energy barriers, absolute minima of each free energy are set equal to zero.

The obtained free energy profiles F(d 1) and F(d 2) are rather smooth, with less features than the MFEP free energy profiles shown in Figure S8. The only remarkable feature of these profiles is seen in the case of F(d 1) in the region connecting basins A and D. It should be noted that this is the transition which involves the largest free energy and largest transition half-life, as identified by the MFEP (see Table). The F(d 1) profile between A and D states presents a small local maxima close to D and a deep energy minimum around A. Interestingly, the profile of F(d 1) in this region (Figure) is close to the MFEP profile between these same states (Figure S8). The free energy difference involved in the transition between A and D following the MFEP path is 17.8 kJ/mol (see Table) and the prediction from of F(d 1) is about 15 kJ/mol (Figure). This comparison demonstrates that although the transition between A and D involves changes in both CV (see Figure), the free energy changes are essentially due to only one of the CV (in this case d 1). This is not the case for most of the other transitions shown in Table, as seen by comparing Figure with Figure S8, which involve more complex free energy contributions from changes in both CV and that need to be described by the MFEP connecting them.

Overall, the results for the free energy barriers found in the calculations reported here show many possible transitions between free energy basins involving free energies much larger than the thermal energy at 310 K (∼2.6 kJ/mol). They can be compared to other biological bottlenecks, such as those related to the RAS family of oncogenes, where barriers of ∼20 kJ/mol due to changes in specific distances and orientations were estimated for the activation of mutated RAS proteins by Fourier Transform Infrared Spectroscopy? and simulations? or to the case of KRAS-G12D, when bound to a cell membrane, where angular barriers in the range 8-17 kJ/mol were estimated.?

Location of Potential Targeting Pockets at

the Surface of MAX

After characterizing the possible configurations of MAX, our interest now is to identify which of these configurations are of interest for in silico drug design efforts. A possible therapeutic strategy could be to target unbound MAX, limiting its availability for MYC binding, as described in the Section “Introduction”. We should note here that the crystal structure employed as initial condition in our simulations corresponds to a MAX protein in a MYC-MAX dimer. The unbound MAX protein in solution can be expected to be in one of the free energy basins described in the previous subsections. Therefore, to identify druggable targets from MAX configurations we should explore the configurations and free energy landscape connecting the minimum free energy state identified in our calculations (state A in Table) with the CS structure.

By following the MFEP from the most stable state (basin A) to the CS, it is clear that the MFEP, depicted in red in FigureA, passes through basin D, leading to a high free energy barrier of approximately 75 kJ/mol. Interestingly, structure G shows a high degree of similarity to the CS within the helical domains when compared to other basins (highlighted in bold characters in Table). Upon aligning the structures of basin G and CS, it becomes evident that MAX only needs to rotate residues 202-225 by 164.2° and residues 264-284 by 133.1° to adopt the conformation of CS (see FigureB). This suggests that state G may represent a transition state to adopt the conformation of the CS during the process of the MAX-MAX/MAX-MYC dimerization. Given the higher population of state G, reflected by the greater number of conformers compared to other basins in Figure, we further examined the MFEP from basin A to G, which proceeds through basin D and involves a free energy barrier of 32.01 kJ/mol. Additionally, the MFEP from basin G to the CS revealed a free energy difference of 49.60 kJ/mol, which is smaller than that from basin A to the CS. These findings highlight states A as the most stable structure of MAX in solution, D and G as key intermediates, suggesting their potential as promising therapeutic targets in MAX-associated diseases. An interesting feature of these structures is the fact that they expose positively charged regions. In FigureC we show the electrostatic potential of the representative structures of the A, D and G states calculated at the Poisson-Boltzmann level of theory using the APBS? software. As seen in that figure, there are large pockets of these structures at potentials around ∼80 mV. This suggests electrostatic interaction as a suitable mechanism for binding to MAX. FigureD shows the charge distribution of the small-molecule c-MYC inhibitor 10058-F4, which has been proposed to interact with both MAX and MYC and inhibit their heterodimerization.?

(A) MFEPs between A to CS (red), A to G (purple), and G to CS (orange) on the free energy surface F(d1, d2) (kJ/mol). (B) Aligned structures of state G and the CS. Angles between vectors along the helical domain are shown here. (C) Electrostatic potential of MAX in states A, D, and G in units of k B T/e (k B T/e ≈ 26.7 mV). The potential drug-binding pocket of each state is marked with a pink star, and positive side chains of these pockets are highlighted here. (D) Charge distribution of small-molecule c-MYC inhibitor 10058-F4.

Our results suggest that the binding of this molecule to states A, D, and G of MAX should be possible, driven by electrostatic complementarity. To verify this hypothesis, we conducted additional docking analyses between MAX and the 10058-F4 inhibitor. Full details are reported in the SI. The docking outcomes, presented in Figure, show the most favorable complex structures for states A, D, and G and align well with our preliminary predictions in FigureC. Coordinates of these complexes are available in our GitHub repository.

Most favorable structures of MAX-inhibitor complex predicted by HADDOCK2.4. The inhibitor 10058-F4 is displayed in van der Waals representation in all cases. Figures were generated using VMD.

Conclusions

The importance of the MYC-MAX dimer lies in its ability to associate with DNA in the cell nucleus and, more critically, in the essential role this complex plays in the dysregulation of gene expression in roughly 70% of all known cancers. Despite extensive research, no effective strategies for directly targeting MYC have yet been developed. Because direct inhibition of MYC has proven extremely challenging, MAX has emerged as an attractive alternative target for novel therapeutic approaches. ?,? In this context, the present work reports a detailed study of the FES of MAX in solution, aiming to identify previously uncharacterized druggable pockets associated with the most relevant minima of the FES. We have observed that the monomeric, unbound MAX motif exhibits a rich conformational ensemble, as revealed by extensive classical MD (4 μs) and WT-MetaD simulations (3.5 μs). Classical MD runs revealed five dominant structural clusters, which guided the selection of two distances as CV for the subsequent WT-MetaD enhanced sampling to further explore the structural diversity of MAX.

The 2D FES constructed from WT-MetaD simulations revealed a broad conformational space, featuring numerous low-energy basins corresponding to metastable states. The most stable basin (A) lies ∼75 kJ/mol below the experimentally resolved crystal structure, indicating that the latter is a low-populated state in solution. Basins A and B form a compact cluster of highly stable conformations with tightly packed helices, while D-G represent moderately stable states. By contrast, higher-energy basins (H-L) are associated with partially unfolded or extended conformations, underscoring the intrinsic flexibility of MAX. The Minimum free energy Path analysis using the package R-Metadynminer indicated that the transition from A to D exhibits the largest energy barrier (ΔF ∼ 18 kJ/mol), suggesting a significant kinetic bottleneck. In contrast, transitions among other basins (e.g., A to B or B to D) have lower energy barriers and shorter lifetimes. Certain basins such as C and I were identified as highly transient (lifetime ≤ 0.5 ps), leading to their absence in the WT-MetaD output frames. This indicates that not all minima contribute equally to the long-lived conformational ensemble of MAX.

Importantly, the substantial energy gap between the stable solution-phase conformations and the crystal structure highlights the need to target MAX in its native, lowest-energy states. The stable conformations identified in the present study may therefore represent promising candidates for the development of effective inhibitors, extending beyond the long-studied crystal structure of the MYC-MAX complex. In particular, the identification of persistent pockets within low-energy basins (particularly A, D, and G) provides promising templates for rational structure-based drug design against the MYC-MAX interaction interface. Furthermore, electrostatic potential maps of states A, D, and G reveal positively charged pockets, while the small-molecule MYC inhibitor 10058-F4 exhibits a negatively charged electrostatic profile, suggesting potential electrostatic complementarity that could be exploited for binding. The resulting complex structures predicted by additional docking analyses between MAX and the 10058-F4 inhibitor align well with our predictions.

Supplementary Material