Machine Learning Based Dimensionality Reduction Facilitates Ligand

Article pubs.acs.org/JCTC

Machine Learning Based Dimensionality Reduction Facilitates Ligand Diffusion Paths Assessment: A Case of Cytochrome P450cam J. Rydzewski* and W. Nowak* Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland S Supporting Information *

ABSTRACT: In this work we propose an application of a nonlinear dimensionality reduction method to represent the high-dimensional configuration space of the ligand−protein dissociation process in a manner facilitating interpretation. Rugged ligand expulsion paths are mapped into 2dimensional space. The mapping retains the main structural changes occurring during the dissociation. The topological similarity of the reduced paths may be easily studied using the Fréchet distances, and we show that this measure facilitates machine learning classification of the diffusion pathways. Further, low-dimensional configuration space allows for identification of residues active in transport during the ligand diffusion from a protein. The utility of this approach is illustrated by examination of the configuration space of cytochrome P450cam involved in expulsing camphor by means of enhanced all-atom molecular dynamics simulations. The expulsion trajectories are sampled and constructed on-the-fly during molecular dynamics simulations using the recently developed memetic algorithms [Rydzewski, J.; Nowak, W. J. Chem. Phys. 2015, 143 (12), 124101]. We show that the memetic algorithms are effective for enforcing the ligand diffusion and cavity exploration in the P450cam−camphor complex. Furthermore, we demonstrate that machine learning techniques are helpful in inspecting ligand diffusion landscapes and provide useful tools to examine structural changes accompanying rare events.

1. INTRODUCTION The process of ligand recognition by a protein is undoubtedly critical in biological signaling.1,2 Passing a signal entails the ligand binding into the protein docking site. The entrance and egress routes usually engage in labyrinthine migration through the channels or accessible cavities.3 The process of ligand dissociation is difficult to assess experimentally due to the complexity of protein structures, and in the absence of timeresolved crystallography experiments on ligand intermediates, the actual ligand expulsion pathways remain undetermined.4,5 Therefore, uncovering the distributions of transport routes has been studied by many enhanced molecular dynamics (MD) methods, i.e. locally enhanced sampling,6−12 targeted MD,13 steered MD,14 temperature accelerated MD,15,16 hyperdynamics,17,18 random acceleration MD (RAMD),19 supervised MD,20 metadynamics,21,22 and, recently, the memetic algorithms.23 In the memetic algorithms, the process of the ligand expulsion from the protein is reconstructed on-the-fly by the efficient conformational space sampling of the protein interior using immune computation techniques in such a way that the obtained trajectory minimizes the functional of the interaction free-energy between the ligand and protein.23 In the explicit all-atoms MD simulations of proteins, the changes in the positions of all the atoms in the protein as a function of time are calculated. Thus, N-atomic trajectories consist of an ordered set of 3N-dimensional vectors. Very often, © 2016 American Chemical Society

observing a particular phenomenon of interest is not possible, because the information on the interesting long-time events (e.g., nucleation, protein folding, dissociation) is buried under short-time events that obscure quantitative observation.24 The reduced probability distribution is a useful concept in the dynamics of the ligand diffusion process and, more generally, any rare event in statistical mechanics. Namely, instead of monitoring the full system trajectory, only a set of collective variables (CVs) is monitored.25 CV refers to any multidimensional function θ of 3N-dimensional atomic configuration x ≡ (xi|i = 1, ..., 3N). The functions θ1(x), θ2(x), ..., θM(x) map configuration x onto an M-dimensional CV space z ≡ (zj|j = 1, ..., M), where M ≪ 3N. Here, the CV space is referred to as a high-dimensional configuration space (HDCS). Typically, dimensionality reduction algorithms construct a low-dimensional configuration space (LDCS) representation of a set of data points distributed in a HDCS, by embedding objects in LDCS in a way that reproduces the pairwise distances between the objects in the HDCS. There are many techniques suitable for dimensionality reduction and manifold learning, i.e. multidimensional scaling,26 isomap,27,28 kernel principal component analysis,29 locally linear embedding,30 Hessian eigenmaps,31 diffusion maps,32 or t-distributed stochastic Received: March 18, 2016 Published: March 18, 2016 2110

