Finding multiple reaction pathways of ligand unbinding

Molecular dynamics (MD) simulations provide sufficient temporal and spatial resolution to study physical processes. Unfortunately, MD fails to reach high energy barriers (≫k_BT) that dictate mechanisms and kinetics of rare events. Transport in heterogeneous media, such as ligand unbinding, cannot be simulated directly, and even biased MD methods often fail to find possible reaction pathways along complex transient tunnels of proteins that form spontaneously during dynamics.^1–3 Although many general purpose methods have been developed to sample rare events,^4–6 finding multiple reaction pathways of ligand unbinding is especially difficult. Also, experimental methods used currently to quantify ligand binding, e.g., time-resolved crystallography and xenon binding, focus primarily on gaseous species, providing indirect evidence for the migration of larger ligands, which makes most details of reaction pathways unresolved.

The main computational limitations that render the reconstruction of reaction pathways for ligand unbinding difficult stem from accounting for internal topological features of proteins (e.g., tunnels), which is related to the degree of coupling between protein dynamics and ligand conformational states. The structural flexibility of protein tunnels allows proteins to facilitate binding by adapting to binding partners along possible multiple pathways to the binding site. This intrinsic dynamics poses a severe challenge to straightforward biased MD methods that have been used to sample reaction pathways in ligand unbinding.^6–9 Typically, such methods either approximate reaction pathways by linear Cartesian coordinates¹⁰ or probe protein tunnels randomly.^7,11

An additional and ubiquitous obstacle in describing ligand unbinding is the overestimation of energy barriers and thus the underestimation of exponentially dependent kinetic rates arising from the sampling of crude reaction pathways. In other words, an inadequate initial guess of reaction pathways leads to false thermodynamics and kinetics. Another problem is related to the degeneracy of microscopic configurations originating from the inability to capture intrinsic degrees of freedom, which is likely to shadow hidden energy barriers.^12,13 As emphasized by Elber and Gibson,¹⁴ sampling should not overestimate preference to more direct and geometrically shorter reaction pathways. Producing and exploring multiple reaction pathways of a complex system remains a huge challenge.¹⁵ 

In this letter, we consider a specific part of ligand binding/unbinding problem that is very relevant and not yet fully solved.¹⁶ To our knowledge, this is the first work to show that sampling multiple transient ligand tunnels in proteins leads to heterogeneous mechanisms of unbinding between the sampled reaction pathways. We present a general enhanced sampling MD method to find multiple diverse reaction pathways of ligand unbinding along transient protein tunnels. The method does not require many parameters and does not require initial guesses of intermediate states,^17,18 which is a major challenge for existing methods. Its only prerequisite is knowledge of the initial bound state without requiring the initial reactive trajectory. The method also takes into account protein dynamics, which is important to observe transient tunnels.

To estimate ligand-protein interaction, we introduce the concept of a loss function. The method minimizes a loss function $sx,y$ during MD simulations of a 3X set of ligand coordinates x ≡ (x₁, …, x_3X) and a 3Y set of protein coordinates y ≡ (y₁, …, y_3Y), where X and Y are the numbers of ligand and protein atoms, respectively. To this aim, we propose that the loss function must fulfill three important criteria, e.g., (i) describe physical interactions in a ligand-protein system, (ii) tend to infinity as the ligand moves too close to the protein, and (iii) decrease as the ligand unbinds from the protein; (ii) prohibits the method from sampling ligand configurations that clash with a protein, and (iii) provides a coarse estimate of how ligand conformations are buried within a protein tunnel.

For a schematic depiction of the method, see Fig. 1. The method follows a procedure: (i) it finds a minimum of the chosen loss function in the neighborhood of the current position of the ligand, and (ii) the position of the ligand is biased in the direction of the localized minimum of the loss function. The minimization is repeated during the MD simulation every Δt MD steps, and the biasing is performed until a new solution in the neighborhood of the current position is calculated. In what follows, we explain in detail the above general outline.

We start off by describing the loss function and the minimization procedure which provides a possible ligand configuration sampled in the proximity of the current ligand conformation from the MD simulation, which corresponds to the lowest loss function score, and by explaining how the neighborhood is defined for such an optimization problem. Next, we move on to the adaptive biasing procedure which enforces the ligand conformation to dissociate toward the conformation selected by the minimization procedure. The method is then summarized and used to model the benzene unbinding reaction pathways from the T4 lysozyme L99A mutant that is often used as a model system to study ligand unbinding processes.

