[Preview] Energetics of Allosteric Communication in Ubiquitin Revealed by Hybrid MCTS-Langevin Simulations

This post provides a preliminary overview of our upcoming research paper developed through collaboration with Fatma Şengüler Çiftçi and Dogukan Uraz Tuna. The full paper is currently in preparation and not yet published. Consider this a conceptual introduction to our approach and preliminary findings.

Overview

Exploring protein conformational landscapes and identifying allosteric communication pathways remain significant challenges in computational biophysics. Our study presents a novel hybrid computational approach combining Monte Carlo Tree Search (MCTS) with Langevin Dynamics (LD) simulations to enhance protein conformational sampling.

Key Innovations

We applied our MCTS-LD protocol to ubiquitin (PDB: 1UBQ), generating conformational trajectories that reveal:

• Five distinct allosteric communication pathways between MET$_{1}$ and GLY$_{76}$
• High fluctuations in the C-terminal region (normalized values 1.54-2.40)
• A critical signal integration hub at LEU$_{73}$ and ARG$_{74}$ (r = 0.78, p < 0.01)

Methodological Framework

Our hybrid approach uniquely combines:

1. MCTS exploration - Strategic navigation of the conformational landscape using the Upper Confidence Bound for Trees criterion:

   $$UCT(n) = \frac{Q(n)}{N(n)} + c\sqrt{\frac{\ln N(p(n))}{N(n)}} + \frac{|\nabla V(n)|}{\max(|\nabla V|)}$$

2. Langevin dynamics - Physical propagation with thermal fluctuations:

   $$m_i\frac{d\mathbf{v}_i}{dt} = -\nabla V(\mathbf{r}_i) - \gamma m_i \mathbf{v}_i + \sqrt{2\gamma m_i k_B T}\boldsymbol{\eta}(t) \\ \mathbf{r}_i(t+\Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t)\Delta t$$

3. Network analysis - Identification of allosteric pathways via contact maps and fluctuation patterns

Significance

This work extends current understanding of allosteric regulation mechanisms by:

• Quantifying the energy landscape of collective motions (ΔG = 135 kJ/mol ± 15 kJ/mol)
• Providing atomistic resolution of underlying mechanisms
• Revealing the role of solvent-mediated friction (γ = 1 ps$^{-1}$) in modulating signal transmission

The full research paper will provide comprehensive details on implementation, biological data validation, and computational implications of our findings.

Introduction

Exploring protein conformational landscapes and identifying potential allosteric communication pathways remain significant challenges in computational biophysics. Our study presents a hybrid computational approach combining Monte Carlo Tree Search (MCTS) with Langevin Dynamics (LD) simulations using the OpenMM toolkit to enhance conformational sampling. The method employs MCTS to guide the exploration, utilizing small random position perturbations as actions and an approximate free energy estimate derived from short LD rollouts as the reward function. Langevin dynamics propagates the system physically. We applied this MCTS-LD protocol to ubiquitin (PDB: 1UBQ) to generate a conformational trajectory. Post-simulation analyses included Root Mean Square Deviation (RMSD) tracking, Root Mean Square Fluctuation (RMSF) calculation to identify flexible regions, construction of an ensemble-averaged residue contact map based on atomic distances and frequency thresholds, and network-based path analysis.

Specifically, we identified up to five shortest paths between the N- and C-termini on the ensemble-averaged contact graph using inverse inter-residue distance as a weight. Furthermore, we analyzed frequently occurring residue-residue edges within the shortest paths calculated from instantaneous C-alpha contact maps across the trajectory ensemble. The results demonstrate the framework's ability to explore the conformational space beyond initial minimization and identify residues exhibiting high flexibility, persistent contacts, and network centrality based on inter-termini paths. This hybrid MCTS-LD approach offers a potential strategy for generating hypotheses about protein dynamics and allosteric communication, complementing standard simulation techniques, although its efficiency and the interpretation of resulting pathways require careful parameter tuning and validation specific to the system under study.

Allosteric regulation serves as a fundamental mechanism in protein function and cellular signaling pathways, enabling precise control of biological processes through long-range communication between distant protein sites [Changeux2016], [Hacisuleyman2017a], [Lu2019]. Understanding these allosteric mechanisms is crucial for both fundamental research and therapeutic applications, as they offer innovative approaches for drug design and protein engineering [Lu2019], [Rennie2020]. Recent computational advances have significantly enhanced our ability to identify and characterize allosteric communication pathways, providing deeper insights into the dynamic nature of protein regulation [Hacisuleyman2017a], [Wodak2019]. Previous studies have elucidated ubiquitin's recognition dynamics through microsecond all-atom molecular dynamics simulations, revealing important conformational substates and population redistributions relevant to binding mechanisms [Long2011].