DOI: 10.1021/acs.jctc.6b00212 J. Chem. Theory Comput. 2016, 12, 2110−2120

Article

Journal of Chemical Theory and Computation

Figure 1. Atoms used to calculate collective variables in (a) P450cam and (b) camphor. All the 12 atoms used to define the high-dimensional configuration space are colored in red (P450cam and camphor) and in orange (heme), and their PDB indices are available in the Supporting Information.

neighbor embedding (t-SNE).33 Therefore, providing a map of reduced dimensionality spanning the configuration space of a given process is of high importance; for instance, this idea was applied in calculating free-energy landscapes of protein folding.34 Recently, the sketch-map algorithm has been proposed to surpass difficulties with interpretation and choose an appropriate set of coordinates embedded in LDCS in the molecular dynamics.24,35−37 In this paper, we propose for the first time to adopt the tSNE nonlinear dimensionality reduction method to represent the protein−ligand dissociation process. To illustrate our approach, we applied this method to camphor migration in cytochrome P450cam. We focused on the exploration of the configuration space of P450cam by probing camphor positions using the memetic algorithm during molecular dynamics simulations. The resulting process of ligand diffusion is not easy to visualize and interpret, even in the CV space, since the configuration space of the P450cam−camphor system is highdimensional. To facilitate analysis, we propose to use the t-SNE manifold learning algorithm to generate LDCS as a 2dimensional map of configuration space. In contrast to the standard representation of diffusion paths, the new representation accounts for both ligand and protein conformational changes in time. The proposed approach enables identification of configuration similarities between the resulting pathway classes and the main attraction basins to which the expulsion trajectories sink. The 2-dimensional map offers numerous possibilities of extended analysis appropriate for molecular systems. There are several physical insights that can be drawn from the presented form of analysis: (a) the ligand diffusion process incorporates the configuration space of both the ligand and protein, because their flexibility often impacts diffusion routes via, for instance, gating mechanisms; (b) the Fréchet metric is particularly helpful in clustering exit routes of a ligand in a fully automatic manner, thus, computer graphics analysis is not necessary; (c) from the low-dimensional configuration space of the ligand−protein dissociation process, it is possible to produce a network of the most perturbed protein residues by the diffusing ligand in terms of root-mean-square distance, and thus exploit information about the gating mechanisms important in the active and inactive conformations of the studied protein; (d) the percentage of protein residues which actively participate in the ligand diffusion can be determined.

2. METHODS 2.1. Construction of the low-dimensional configuration space. The ligand expulsion trajectories were found using on-the-fly minimization of the free-energy interaction functional ΔG between P450cam and camphor during molecular dynamics simulations. Finding the minimum of ΔG was performed for an each intermediate binding position, beginning from the initial docking pocket of camphor in P450cam. Camphor was pulled along the minimum ΔG path using the constant force steered MD scheme.23 In the latest work we used multiple sampling algorithms capable of minimizing ΔG for an each intermediate binding position of camphor to P450cam along the diffusion route. But here, we used the best algorithm from ref 23, the memetic algorithm with immune evolution scheme and random mutation hill climbing38 (IA-RMHC) as a sampling procedure. This algorithm is described in the Supporting Information (Section S1: The memetic algorithm). To reconstruct the main features of the high-dimensional landscape of P450cam, we gathered a large set of CVs from 30 expulsion simulations given by the memetic algorithm. We selected CVs that describe not only the position and orientation of camphor relative to P450cam, but also the P450cam flexibility during MD simulations.39 The chosen CVs are based on the distances, di, between the atoms of the ligand (Figure 1b), the Fe atom in the heme group, and 7 arbitrary spaced points on the protein surface (Figure 1a), designated based on the experimental data on the ligand binding intermediates in P450cam,40,41 experiments of xenon binding to P450cam at high pressure,42 and theoretical work on finding expulsion pathways in the cytochrome P450 family using RAMD.19,42−45 The list of selected atoms is available in the Supporting Information. The distances, di, were transformed by a switching function, si, in order to suppress a possible influence of camphor being expulsed far from P450cam on CVs. The following switching function was used: 4

() s = 1−( ) 1−

i

