Effects of Colored Noise in the Dynamic Motions and Conformational Exploration of Enzymes

Abstract

The intracellular environment displays complex dynamics influenced by factors such as molecular crowding and the low Reynolds number of the cytoplasm. Enzymes exhibiting active matter properties further heighten this complexity which can lead to memory effects. Molecular simulations often neglect these factors, treating the environment as a “thermal bath” using the Langevin equation (LE) with white noise. One way to consider these factors is by using colored noise instead within the generalized Langevin equation (GLE) framework, which allows for the incorporation of memory effects that have been observed in experimental data. We investigated the structural and dynamic differences in Shikimate kinase (SK) using LE and GLE simulations. Our results suggest that GLE simulations, which reveal significant changes, could be utilized for assessing conformational motions’ impact on catalytic reactions.

1. Introduction

The intracellular physical environment exhibits complex dynamics due to various factors, including crowding components and the low Reynolds number of the cytoplasm, among others [1,2]. Furthermore, this complexity is heightened by certain components, such as enzymes [3,4], which display active matter properties, meaning they possess internal driving forces [5]. Within this environment, a single molecule, for instance, a protein chain, is subject to continuous collisions with other components. Collectively, these factors can induce correlations between the motions of any molecule and its surroundings, resulting in memory effects [3,6].

This physical environment is generally not fully captured by molecular simulations in which some of the aforementioned factors are often neglected, such as the effects of the surrounding molecules, except for possibly water and ions. A standard approach to simulate the effects of the surrounding environment, viewed as a “thermal bath”, is through the Langevin equation (LE) with an external random force [7,8,9,10]. This random force is typically modeled as “white noise” where it is assumed that the interactions of the system and thermal bath are uncorrelated over time [10,11]. However, experimental observations indicate that the true physical environment plays a more complex role, for instance, in proteins, in which correlations exist among certain degrees of freedom, which are reflected in the behavior of time-correlation functions [12,13,14].

To incorporate memory effects in molecular simulations, an alternative to white noise can be employed, known as “colored noise” [10]. Colored noise has been utilized in various contexts, including thermostats for molecular dynamics (MD) [11,15], in the study of magnetic relaxation processes [16,17], enzyme cycle analysis [18], and simulations of biomolecules [19,20,21,22,23,24]. The introduction of colored noise leads to a generalized Langevin equation (GLE) [10,25].

The investigation of protein dynamics, particularly in enzymes, is critically important due to its potential role in catalysis—a topic that remains elusive [26,27]. The influence of the fluctuating environment on dynamics is typically incorporated into simulation studies via white noise [28,29], wherein, as mentioned above, the correlations are neglected a priori. Correlations can be crucial for enzymatic reactions as time-dependent observables, such as rate constants, display a wide range of variation [12]. For instance, these correlations can influence the hierarchy of functional motions [27] relevant during catalysis. Previous work reported the uncoupling of dynamic motions and catalysis [28]; however, the results were based on a specific type of noise satisfying the fluctuation–dissipation theorem (FDT).

In this work, we investigated the structural and dynamic differences in the Shikimate kinase (SK) enzyme when LE and GLE simulations are employed. SK is an enzyme involved in the production of chorismate, an essential compound for the functioning of pathogenic bacteria [30]. The simulation lengths considered in this work align with the timescales captured in current quantum mechanical and molecular mechanical (QM/MM) simulations [31]. Given that the LE and GLE simulations displayed structural and dynamic changes in SK, we propose that GLE simulations could be utilized to assess the influence of conformational motions on catalytic reactions in QM studies.

2. Materials and Methods

In this section, we outline the methodology employed in this work, specifically the Langevin equation (Section 2.1) and the generalized Langevin equation (Section 2.2), as well as the protocol followed for conducting the simulations (Section 2.3).

2.1. Langevin Equation (LE)

The LE is used to simulate the dynamics of a particle i, with a moment ${\mathbf{p}}_i$, in an environment modeled by a Stokes term with a friction coefficient $\gamma$. ${\mathbf{f}}_i$ and ${\zeta }_i$ are the deterministic and random forces acting on the particle, respectively, and expressed as follows:

$${\dot{\mathbf{p}}}_i={\mathbf{f}}_i-\gamma {\mathbf{p}}_i+{\zeta }_i.$$

(1)

In the LE, it is assumed that fluctuations from the heat bath, which is kept at a temperature $T$, occur on a shorter time scale than those of the particle itself [10,11,25]. Because of this, the former fluctuations can be modeled as white noise. This noise satisfies the following expression, which agrees with the FDT [9,32,33]:

$$〈{\zeta }_i\left(t\right){\zeta }_j\left(t^{\prime }\right)〉=2{\delta }_{ij}{m}_i{k}_{B}T\gamma \delta (t-t^{\prime }),$$

(2)

where $k_B$ is the Boltzmann constant, $m_i$ is the particle’s mass, and ${\delta }_{ij}$ is the Kronecker delta. The Dirac delta function $\delta (t-t^{\prime })$ reflects the fact that the impacts from the environment are almost instantaneous, and therefore they are not correlated in time.

2.2. Generalized Langevin Equation (GLE)

When the motions of the molecules in the system are correlated in time with the thermal fluctuations, the dynamics is better described using the GLE [25,34] as follows:

$${\dot{\mathbf{p}}}_i={\mathbf{f}}_i-\gamma {\int }_{-\infty }^{t}d\tau K\left(t-\tau \right){\mathbf{p}}_i\left(\tau \right)+{\eta }_i.$$

(3)

A key distinction between this expression and that of the LE (Equation (1)) lies in the friction term, as the former incorporates the history of the momentum weighted by a kernel function, as detailed in Ref. [35].

$$K\left(t-\tau \right)=2\mathsf{\delta }\left(t-\tau \right)-\frac{\lambda }{t_L}{e}^{-\frac{t-\tau }{t_L}}.$$

(4)

A similar kernel was used previously by Ceriotti et al. in the context of thermostats for MD [11]. Here, $t_L$ is the local average memory time that is used to smooth out high frequency while increasing low frequency fluctuations in the system. λ denotes the intensity of the guiding frictional force. Notably, this kernel contains the Ornstein–Uhlenbeck term (second term in Equation (4)), which vanishes when one considers long-time averages ($t_L\to \infty$) or minimal guiding force contributions ($\lambda \to 0$). The colored noise, ${\eta }_i$, also satisfies the FDT, akin to Equation (2) for LE, but it is related to the kernel in Equation (4) instead.

$$〈{\eta }_i(t){\eta }_j(t^{\prime })〉={\delta }_{ij}{m}_i{k}_{B}T\gamma K\left(t-t^{\prime }\right).$$

(5)

The fact that noise types described by Equations (2) and (5) satisfy the FDT guarantees that the systems are in equilibrium. Previous studies have employed various noise models to elucidate single-molecule experimental results [12,13,14]. However, it is important to note that colored noise is not limited to the type described here (Equations (4) and (5)); other types can also be explored that may lead to systems out of equilibrium.

2.3. Simulation Protocol and Analysis

2.3.1. Setting up the Simulations

SK from Helicobacter pylory (PDB code 3MUF, structure resolved at 2.3 Å) was used as a test case. It was studied previously through extensive MD simulations by one of the authors [36,37]. The protein structure was solvated with TIP3P [38] water molecules and salt, consisting of Na⁺ and Cl⁻ ions, at a concentration of 150 mM, resulting in a 66 ų cubic box. The setup of the simulation box and the input files’ generation for MD were carried out on the web-based graphical user interface for CHARMM [39] using the CHARMM-27 force field set of parameters [40,41,42]. Short-range interactions for vdW and electrostatic terms were computed explicitly up to a 12 Å cutoff distance. The former term was truncated by using a force-switching approach [43] starting at a 10 Å distance, while in the latter, the long-range interactions were solved with the fourth-order interpolation particle mesh Ewald (PME) method [44,45] with an FFT grid size of 72 in all directions. Simulations were conducted with the AMBER package (version 2022, San Francisco, CA, USA) [46].

The system was initially minimized during $10\times {10}^4$ steps, with the first 2500 cycles utilizing the steepest-descent method and the remaining steps employing the conjugate gradient method. After minimization, an equilibration procedure was performed in the NVT ensemble at 303.15 K using the Langevin dynamics thermostat with a friction coefficient of 1 ps⁻¹ during $1.25\times {10}^6$ steps (1 fs time step). Hydrogen bonds were constraint with the SHAKE algorithm [47] for the protein structure and the SETTLE algorithm [48] for the water molecules. Harmonic positional restraints of 1 kcal mol⁻¹ Å⁻² on protein atoms were applied during the minimization and equilibration steps.