Loss Function—Because empty space in proteins and its intrinsic fluctuations constitute a key feature of tunnels,^1,3 we use a coarse physical model for ligand-protein interactions, which accounts for steric effects only. We motivated our decision by the simplicity of this approach. For the ith pair of ligand-protein atoms, we define a partial loss function as $e−riri$⁠, where the rescaled distance between the atoms is given by r_i = λ∥x_k − y_l∥. The λ constant sets length scale in the loss function and is equal to 1 when using angstroms, i.e., λ = 1 Å⁻¹. Hence, we used the loss function of the following form:

where P_l is the number of ligand-protein atom pairs in the local neighborhood of the ligand (see the supplementary material for details). The sum over all pairs meets the criteria of the loss function for ligand unbinding presented here. The aim of the proposed method is to efficiently sample the configurational space of the ligand-protein complex and optimize Eq. (1) during MD simulations so that the reconstructed reaction pathways of ligand unbinding minimize the loss function along multiple tunnels. We also allow the ligand to be flexible during the unbinding simulations. In this method, many MD simulations are required to sample multiple reaction pathways.

Minimization—Such an optimization problem can be solved by any method suitable for nonconvex loss functions.^19–21 Here, for the sake of simplicity, the minimization of the loss function is performed using simulated annealing.²² To this end, the method checks if a randomly chosen neighboring position of the ligand x′ is preferred in terms of the loss function. The neighbor is selected as a next solution according to the Metropolis-Hastings algorithm²³ with the Boltzmann factors (we omit protein coordinates y only in the notation),

where β_j = 1/T_j is a parameter introduced to decrease the probability of acceptance of a worse solution as the minimization scheme proceeds. T_j is reduced according to the recursive formula, T_j = kT_j−1, where j is the iteration number during the optimization phase to promote convergence to an optimum.²⁴ The minimization procedure is reiterated to find an optimal solution. For details concerning the parameters for simulated annealing, see the supplementary material.

Next, we describe how the neighborhood is defined in our method. The minimization procedure needs constraints to optimize the loss function locally (in the current neighborhood). In our method, intermediate ligand unbinding states are searched sequentially to get an optimal transition between the X-ray structure and the unbound state. A global minimization of the loss function without a specific definition of the neighborhood would identify only the final state of ligand unbinding. A naive approach^20,25,26 is to sample ligand conformations constrained to a sphere with a constant radius and positioned at the center of mass of the ligand, but this requires an estimate of the radius, which is clearly system-dependent and should change as protein dynamics is simulated. To alleviate this issue, we take the sampling radius equal to the minimal distance between the ligand-protein atom pairs, e.g., r_s = min_ir_i. By doing so, the method dynamically adjusts the conformational space available for the sampling. We underline that the protein neighborhood of the ligand changes as the ligand dissociates during the simulation and so does the number of ligand-protein atom pairs P_l, which makes the identification of the next minimum possible.

Adaptive Biasing—Once the optimal ligand-protein conformation x′ = min_xs(x) is calculated in the minimization scheme, the conformation of the ligand is biased in the direction of x′ along transient protein tunnels. This stage is performed by biasing the conformation of the ligand using an adaptive harmonic potential,

where $xi′$ is the ith optimal solution, $v$ is the biasing rate, Δt is the MD time between subsequent loss function minimizations, and α is the force constant. The bias potential [Eq. (3)] is a generalization of the harmonic biasing potential introduced by Heymann and Grubmüller¹⁰ to curvilinear reaction pathways. The biasing potential from Ref. 10 uses a constant direction of biasing, but in Eq. (3), this direction is approximated as the normalized difference between the subsequent minima of the loss function. In contrast to this method, several recently introduced approaches used a constant bias to sample complex reaction pathways.^11,20,25 The bias potential shown by Eq. (3) is adaptive and dependent on the optimal reaction pathways calculated by minimizing Eq. (1).

The harmonic bias potential used in this study is selected to be simple as the potential energy in MD simulations already includes bonded terms for interaction of atoms that are linked by covalent bonds and nonbonded terms that describe long-range electrostatics and van der Waals forces. Clearly, the bias potential should not be expected to be quantitative as a method to calculate energy barriers along the reaction pathways, but it was employed to enforce the process of ligand unbinding with a constant velocity as in Ref. 10. The bias potential, however, may serve as a means to shorten the time scale of ligand unbinding and as a qualitative measure to estimate relative differences of bias between the reaction pathways.

Our enhanced sampling method is outlined as follows:

1. Initialize the MD simulation,