di d0

di d0

8

(1)

where di is the ith distance and d0 = 20 Å. After the application of the switching function, CVs span 66-dimensional config2111


Article

Journal of Chemical Theory and Computation uration space representing MD trajectories of the expulsion of camphor from P450cam (HDCS). Interpreting a high-dimensional, maze-like landscape is not an easy task. In order to understand better the nature of expulsion pathways, we aimed at reduction of the dimension of configuration space. We used a machine learning algorithm for nonlinear dimensionality reduction developed by van der Maaten and Hinton: tdistributed stochastic neighbor embedding.33 This technique is particularly helpful for embedding high-dimensional data into a space of 2 dimensions, which can then be visualized, for example, in a scatter plot. In other words, t-SNE reduces the dimensionality of a space by representing similar objects by nearby points and dissimilar objects by distant points in a lowdimensional space. The t-SNE algorithm has two steps. First, tSNE constructs a probability distribution over pairs of highdimensional objects, pij, in such a way that similar and dissimilar objects have a high probability and low probability, respectively, of being chosen. The probability of the high-dimensional objects, zi and zj is the following:

pj|i =

⎛ − || zi − z j ||2 ⎞ ⎟ exp⎜ ⎝ 2σi2 ⎠

(

∑k ≠ i exp

− || zi − z k ||2 2σi2

and

)

pij =

ΔG ∝

i≠j

Bij ⎤ ⎥+ rij6 ⎥⎦

qiqj ϵ(rij)rij

⎡B

∑⎢

ij 12 ⎢ r i < j ⎣ ij

−

Cij ⎤ ⎥ rij10 ⎥⎦

+ ∑ (SiVj + SjVi ) exp( −rij2/2σ 2) i1 kcal/mol while camphor is leaving the protein), which means that the percentage of P450cam atoms surrounding camphor is smaller than in the binding pocket. In other words, a high value of Γ suggests a successful expulsion of camphor to the P450cam exterior. On the other hand, a small value of Γ indicates either an ongoing expulsion simulation or a trajectory which ends with camphor not being expulsed from P450cam. The choice of Γ seems to be somewhat arbitrary, but calculating this quantity gives us the most helpful results in determining the current stage of the diffusion path. Neither RMSD nor ΔG alone provides a strong distinction of the diffusion stage. A comparison of Γ with RMSD or ΔG is given in the Supporting Information (Figure S1). 2.2. Pathways similarity analysis. To quantify the topological similarities between the P450cam−camphor system trajectories Ti and Tj represented on the LDCS Euclidean metric S, we calculated the Fréchet distances F(Ti,Tj) from the following formula:46

2N (2)

(3)

The t-SNE algorithm minimizes the nonsymmetric Kullback−Leibler divergence, KL(P∥Q), between the low-dimensional probability distribution, Q, and the high-dimensional probability distribution, P, with respect to the locations of the points in the space:

∑ pij log

−

(6)

(1 + || yi − yj ||2 )−1

KL(P || Q ) =

∑ i1.5 Å mol/kcal). Expulsion paths which do not satisfy this criterion were grouped into one cluster describing trajectories trapped in the interior of P450cam (Figure S2, 12 trajectories in aqua). It is necessary to mention that these interior trajectories were not artifacts, since the memetic algorithm was used to probe the protein interior. The rest of the trajectories were used to perform the agglomerative clustering in LDCS. Since the last stage of each MD expulsion trajectory was the most important, we used 15 last structures from each expulsion to calculate the Fréchet distances. The paths were clustered automatically into 4 groups shown in Figure 3: PW1 (5 trajectories in blue), PW2 (6 trajectories in orange), PW3 (5 trajectories in green), PW4 (2 trajectories in yellow). Some of these pathways (PW1−3) were already analyzed by means of theoretical and experimental methods.19,41−45 The Pearson product-moment correlation coefficient showed total positive correlation (r = 1) between