The final structure of the equilibration step was used as the initial structure for data production for both LE and GLE simulations. In the former, Langevin dynamics was assessed in the NVT ensemble during 600 ns with similar MD parameters as in the equilibration step but with a larger time step of 2fs. Regarding GLE simulations, a local averaging time ($t_L$) of 0.2 ps and a momentum guiding factor ($\lambda$) of 1.0 were used. These simulations were also conducted in the NVT ensemble with similar parameters for the Langevin dynamics as in the LE case. Frames were saved every 0.02 ns for analysis with the first 100 ns being skipped and considered part of the equilibration step. In total, 3× repetitions for each case, LE and GLE, were run to obtain standard deviation bars.

2.3.2. Analysis of Simulations

Principal component analysis (PCA) [49] was performed by selecting the H, C, O, N, and C^α atoms with the AmberTools 2022 [50]. The first two PC values were used to generate 2D histograms of counts (heatmaps) that provide information on the explored conformational space using a bin size of 1 Ų. The ratios of areas were calculated by taking the number of bins explored divided by the total displayed area of $125\times 97$Ų. The root mean square fluctuation (RMSF) was computed with in-house Tcl scripts for VMD software (version 1.9.4, Urbana-Champaign, IL, USA) [51] based on the C^α atoms.

An interaction network was obtained, for both LE and GLE simulations, using an in-house Julia script (Code S1 in Supplementary Material), where C^α atoms formed a link if they lay within a cutoff distance of 7 Å, as suggested in previous works [52,53,54]. The weights for the interacting links were computed by taking the ratio of the frequency of link formation during the trajectory and the total number of frames. Based on these two computed networks an alluvial diagram was generated to compare the structural differences between them with Infomap (version 2.6.1, Umeå University, Umeå, Sweden), for community detection and then applying the mapping change algorithm [55,56,57,58]. In total, 100 trials with the two-level algorithm were used to obtain network communities [59]. Alluvial diagrams facilitate the visual comparison of the community structures across different networks. Infomap communities and alluvial diagrams were generated by using the online alluvial generator [60].

The normalized autocorrelation function (NACF) was computed according to the following equation (Code S2 in Supplementary Material):

$$\mathrm{N}\mathrm{A}\mathrm{C}\mathrm{F}\left(t\right)=\frac{〈\delta d\left(\tau \right)\delta d(\tau +t)〉}{〈\delta d\left(0\right)\delta d\left(0\right)〉},$$

(6)

Here, we consider the distance between the center of masses (COMs) of the LID (residues numbers 109–123), $r_{COM}^{LID}$, and SB (residues numbers 32–60), $r_{COM}^{SBD}$, regions, where $d\left(t\right)=\left|r_{com}^{LID}\left(t\right)-r_{com}^{SBD}\left(t\right)\right|$, $\delta d\left(t\right)=d\left(t\right)-〈d〉$, and $〈\cdots 〉$ denote time averages. The LID and SB domains, together with the link d connecting their COMs, are shown in the inset of Figure 4.

3. Results and Discussion

In the following subsections, we describe the structural and dynamic differences that were observed in the simulations for the LE and GLE simulations of SK.

3.1. Structural Modifications of SK in the Simulations with White and Colored Noise

Previous studies have indicated that the SB and LID domains of Shikimate kinase exhibit a high degree of fluctuations with and without substrates [36,37]. Changes in the fluctuations in these domains are important, as they may be linked to catalysis not only in Shikimate kinase [61] but also in other protein kinases such as the topologically similar [62] adenylate kinase enzyme [63].

Our observations indicate an increase in RMSFs in the SB region of SK for the GLE simulations compared to the LE simulations (see Figure 1). Conversely, the LID domain exhibited a slight decrease in fluctuations. Because both domains are involved in the catalytic step, through the binding of the substrates, we argue that the changes in their motion can be used to assess the influence of dynamics on the catalysis through QM-free energy simulations. This will be discussed in more detail in Section 3.3.

3.2. Differences in the Explored of the Conformational Space with the Two Types of Noise

The exploration of the conformational space was monitored through heatmaps, which were built upon the first two PCs; the results are shown in Figure 2. It was noticed that, in the LE case, the sampled conformational space was smaller (39.6% area ratio) than in the GLE case (48.9% area ratio). This was expected as the original purpose of the GLE method was to enhance the exploration of the conformational space.