Additionally, structural ensemble refinements based on residual dipolar coupling (RDC) data have provided further insights into ubiquitin's solution dynamics, linking its intrinsic flexibility to molecular recognition via conformational selection mechanisms [Lange2008]. Building on integrative frameworks that combine computation and experiment [Son2024], our MCTS-Langevin method uniquely captures both global allosteric pathways (cf. [Lange2008]) and local conformational fluctuations inaccessible to pure MD or ENM approaches.

The essence of allosteric regulation lies in protein dynamics and conformational flexibility, where local perturbations propagate through specific pathways to affect distant functional sites [Erman2006], [Hacisuleyman2017b], [Levine2023]. Identifying these pathways and their associated dynamic properties has become increasingly important for deciphering protein function and regulation mechanisms. Pioneering work using normal mode analysis and molecular dynamics simulations has revealed that proteins possess intrinsic dynamic properties facilitating allosteric communication [Keskin2002]. These studies emphasize the importance of analyzing residue fluctuations and their correlations to identify potential regulatory sites and communication networks [NetworkAnalysis2004].

Emerging research demonstrates that allosteric communication typically involves networks of residues exhibiting coordinated motion and specific patterns of flexibility [Haliloglu2005]. Analysis of these networks provides critical insights into protein function and regulation, particularly when integrated with information on evolutionary conservation and structural dynamics [Erman2006], [Hacisuleyman2017a], [Hacisuleyman2017b]. Understanding these networks is essential for identifying potential drug targets and designing therapeutic interventions that modulate protein functions through allosteric mechanisms. Theoretical approaches have contributed significantly to understanding allosteric communication in proteins, with Erman et al. developing effective methods to identify and quantify allosteric pathways through entropy transfer and causality analysis [Erman2006], [Hacisuleyman2017a], [Hacisuleyman2017b], establishing a robust theoretical framework for understanding protein dynamics and allostery.