4. RESULTS AND DISCUSSION 4.1. Configuration space of P450cam. The low-dimensional configuration space of P450cam (Figure 2) consists of more than 3000 points acquired using MD simulations with the memetic algorithm for enhanced expulsion of camphor. These points are colored due to Γ, which indicates the current stage of a simulation and the result, if the simulation is finished. The shape of the objects distribution in LDCS is circular, since we used the switching function to transform CVs. This transformation was applied in order to threshold the impact of camphor migrating in the P450cam exterior, since the main focus of this study was to examine the internal maze-like configuration space of P450cam and the expulsion routes of camphor. All the expulsion trajectories of camphor start in the initial docking pocket of P450cam, which is indicated in Figure 2. We call it the kernel of P450cam LDCS. The memetic 2114


Article


other steric barriers caused by, i.e., heme and helices B′, C, F, G, and I (Figure 4). The characteristics of the pathways are as follows. Camphor on PW1 leaves the docking pocket of P450cam approximately on the distal side of the heme plane (Figure 2). PW1 involves overcoming an energy barrier corresponding to distortion of the heme group and crawling through a region consisting of helix I characterized by low flexibility40 (Figure 6). This region of PW1 is followed by a region of high flexibility (helix C′). PW1 was also predicted by random acceleration molecular dynamics.19,42 There is evidence that PW1 may have a functional role in transport. Opening of PW1 was observed in CYP51.42 Its location near a high-flexibility region in P450cam also suggests a potential transport pathway for small species. Moreover, a xenon binding site was observed near PW1 which suggests a possible egress pathway for gaseous species.42 We observed that all the ending points on LDCS belonging to PW1 are not monotonically increasing in Γ, which indicates large fluctuations of the interaction free-energy between camphor and P450cam, ΔG, along PW1. Thus, the expulsion through a helix I of small flexibility on PW1 causes a sudden rise in ΔG. Perhaps, by applying a dragging force constant smaller than 5 kcal/mol Å in the memetic algorithm, the probability of hitting PW1 would decrease. However, we wanted to probe all the possible conformational exit routes for camphor, and by decreasing the dragging force constant, we would bias the results to the most probable exit pathways. That was not our intention. The residues approximate to PW1 are the most highly perturbed by expulsing camphor in terms of RMSD. Our analysis showed that there are a number of residues that collide with camphor, giving RMSD of more than 4 Å, which is a high value compared to the other pathways (Figure 5). PW1 incorporates 3% of all connections (Figure 5) between all the residues active in the transport. They are depicted in a classical way in Figure 6a (Thr252, Leu362, Leu204, Leu166. Ile160, Ile162, Phe163 that belong to C′, I, G, L helices). Taking into consideration the results, PW1 is perhaps not preferred for camphor either energetically or topologically. The region involving migration on PW2 is variable in sequence and structure between different cytochromes in the P450 family.42,44 Its high flexibility in all P450s suggests a role

the ligand diffusion pathways (PW1-4) observed on graphics and grouped by the agglomerative clustering (Figure 3). After clustering, we calculated a mean distance within the each cluster and a mean distance between all the cluster pairs (Table 1, Table 1. Average Fréchet Distances [a.u.] Calculated for Pathway Classes PW1−4 Determined by the Agglomerative Clustering: within each Pathway Class ⟨dii⟩ and between All Pairs of Classes ⟨dij⟩ Pathway class

⟨dii⟩

⟨dij⟩

PW1 PW2 PW3 PW4

5.3 7.8 5.6 3.9

17.0 18.7 12.4 12.3

Figure 3). The cluster with PW2 trajectories was the most distant from the other clusters. The results show that clustering on LDCS with the Fréchet measure can be used as a tool for comparing MD trajectories. Therefore, the nonlinear dimensionality reduction of the maze-like configuration space of P450cam using the t-SNE method opens new possibilities of the MD analysis. More importantly, the classified main expulsion pathways are distinguishable by means of clustering and relatively easy to interpret from the LDCS map in Figures 2 and 3. 4.2. Camphor expulsion trajectories. The crystal structure of P450cam has 13 α-helices (A, B, B′, and C−L) and five β-sheets (β1−β5). We have taken a departure from the structural nomenclature used in ref 38. Figure 4 depicts αhelices which are in the close proximity of ligand diffusion pathways PW1−4 (B′, C, F, G, and I). There are 4 main camphor expulsion pathways from P450cam (PW1−4) that are distinguishable by means of the graphical analysis of MD simulations, the nonlinear dimensionality reduction to LDCS, and the agglomerative clustering (Figures 2 and 3). Each pathway class is defined by a cluster of trajectories (swarm) that leads to the P450cam exterior, showing varying extents of diversity (Figure 3, S2). The swarm of trajectories sinks to different attraction basins, which are separated due to the secondary structure of P450cam, the shape of camphor, and the