2. Sample ligand conformations within the protein tunnel using constraints defined as the minimal distance between the ligand-protein atom pairs,

3. Minimize the loss function using a nonconvex optimization algorithm, and set the biasing direction toward the found minimum,

4. Bias the ligand conformation using Eq. (3) during Δt steps of the MD simulation,

5. Repeat steps 2–4 during the MD simulation until the loss function reaches zero,

6. Stop the MD simulation,

which concludes the introduction of the method components, e.g., loss function, minimization, and adaptive biasing, needed to sample ligand unbinding reaction pathways.

Unbinding Benzene from T4 Lysozyme L99A—We illustrate the method on T4 lysozyme L99A (T4L) with bound benzene, which is considered as a model system to study ligand unbinding from proteins. In this example, 300 10-ns trajectories were run to reconstruct the reaction pathways of benzene unbinding from the protein. We used the biasing rate $v$ = 0.02 Å/ps with the force constant in the stiff-spring regime,¹⁰ α = 3.6 kcal/(mol Å). The optimal position of the ligand was recalculated by minimizing the loss function every Δt = 200 ps. We found that for lower biasing rates, the method is unable to find the reaction pathways in the desired span of 10 ns for a single simulation. This is, however, only a technical nuisance that can be overcome by sampling longer MD trajectories. The method is implemented in the official Plumed-2.5 repository²⁷ which is available on Github²⁸ and described in Ref. 29. For details concerning the model of T4L with bound benzene and the MD simulations, see the supplementary material.

We directly compared our results with reaction pathways found in previous studies. The method identified five reaction pathways for benzene exit from the binding cavity buried in T4L. These reaction pathways correspond to five tunnels of T4L, named pwa–D/F/G (tunnel through helices D, F, and G), pwb–C/D, pwc–F/G/H, pwd–H/J, and pwe–D/G (Fig. 2). The reconstructed reaction pathways pwa–d are mostly in agreement with a recent study by Nunes-Alves et al.³⁰ in which the reaction pathways of benzene unbinding were sampled using temperature-accelerated MD simulations.³⁰ Other studies also reported pwa⁸ and pwc.⁹ To our knowledge, the benzene unbinding via pwe is first identified in this study.

Apart from the work of Nunes-Alves et al.,³⁰ other studies found only one reaction pathway, probably because of the employed biased MD methods. Biased MD methods employed by Wang et al.⁹ and Miao et al.⁸ may limit the search in configuration space to a most optimal solution. Wang et al. used metadynamics³¹ to bias a reaction pathway identified initially by self-penalty walk,³² which agrees with the observation that such methods strongly rely on the initial guess of a pathway.^33,34 Thus, it may not be possible to identify all possible reaction pathways that exist in the form of transient sparse tunnels in the studied ligand-protein complex. Interestingly, the reaction pathways identified here agree mostly with exit tunnels retrieved by Nunes-Alves et al.,³⁰ where MD simulations with elevated temperature were used to overcome large energy barriers along reaction pathways and increase the probability of the rare event. In both temperature-accelerated MD³⁵ and the method presented in this article, there is no need for initial guess of trajectories, which clearly improves sampling of diverse pathways.

The detailed characteristics of the reaction pathways for benzene unbinding from T4L are shown in Table I. We found an additional reaction pathway that, to the best of our knowledge, was not identified previously. This pathway corresponds to the benzene unbinding along the T4L tunnel between helices D and G. The method is able to provide the atomistic characterization of unbinding pathways. If, however, one is interested in knowing estimates of the energy barriers along unbinding pathways, a postprocessing procedure is needed to analyze the data. For instance, one way to further understand the results is to look at the averaged bias potential V(s). We followed this approach and averaged the bias potential along each pathway projecting it on the loss function (Fig. 3) using the following relation:

where the average $⋅p$ is taken over all trajectories classified as a particular reaction pathway p and s is the loss function defined in Eq. (1). This way, we were able to identify energy bottlenecks in tunnels indicated by high values of the averaged bias potential V(s) or the sparsity of conformational space available for sampling. Treating the loss function as a collective variable, although may be not intuitive, provides a simple formula to check which transient tunnels are biased the most. Despite roughly the same level of the bias along each reaction pathway (Fig. 3), it is clear that the pathways employ different mechanisms of unbinding without the need to reconstruct free energies. This is underlined by the bias barriers along pwd and pwe and a rather smooth decrease in the bias along pwa, pwb, and pwc. Moreover, it is perhaps possible to explain the bias barriers by inspecting the average radius of each tunnel which is used as the sampling radius r_s. For instance, pwd and pwe have r_s at about 2.36 Å and 2.34 Å, respectively, and the highest barriers among pathways. This is an indication that the reaction pathways are heterogeneous with respect to each other, and their specific atomistic mechanism of unbinding would not be obvious by calculating averages of the full ensemble of the unbinding trajectories without decomposition into classes first.