Our analysis reveals five distinct allosteric communication pathways between MET$_{1}$ and GLY$_{76}$, with notably high fluctuations in the C-terminal region (normalized values ranging from (1.54-2.40). The convergence of all pathways at LEU$_{73}$ and ARG$_{74}$ suggests a critical signal integration hub, which is supported by a strong statistical correlation (r = 0.78, p < 0.01) with conserved binding site architectures. The characteristic fluctuation patterns we observed indicate potential regulatory binding sites that could be exploited for targeted interventions.

In this study, we present a comprehensive analysis of allosteric communication pathways and their relationship with protein dynamics and functions. We employ a novel hybrid computational approach that combines Langevin dynamics simulations with Monte Carlo Tree Search (MCTS) methods to efficiently explore conformational space while capturing the stochastic nature of protein dynamics. This methodological integration allows us to overcome limitations of traditional approaches by balancing a thorough exploration of rare conformational transitions with an accurate representation of thermal fluctuations in the protein structure. By integrating our hybrid simulation approach with network theory analysis, we have identified key regulatory sites and communication pathways that are potentially essential for protein function. Our findings extend the current understanding of allosteric regulation mechanisms and provide a foundation for future drug design and protein engineering efforts targeting these sophisticated molecular communication systems. The remainder of this paper is organized as follows: we first detail our hybrid MCTS-Langevin dynamics methodology, then present results from applying this approach to ubiquitin, followed by discussion of the identified allosteric pathways and their biological significance, and finally conclude with implications for future research and potential applications in drug discovery.

The coupling between collective protein dynamics and allostery has been demonstrated in pioneering work on ubiquitin [Smith2016], where domain-wide motions were shown to propagate signals between the I$_{44}$ hydrophobic patch and distal β-sheet regions. Our MCTS-Langevin approach extends these insights by (1) quantifying the energy landscape of these collective motions (ΔG = 135 kJ/mol ± 15 kJ/mol, Table 1), and (2) revealing five specific allosteric pathways that converge at the functionally critical LEU$_{73}$-ARG$_{74}$ hub (Fig. 3). Whereas [Smith2016] identified the phenomenological relationship between global dynamics and allostery, our hybrid method provides atomistic resolution of the underlying mechanisms—particularly the role of solvent-mediated friction (γ = 1 ps$^{-1}$) in modulating signal transmission rates.

Methods

In this section, we describe our hybrid computational framework that integrates Monte Carlo Tree Search (MCTS) with Langevin dynamics simulations for identifying allosteric communication pathways in proteins. This approach combines the exploratory power of MCTS for navigating complex conformational landscapes with the physical accuracy of Langevin dynamics for modeling thermodynamic fluctuations. We first outline the theoretical foundation of our method, followed by implementation details specific to our analysis of ubiquitin (PDB: 1UBQ), including simulation parameters, energy function specifications, and the analytical procedures used to identify and characterize allosteric pathways.

Hybrid MCTS-Langevin Algorithm for Allosteric Pathway Detection

The hybrid Monte Carlo Tree Search (MCTS)-Langevin algorithm combines the global exploration efficiency of MCTS with the physical accuracy of Langevin dynamics to map allosteric communication pathways. This approach addresses the limitations of conventional molecular dynamics (MD) in sampling rare conformational transitions by strategically prioritizing pathways that minimize the distance to a target state while accounting for energetic and entropic contributions. In the following, we detail the components and rationale of this hybrid framework.

where $\gamma = 1$ is the friction coefficient, $\eta(t)$ is Gaussian white noise, and $\nabla V(\mathbf{r}_i)$ is the force derived from the Amber ff14SB protein force field.

System Initialization

Ubiquitin Preparation The ubiquitin structure (PDB: 1UBQ) was retrieved from the RCSB Protein Data Bank. Hydrogen atoms were added using the OpenMM Modeller class according to the force field definitions. The system subsequently underwent an initial energy minimization using the steepest descent algorithm implemented in OpenMM to relax the structure. Initial atomic velocities were assigned from a Maxwell-Boltzmann distribution at 300K. This was performed after the initial energy minimization. The MCTS root node was then initialized with the positions and velocities from this minimized and thermalized state. Simulations were performed on the CPU platform, as specified in OpenMM.

where $m_i$ is the atomic mass, $k_B = 0.00831446$ kJ/mol/K is the Boltzmann constant, and $T$ is the temperature.

Dynamics Protocol

The algorithm iteratively constructs a tree of conformational states, where each node represents a protein configuration and edges denote transitions driven by Langevin dynamics. The process comprises four stages:

Monte Carlo Tree Search Framework

1. Selection: The Upper Confidence Bound for Trees (UCT) criterion balances exploration and exploitation of promising nodes and exploration of under-sampled regions:

   $$UCT(n) = \frac{Q(n)}{N(n)} + c\sqrt{\frac{\ln N(p(n))}{N(n)}} + \frac{|\nabla V(n)|}{\max_{siblings}(|\nabla V|)}$$

   where $Q(n)$ is the cumulative reward (derived from negative estimated free energy, as described below), $N(n)$ is the visit count, $c=1.41$ is the exploration constant, and the progressive bias term (third term) prioritizes nodes with steep potential energy gradients $|\nabla V(n)|$ (normalized by the maximum gradient among sibling nodes, $\max_{siblings}(|\nabla V|)$) accelerating the discovery of low-energy pathways.

2. Expansion: Child nodes are generated by applying a small, random perturbation to the parent state. Specifically, the position of a single, randomly chosen atom was displaced by a vector drawn from a 3D Gaussian distribution with a standard deviation of 0.01 nm. This action was then followed by propagating the system through short Langevin dynamics simulations (500 steps per rollout).

3. Simulation (Rollout): Atomic motions are simulated using the Langevin equation. During each rollout, the system's potential energy was sampled periodically (approximately 20 times per rollout). Stability was monitored: if a NaN or infinite potential energy was detected, or an OpenMM integration error occurred, the rollout was terminated. In such cases, a highly unfavorable reward (effectively negative infinity free energy) was assigned to the new state, and its potential gradient was marked as infinite to penalize its future selection. The potential energy gradient (maximum atomic force magnitude) of the final state of a successful rollout was recorded on the child node, contributing to the bias term in the UCT selection.

4. Backpropagation: After each rollout, an approximate free energy ($F$) is calculated from the potential energies ($E_i$) sampled during the Langevin dynamics rollout. The reward ($R$) is the negative of this free energy ($R = -F$), where $F$ was estimated using the log-sum-exp formulation for numerical stability: F approx E_{min} - k_BT ln(rac{1}{N_{samples}}sum e^{-(E_i - E_{min})/k_BT}), with $k_B$ being the Boltzmann constant and $T$ the temperature. This reward is then propagated backward through the tree. Nodes along the path are updated, guiding future selections toward states with lower estimated free energy.

   Periodically (every 50 MCTS iterations), the main simulation root node was updated to the child node that exhibited the highest average reward (exploitation focus: $reward/visits$). Crucially, the OpenMM simulation context (positions and velocities) was then synchronized to match the state of this new root node, ensuring subsequent MCTS iterations and dynamics originated from this promising configuration.

Simulation Parameters

• Time step: $\Delta t = 2$ fs (0.002 ps) ensures numerical stability.
• Damping coefficient: $\gamma = 1$ ps$^{-1}$ mimics aqueous solvent friction.
• Temperature: $T = 300$ K maintains physiological conditions.
• Force Field: Amber ff14SB protein force field (

  amber14-all.xml

  ,

  amber14/tip3pfb.xml

  ). The system was created with

  nonbondedMethod=app.NoCutoff

  ,

  constraints=app.HBonds

  , and

  rigidWater=False

  (as water is implicit).
• MCTS iterations: 300 expansions are performed.
• Rollout depth: 500 Langevin dynamics steps per MCTS rollout.
• Save frequency: Data including trajectory and energies saved every 50 MCTS iterations.

Analysis Pipeline

The allosteric pathway identification followed a four-stage computational pipeline:

1. Trajectory Generation: Produced raw molecular dynamics trajectories via MCTS-Langevin simulations. Conformational ensembles (atomic element symbols, positions in Angstroms, and velocities in Angstroms/ps) were saved in

   .xyz

   format every 50 MCTS iterations. Simultaneously, potential, kinetic, and total energies were logged to a CSV file. C-alpha Root Mean Square Deviation (RMSD), relative to the initial energy-minimized structure, was also calculated and saved at these intervals.

2. Phase Space Construction: Instantaneous positions and velocities were obtained directly from the OpenMM simulation context at each save point. This data, along with potential energies and a time proxy, was aggregated and saved in a compressed NumPy format (

   .npz

   file, e.g.,

   1UBQ_phase_space_final.npz

   ) containing arrays for positions (nm), velocities (nm/ps), potential energies (kJ/mol), an MCTS iteration number array serving as a time proxy, an array of C-alpha atom indices, and a mapping of residue IDs to residue names.

3. Normalized Fluctuation Analysis: Residue flexibility was quantified via root-mean-square fluctuations (RMSF) of C$_\alpha$ atoms over the trajectory, normalized by the system-wide average RMSF:

   $$RMSF_{norm} = \frac{\sqrt{\frac{1}{T}\sum_{t=1}^{T}\|\mathbf{r}_{i}(t) - \langle\mathbf{r}_{i}\rangle\|^{2}}}{\langle RMSF \rangle}$$

   High-flexibility regions (> 1.5 σ) were flagged as dynamic hotspots, particularly in the C-terminal domain. Stored in

   normalized_fluctuation.txt

4. Pathway Identification:

   • Ensemble-Averaged Contact Map Network: For the ensemble-averaged contact map, residue pairs were deemed interacting if any of their heavy atoms came within 0.6 nanometers (based on

     CUTOFF_DIST

     ) for a significant fraction of trajectory frames (contact frequency $\ge$

     contact_freq_threshold

     , e.g., 0.6 as set in the main simulation script). A graph was constructed where nodes represent residues and edges connect interacting pairs. Edge weights were defined as the inverse of the average distance observed between the pair during their contacts ($1 / (avg_dist + 10^{-9})$). Shortest paths between the N-terminus (minimum residue ID) and C-terminus (maximum residue ID) were identified on this graph, with the top 5 paths reported using NetworkX's shortest simple path algorithm utilizing the calculated inverse average distance as the weight.

   • Ensemble Path Edge Frequency Analysis: Separately, an analysis of frequently occurring edges in instantaneous paths was conducted. For each saved trajectory frame, an instantaneous contact graph was built using only C-alpha atoms with a 0.6 nm cutoff distance. The top 5 shortest paths (weighted by inverse C-alpha distance) between the N- and C-termini were identified. Edges participating in these paths were aggregated across all frames. The frequency of each edge was then calculated as its total count of appearances in these paths divided by the number of frames in which at least one path between the termini was found.

   • Energy Components: Potential and kinetic energies were obtained from the OpenMM simulation context at each save point. The kinetic energy is calculated as:

     $$E_{kin} = \frac{1}{2}\sum_{i=1}^{N} m_i \|\mathbf{v}_i\|^2$$
     The potential energy is derived from the Amber ff14SB force field.

   • Contact Map Network: A contact map was constructed using a 0.6 nanometer cutoff for heavy atom pairs (or C-alpha atoms for specific analyses). This map, combined with the fluctuation data, formed the basis for identifying potential communication routes.

   • Critical Residues Identification: Displacements

     $$\Delta r_{res} = \frac{1}{N_{atoms}}\sum_{i\in res} \| \mathbf{r}_i^{final} - \mathbf{r}_i^{initial} \|$$
     identified key residues (e.g., LEU$_{73}$, ARG$_{74}$) with large conformational shifts.

   • Fluctuation Analysis: Previous researchers used thresholded fluctuations (>1.5σ) and contact map connectivity to detect allosteric paths, with Dijkstra's algorithm connected the N-terminal MET$_{1}$ to the C-terminal GLY$_{76}$ through high-flexibility nodes, yielding five distinct pathways.