Figure 4. Cytochrome P450cam fold. The secondary structure elements, α-helices, in the close proximity of ligand diffusion pathways PW1−4 are labeled with letters. The heme group is omitted for clarity. 2115


Article


Figure 5. Network of P450cam 238 residue perturbations that frequently occur during camphor diffusion. The size and color of nodes correspond to the average RMSD of P450cam residues calculated from 30 enhanced MD simulations with respect to the initial docking pocket of camphor (S). The edges are colored corresponding to the 4 clustered trajectory groups: PW1, blue; PW2, orange; PW3, green; PW4, yellow. Distances between the nodes correspond to low-dimensional configuration space. All the residues perturbed by more than 2 Å are mapped from LDCS to the network, and those with the average RMSD bigger than 4 Å are additionally labeled by a residue identification number. The average RMSD of the nodes belonging to each path is shown in the legend. Moreover, the percentage of P450cam residue connections used in each diffusion pathway is also present in the legend.

results are in agreement with the identification of the PW2 exit region (BC loop/B′ helix) as a putative substrate recognition site.42 Moreover, our results showed also topological indications that PW2 is preferred in substrate transport in P450cam. The residues accompanying transport along PW2 are perturbed by average RMSD approximately of 2.6 Å, which is a small value compared to PW1 and PW3 (Figure 5). Also, there is no single amino acid which is affected in terms of RMSD by more than 4 Å. The cluster formed by PW2 showed that camphor can easily migrate in P450cam along this pathway, since there is plenty of internal space leading to the PW2 exit. Also, PW2 uses approximately 5% of connections between active residues (Figure 5). Other pathways are heavily guarded by two amino acids, Leu244 and Val247, which belong to the region of small flexibility in P450cam−helix I (Figure 6). Our results are consistent with experimental records, because PW2 was designated to be preferable both energetically and topologically.43 The third pathway in P450cam corresponds to the so-called water channel (Figure 6c). PW3 is a hydrophilic channel located near the heme propionic groups. The existence of PW3 has been already suggested as the exit channel for the product, 5-hydroxy-camphor, which is more hydrophilic than camphor.40,41 PW3 was not observed in P450s cytochromes in standard RAMD calculations.19,42,43 PW3 opens toward the

in selectivity of a substrate. Moreover, the end of PW2 is known to have the most wide cavity in P450cam, since the other channels are blocked by the regions of small thermal motions (helix I) and the heme group. Because of that, there are many theoretical results suggesting that PW2 is most likely subdivided into classes (based on the proximity to F/G loop and BC loop/B′ helix), that vary from 2 in CYP102A1 to 3 in P450cam.42 Such splitting of PW2 into subclasses is perhaps related to the nature of the RAMD algorithm used in the literature previously. RAMD is a random walking technique and, because of that, a higher diversity in the resulting expulsion trajectories is observed. This is particularly troublesome in telling which PW2 subclasses are more probable. The sampling method used in this study provided the most probable expulsion trajectory in PW2, since, in contrast to RAMD, it is based on the energetic measure of the preferable ligand position in the protein. Thus, the resulting expulsion trajectories that are clustered in PW2 (Figure 3) are convergent and similar.23 The expulsion trajectories acquired during MD from PW2 are characterized by the highest Γ in comparison with the other pathways. For a globular protein, this means that PW2 is the most energetically preferred expulsion pathway in P450cam, which confirms the postulated functional role of this pathway in a substrate recognition between different P450s. In Figure 2 it is shown that all the expulsion trajectories clustered to PW2 converge to the biggest attraction basin on LDCS. Our 2116


Article


Figure 6. Detailed scheme showing camphor diffusion pathways from P450cam found by the memetic algorithm and agglomerative clustering in the low-dimensional configuration space. The active residues (shown as sticks) identified within the network are shown for (a) PW1, (b) PW2, (c) PW3, and (d) PW4.