In the LE simulations (see Figure 2a), four conformational basins can be distinguished, where the protein stays trapped most of the time. However, in GLE simulations (see Figure 2b), only one basin, around the (25,15) bin, is observed, which corresponds to the initial conformation. In this case, after the GLE gathers enough information on the low-frequency fluctuations, it escapes this basin and explores the conformational space more efficiently than the LE simulation. For systems with a few degrees of freedom, the sampled distributions for both LE and GLE simulations should converge to the canonical distribution, but in more complex systems, such as proteins, the relaxation process could extend over long time scales, and the distribution could differ in practical MD [64] and even more in QM simulations.

The difference in the conformational ensemble was also monitored through the analysis of residue network interactions by using alluvial diagrams. They allow for the detection of structural changes between networks [56,60]. In the alluvial diagram, the two interaction networks from LE and GLE simulations were compared, and the results can be seen in Figure 3.

In this diagram, communities are depicted as blocks (in different colors). When a community block from one network maps to another, it indicates that the two communities share the same network nodes. Otherwise, when a community block in the first network splits into two (or even more) community blocks in the second network, it means that some nodes from the original community interact more in the newly formed community from the second network. Thus, splitting reveals changes in the structure of the networks.

We noticed that the LE community formed by the residues ALA124, ARG112, and LEU118 from the LID domain was split into two community blocks in the GLE case. This was also the situation for the LE community of residues SER36, ASP30, and ILE35, which were in the SB domain. This reveals that the interaction network of residues in the LID and SB domains is modified in the GLE simulation compared to the LE case. Taken together, the alluvial diagram allowed us to find out where new interactions appeared in the protein networks and to visually present the differences in both networks.

3.3. Behavior of Time-Dependent Observables in White and Colored Noise Simulations

Besides the structural analysis, dynamic properties (in terms of NACFs) were computed. Because the LID and SB domains showed the largest RMSFs (see Figure 1) in the LE and GLE simulations, we computed the NACF of the distance between the COMs of these domains; the results are plotted in Figure 4. The plot shows that up to 10 ns, both systems behaved in a similar manner, with a slightly slower decaying behavior in GLE simulations. However, for longer times, they started diverging. It can also be noticed that beyond 100 ns, the NACF of the LE case decayed faster than that of the GLE case, indicating that a higher degree of correlation was maintained in the latter. We propose that the difference in these correlation rates between LE and GLE simulations may be used to assess the coupling between enzyme dynamics and the catalytic reaction in QM simulations. This may be achieved, for instance, by running both simulations and noticing the differences between energetics (through free energy differences) and dynamics (through the computation of correlation functions).

In these simulations, we observed differences in the NACFs between the LE and GLE cases, where the noise changed from white to colored. Now the question is whether there is any physical basis for colored noise. In experiments, colored noise was found to explain time-dependent correlation functions in enzymatic turnovers [12,13,14]. Also, computational studies suggest that the time-dependent transmission coefficient of the unfolded–misfolded amyloid-$\beta$ [23] and the viscosity dependence of folding–unfolding rates in proteins [21] are better described by using a colored noise in the GLE. Thus, there is already evidence supporting the idea that colored noise is present in physical systems such as enzymes.

4. Conclusions

In this study, we investigated the structural and dynamic differences in the Shikimate kinase (SK) enzyme when simulated using both the traditional Langevin equation (LE) with white noise and the generalized Langevin equation (GLE) with colored noise. The primary objective was to understand how incorporating memory effects, as opposed to treating the environment as a thermal bath with uncorrelated interactions, influences the structural and dynamic behavior of enzyme systems. By employing colored noise, which accounts for time-correlated interactions, the GLE can provide a more realistic model of the intracellular environment.

The simulations revealed notable structural and dynamic differences in SK between the LE and GLE approaches. Specifically, compared to the LE case, GLE simulations exhibited enhanced RMSFs in the SB domain, which is responsible for substrate binding during catalysis. The difference in the community networks of LE and GLE cases, monitored through alluvial diagrams, showed that mainly the residues in the SB and LID domains exhibited different interactions. Also, the NACF for the distance between the SB and LID domains displayed noticeable differences beyond 10 ns, which is a time scale achievable in QM simulations.

The ability of GLE simulations to capture these differences in the simulations suggests their potential utility in quantum mechanical and molecular mechanical (QM/MM) studies, where understanding the influence of conformational motions on catalytic events is still elusive. As opposed to using restraints on predefined collective variables in a QM/MM simulation to modify the behavior protein motions, one could run GLE simulations where those variables are not necessary.

The structural changes observed in the SK enzyme under GLE conditions indicate that traditional LE simulations may overlook critical dynamic aspects that are essential for enzymatic function.