As recently underlined in Ref. 36, multiple pathways for benzene escape from the T4L crystallographic binding site exist if one of the end states consists of multiple substrates.³⁷ Interestingly, the reaction pathway via the F/G/H tunnel identified by Feher et al.³⁶ is argued to be the most probable in their study. Our results show that this result may be due to the highest tunnel width in comparison with the other pathways (Table I). As it is shown in Table I, the sampled trajectories that yield the same reaction pathways (same tunnels) are similar to each other as it is underlined by the small standard deviations of the unbinding time and the sampling radius estimated using bootstrapping (see the supplementary material).

It should be noted that the method lends itself to use as an optimal initial guess of reaction pathways in other biasing MD methods to estimate thermodynamic and kinetic quantities, i.e., metadynamics^6,38 or variationally enhanced sampling.³⁹ We point out that computing reaction pathways for the T4L-benzene complex is not needed when calculating the mean-first-passage times of binding and unbinding as shown by Wang;⁴⁰ however, it is important in estimating how the mechanisms of binding vary between the calculated reaction pathways, including free energies and conformational changes. Recently, it was shown that some protein-ligand systems can exhibit pathway hopping,^41,42 and the method presented here can be used to quantify this process.

We note that the reaction pathways of ligand unbinding sampled using the method presented here diverge to diverse suboptimal basins. This is the feature that enables sampling multiple heterogeneous reaction pathways and allows us to overcome the problem of the intrinsic dynamics of protein tunnels. This is due to the used sampling which is constrained by the protein structure to provide a local minimum. Also, the probability of selecting a new solution given by the Metropolis-Hastings algorithm is important for the heterogeneity of the reaction pathways. The method searches for an optimal ligand conformation locally to extend the current reaction pathway step by step. This way, the method is able to sample multiple possible unbinding pathways, which for a rare event as with ligand unbinding is necessary to explore configurational space of tunnels exhaustively.

In conclusion, we have presented a general method for finding reaction pathways of ligand unbinding, starting only from available crystallographic information. The method does not need any prerequisite guesses of intermediate states. The introduced approach uses an adaptive bias to drive the ligand to unbind from the fluctuating protein, in the direction effectively calculated by minimizing a simple loss function. The methods adapt to transient tunnels of proteins by estimating the configurational space from which it samples plausible ligand conformations (i.e., it can be also used to determine the tunnel widths). We think that the method should be applicable to proteins in which prominent structural motions on a larger scale are important for ligand unbinding (e.g., trypsin⁴³).

Various enhanced sampling techniques have been tested for characterization of rare events and long-time scale dynamics. The method proposed here was suitable to sample a rare conformational event such as benzene escape that occurs on the millisecond time scale experimentally. We provided a rigorous method to find possible reaction pathways, which can be used as an initial reference trajectory to reconstruct thermodynamic and kinetic data. Overall, our results from studies of ligand unbinding from T4L suggest that the method presented here can improve the reconstruction of reaction pathways along transient tunnels and serve as an optimal choice for other biasing methods, limiting overestimation of hidden free energy barriers. With some adaptations, the method can be also used to study other transport processes, e.g., diffusion through a membrane.

Note. At the submission stage, we became aware of Ref. 44 in which the benzene unbinding pathways from the T4L protein are also studied. Capelli et al. found various reaction pathways that they classified into eight reaction pathways using different criteria than we use here. In particular, the benzene unbinding pathways marked by Capelli et al. as C, F, G, and H⁴⁴ are subclasses of pwc, while another four are the same as the ones we have identified here.

See supplementary material for the model of T4 lysozyme L99A, MD simulations, loss function, neighborhood for the loss function, minimization procedure, adaptive biasing to a loss function minimum, classification of the reaction pathways, biased unbinding times, and software.

We thank H. Grubmüller, W. Nowak, and M. Parrinello for useful discussions and Tristan Bereau and Claudio Perego for critically reading the manuscript. This work was supported by the National Science Centre, Poland (Grant Nos. 2016/20/T/ST3/00488 and 2016/23/B/ST4/01770). The MD simulations were computed using facilities of the Interdisciplinary Centre of Modern Technologies, Nicolaus Copernicus University, Poland.