proximal side of the heme group, which suggests that it may be involved in the electron transport system42 or in the transport of protons through a network of ordered water molecules.45 Another study has shown that PW3 is open if camphor is on the distal side of the heme.41 Since expulsion started from the distal side and the memetic algorithms offer an enhanced conformational sampling, 30 simulations were sufficient to identify and probe PW3. All the expulsion trajectories clustered into PW3 show increasing Γ and RMSD; therefore, the interaction free-energy, ΔG, is decreasing, albeit nonmonotonically. This suggests an energetic steric barrier that is imposed by bending heme to access PW3 (Figure 2). Once the hydrophilic PW3 exit route is accessible, camphor migration through it is not preferred energetically, due to the low hydrophilicity of camphor. PW3 is characterized by the highest average RMSD (2.85 Å) compared to the other pathways (Figure 5), suggesting that transport through a channel near the heme priopionic groups would be preferred by a smaller substrate than camphor. PW3 incorporates several

amino acids, which are highly affected by camphor diffusion, such as Leu244, Leu245 (helix I), Thr349 (L/β1 loop), and Leu200 (helix G). Arg299 has a large RMSD (Figure 5) but is isolated from PW3 on the network. Thus, LDCS prompts us to analyze the role of this residue in PW3 using computer graphics. We observed that camphor migration along PW3 is strongly connected with deviating helix I, either by camphor or by heme. Subsequently, a gating mechanism is triggered by twisting Arg299 (Figures 5 and 6c). Interestingly, Arg299 is not deviated directly by camphor, but through heme. Therefore, there is no an active connection in the LDCS network (Figure 5), as noticed before. This suggests that Arg299 is bended before camphor arrives, in close proximity of the L/β1 loop. Moreover, PW3 uses 3% of all topological connections in P450cam that are active in camphor migration (Figure 5), which suggests that transport on PW3 is a very rarely populated pathway. We found also some expulsion trajectories that were not observed in the previous studies. These were grouped by 2117


Article


We demonstrated, using as a challenging example the classical P450cam−camphor system, that the memetic algorithm can be used as a technique of efficient probing of internal space of a protein by a ligand. Here, the memetic algorithm gave not only the suboptimal exits (in terms of the interaction free-energy) of P450cam as we presented in ref 23, but its internal channels and cavities as well. One cluster of pathways, PW4, was not discussed in the literature before, and this new observation may be used as a hint in further experimental studies. A swarm of trajectories from PW2 sinks into an attraction basin much wider than those of the other PWs. The free-energy differences also suggest that PW2 is energetically the most preferred pathway. PW2 is known to be the main substrate recognition pathway, and the memetic algorithm is in accordance with the previous observations. Nonlinear dimensionality reduction is a general approach; for example, it may be further used in elucidating differences between expulsion trajectories from different P450s. It may help in visualization of the impact of particular amino acid mutations on the dissociation process. We plan to investigate these problems in our laboratory. The memetic algorithm and the proposed nonlinear dimensionality reduction coupled with the proposed visualization of expulsion trajectories can be used not only in ligand dissociation studies. We are confident that our approach is general and may be applied in modeling of any transport phenomena observed in nature.

agglomerative clustering into pathway PW4. Only two trajectories emulated characteristics representative to this pathway. In other words, PW4 is the least populated pathway cluster found in this study and, because of that, results might be not conclusive. These trajectories are characterized by Γ equal to 1.9 Å mol/kcal and 1.1 Å mol/kcal; their shapes on LDCS are more complex in comparison to the most probable PW2, but definitely simpler than PW3 (Figures 2 and 3). The agglomerative clustering revealed that trajectories from PW4 are very similar (⟨d44 = 3.9⟩), although the mean Fréchet distance to the other clusters is small (⟨d4j = 12.3⟩). Our calculations revealed that camphor migrating on PW4, alongside the low-flexibility region (helix I) in P450cam, is not blocked by heme nor any other steric barrier. PW4 is also characterized by small RMSD fluctuations (

Machine Learning Based Dimensionality Reduction Facilitates Ligand

Recommend Documents