Role of Protein Dimeric Interface in Allosteric ... - ACS Publications

Jul 24, 2017 - Search for allosteric sites using the Allosite web server and SiteMap analysis allowed us to identify substrate ... pointed areas at th...
0 downloads 12 Views 9MB Size
Download PDF
Article pubs.acs.org/jcim

Role of Protein Dimeric Interface in Allosteric Inhibition of N‑AcetylAspartate Hydrolysis by Human Aspartoacylase Ekaterina D. Kots,†,‡ Sofya V. Lushchekina,‡ Sergey D. Varfolomeev,†,‡ and Alexander V. Nemukhin*,†,‡ †

Department of Chemistry, Lomonosov Moscow State University, Moscow 119991, Russia Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, Moscow 119334, Russia

‡

S Supporting Information *

ABSTRACT: The results of molecular modeling suggest a mechanism of allosteric inhibition upon hydrolysis of N-acetyl-aspartate (NAA), one of the most abundant amino acid derivatives in brain, by human aspartoacylase (hAsp). Details of this reaction are important to suggest the practical ways to control the enzyme activity. Search for allosteric sites using the Allosite web server and SiteMap analysis allowed us to identify substrate binding pockets located at the interface between the subunits of the hAsp dimer molecule. Molecular docking of NAA to the pointed areas at the dimer interface predicted a specific site, in which the substrate molecule interacts with the Gly237, Arg233, Glu290, and Lys292 residues. Analysis of multiple long-scaled molecular dynamics trajectories (the total simulation time exceeded 1.5 μs) showed that binding of NAA to the identified allosteric site induced significant rigidity to the protein loops with the amino acid side chains forming gates to the enzyme active site. Application of the protein dynamical network algorithms showed that substantial reorganization of the signal propagation pathways of intersubunit communication in the dimer occurred upon allosteric NAA binding to the remote site. The modeling approaches provide an explanation to the observed decrease of the reaction rate of NAA hydrolysis by hAsp at high substrate concentrations.

■

INTRODUCTION Human aspartoacylase (hAsp) is an enzyme responsible for hydrolysis of one of the most abundant amino acid derivatives in the brain, N-acetyl-aspartic acid (NAA). Although a mechanism of regulation of NAA level has not yet been established, the necessity to maintain NAA concentration in relatively narrow limits is commonly accepted: low levels of NAA are observed in cases of multiple sclerosis, Alzheimer’s disease, and epilepsy, whereas high levels are correlated with Canavan disease, a fatal neurological disorder. By these reasons, investigations of the hAsp-catalyzed hydrolysis of NAA have gained considerable attention in recent years.1−4 The chemical reaction catalyzed by hAsp is cleavage of the amide bond in NAA producing acetate and L-aspartate:

could be attributed to presence of additional substrate binding sites.5 Initially, our aim was to propose a molecular mechanism of allosteric regulation of substrate hydrolysis by hAsp. We applied modern modeling tools to prove that the experimental results could be explained by presence of a remote substrate binding site at the interface of hAsp subunits of the protein dimer. In the course of this work, we also found that molecular simulations showed coupling of allosteric regulation, which takes place at the dimeric interface, with dynamical communication of the monomers across the dimeric interface. In our analysis, we rely on several available crystal structures of hAsp6−8 from the Protein Data Bank.9 The results of experimental studies8 show that hAsp appears as a functional homodimer. Figure 1 illustrates the protein dimer composed of two monomers (also called here as subunits) denoted in this paper as ASPA and ASPB. The atomic coordinates are taken from the structure PDB ID: 2O53,6 which corresponds to the apo-form of the enzyme. The active site in each subunit contains the Zn2+ ion, its coordination sphere (His21, Glu24, His116), and several amino acid residues including Glu178, which directly participate in the hydrolysis reaction. In our previous quantum mechanics/molecular mechanics and molecular dynamics based simulations,10 we constructed computationally the reaction energy profile for the complete catalytic

The parameters of the Michaelis−Menten kinetics (kcat = 12.7 ± 0.05 s−1, Km = 0.12 ± 0.03 mM) were reported, in particular, by Le Coq and co-workers.5 These experiments, hypothesizing co-operative effects, also revealed substantial substrate self-inhibition at high NAA concentrations, which © 2017 American Chemical Society

Received: March 4, 2017 Published: July 24, 2017 1999

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Figure 1. Homodimer of hAsp constructed by motifs of PDB ID: 2O53. The monomers, ASPA and ASPB, are shown here and in other figures in different colors, light green and light violet, respectively. Insets show coordination spheres of zinc and major gate-forming residues Arg71 and Glu293. At the bottom, we show the NAA molecule (balls and sticks); the substrate entrance channel to the active site of the ASPA monomer is indicated by the pink arrow.

for the prediction of allosteric sites in proteins available through a web server, Allosite,21 based on the Allosteric Site Database (ASD v.3.0),22 was used for initial screening. Allosite2015 scoring functions are adopted from the geometry-based f pocket algorithm23 and can be used for cavity screening in both holoand apo-protein structures. The SiteMap cavity detection algorithm24 was applied to extend the results obtained with Allosite. The SiteMap program is an energy-based tool aimed to identify binding sites and to range their druggability. SiteMap scoring functions, SiteScore and DScore,24 include contributions accounting for adequate shape and size of the binding pocket, tightness of the site as well as a penalty function imposed on increasing hydrophilicity. The sites predicted to be druggable should have the highest scores. Next, we used the Autodock4.2 program25 to search through the areas predicted by Allosite and SiteMap and to find location of NAA on the protein surface with the largest binding energy. Autodock4.2 algorithm is based on the validated semiempirical force field parameters allowing one to predict association energies of small organic molecules to macromolecule targets. In this application, the dimer protein interface area was split into five rectangular grids with approximate dimensions of 45 × 50 × 45 Å3. The Lamarckian Genetic Algorithm reinforced by Solis and Wets local search procedure26 was applied to each grid with the parameters set as follows: population size 10 000, generations number 2.7 × 105, evaluation number 2.5 × 107, mutation and crossover rates of 0.02 and 0.8 respectively, the total runs number was set to 300. The NAA structure for docking was preoptimized by using quantum chemistry RHF/6-31G* calculations. The docking system consisted of protein dimer, zinc cations in active sites, and NAA. The Gastaiger charges were applied. Polar hydrogen atoms were assigned by Open Babel algorithm.27 The rotatable bonds in NAA were assigned by Autodock4.2 in an automatic mode assuming 4 torsional degrees of freedom.

cycle of hAsp consistent with the available structural and kinetics experimental data. Molecular dynamics simulations10 showed that the enzyme active site could be reached by a substrate molecule through the gates formed by several pairs of charged species, in particular, by the Arg71 and Glu293 side chains. Release of the products also occurred through these gates. Protein Data Bank contains an entry PDB ID: 2Q517 for an ensemble of 16 homodimer models of apo-hAsp, which suggests motion of the loops forming gates to the active site. Figure S1 in the Supporting Information illustrates this feature focusing on conformations of the gate-forming side chains Arg71 and Glu293. Their flexibility is an important issue when studying substrate binding to the protein and enzyme-catalyzed hydrolysis as discussed below. According to the structure PDB ID: 2I3C,8 also of the apo-form of hAsp, the dimeric interface, which is clearly shown in Figure 1, buries ≈1200 Å2 of surface accessible solvent area (SASA) and involves 12 hydrogen bonds and 2 salt bridges. We describe below the results of molecular modeling, including allosteric site search, docking and molecular dynamics (MD) simulations, which allow us to identify an area of substrate binding at the interface between the subunits of hAsp dimer and select a binding site with the highest affinity to NAA. These results explain the observed self-inhibition kinetics for NAA hydrolysis. We also apply dynamical network analysis to study signaling pathways of subunit communication to the remote sites with and without the substrate bound to the allosteric site. To perform this study we were encouraged by successful applications of dynamics correlation network analysis, in particular, those explaining allosteric inhibitor binding, known from the literature.11−17

■

METHODS We apply a variety of modern modeling approaches18−20 to characterize NAA allosteric targets in hAsp. An automatic tool 2000

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling Initial models for MD simulations were constructed by motifs of the crystal structure of native hAsp from human brain (PDB ID: 2O536). This structure was solved at better resolution than that in the PDB ID: 2I3C8 entry. All hydrogen atoms were added by using the Reduce package.28 Model macromolecule systems were solvated in a box of TIP3P water molecules. 6 and 8 sodium ions were added in apo-hAsp, and holo-hAsp models, respectively, to neutralize the systems. The water shells were built at a distance of 15 Å from the surface protein atoms and then equilibrated in 2 ns MD run with all protein atoms and Zn2+ ions frozen. Next, the protein was released and the entire system exceeding 72 000 atoms was equilibrated for 5 ns with subsequent energy minimization by conjugate-gradient algorithm in 5000 steps. Classical MD trajectories were computed using NAMD2.10 29 with CHARMM36 all atom force field.30,31 Each simulation was performed in NPT-ensemble at 300 K and 1 atm. The temperature control was implemented by means of Langevin dynamics and the pressure was maintained by using the Nosé− Hoover Langevin Piston method. For electrostatic, the ParticleMesh-Ewald method was implied. The long-ranged van der Waals interactions were given a cutoff distance of 13.5 Å and a switching distance at 10 Å. The integration step was set to 1 fs. We simulated several series of productive MD runs of 50−100 ns each for apo-hAsp and holo-hAsp model systems. Accelerated MD32−34 was applied to trace interaction of NAA with protein surface. Ten NAA molecules were placed randomly not less than 10 Å apart from hAsp surface. Twentysix sodium ions were added to neutralize the system. The water cell was constructed at a 15 Å distance from the atoms on the protein surfac. Both active site zinc cations were included in the model. After equilibration, 100 ns trajectory of accelerated MD was simulated. Analysis of MD trajectories was elaborated with the tools implemented in the VMD package (version 1.9.2).35 Root mean square deviations (RMSD) were calculated for equilibrated trajectories of 100 ns with a step 0.1 ns. The RMSD values reported in the text were computed for protein heavy atoms. H-bond occupancies were evaluated with the cutoff distance set as 3 Å and the donor-H-acceptor angles less than 20°. For salt bridges distance analysis, interaction were considered to be formed if the distance between any of the oxygen atoms of acidic residues and the nitrogen atoms of basic residues was within 3.2 Å. To analyze coupled motions and possible transition shifts in hAsp upon NAA binding to allosteric pocket site, the dynamical network analysis36 was employed. The method enables one to elaborate allosteric communication pathways from productive trajectories obtained in MD simulations. The method implies calculating of covariance and correlation matrices for timedependent variables defined as nodes of the network. These nodes serve to model a configuration transition of a system in time. In the present study, an α-carbon nodal method was implemented for its confirmed efficiency and reasonable computational time. Calculation of covariance and correlation values for all predefined nodes was performed with Carma 0.8 program35,36 for the last 20 ns of equilibrated MD trajectories of apo- and holo-hAsp model systems. The obtained values were then incorporated into a function that modeled signaling between nodes as a distance. The interaction between nodes is characterized as a length in a network space (d) obtained by functionalizing the correlation value

dij = −log(|Cij|) where Cij =

⟨Δri(t ) ·Δ rj(t )⟩ (⟨Δri(t )2 ⟩⟨Δ rj(t )2 ⟩)1/2

Δri(t ) = ri(t ) − ⟨ri(t )⟩

and ri (t) is the position of C-alpha atom of ith node. The network edges were defined as distances between pairs of nodes whose residues are within a cutoff (4.5 Å) for at least 75% of MD trajectory. Using weighted values of edges, it is possible to estimate putative signaling paths of a protein assuming the most probable the one that minimizes distances in a network space. The overall length of an optimal signaling path in a network space is defined as a sum of edges. Suboptimal paths characterized by accessible user-defined deviation from an optimal one could be evaluated in the same way. Suboptimal paths search was performed by subopt script provided with Dynamical Network Analysis tutorial files.37 The most highly correlated nodes can be organized into aggregates that in our case move in concert with each other. The signaling paths between these community aggregates were shown to correspond to biologically relevant ways of signaling. Information is transferred via particular node pairs characterized by high degree of network degeneracy, the so-called critical nodes. The community analysis was performed with Girvan−Newman algorithm.38 The clustering analysis was done by using gncommunities script also provided with VMD tutorial files.

■

RESULTS AND DISCUSSION Allocation of Putative Allosteric Sites. Application of Allosite allowed us to detect several substrate binding pockets. First, the program successfully found cavities near the active sites in both subunits (see Figure S2 in the Supporting Information) with Allosite score values 0.65 and 0.59, respectively. It should be noted that these scores do not pretend for quantitative characterization of active site availability, because the active sites of hAsp are located deeply inside the protein at the bottom of transport channel (∼12 Å from the surface), and they cannot be compared to the sites located at the surface. Beyond these sites, the program predicted three prospective substrate binding pockets (denoted here AS1, AS2, AS3) located at the dimer interface as shown in Figure 2. These pockets detected as feasibly druggable are located at the hAsp dimeric interface emphasizing the role of intersubunit area in hAsp activity regulation. Partly overlapping, AS1 and AS2 together occupy almost the whole cavity along the interface, whereas AS3 is between them. Practically, the whole area of the interface space is covered by these three sites. Respective scoring values for the detected pockets are listed in Table 1. The Allosite and Pocket feature scores do not allow us to discriminate the putative allosteric sites because their parameters are close enough. However, we can exclude from further consideration AS3 because its solvent accessible area (297 Å2) is less than that of the substrate (360 Å2). Application of the SiteMap algorithm to search allosteric sites in the hAsp dimer led to the results consistent with those 2001

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

discussed, the monomers are not fully symmetrical). The anchor and driven atoms identified as described in refs 39 and 40 are shown in Figure S3 in the Supporting Information. We note that this conformation of the binding site is well separated from other structures with less stable pockets. The pose with the closest score was significantly less populated (11% versus 78% in docking histogram) and had free energy value −2.6 kcal/mol. Therefore, we restricted our analysis by single allosteric site with the largest binding energy (−4.7 kcal/ mol). Molecular Dynamics Simulations. We compare results of MD simulations for the following model systems: hAsp monomers, hAsp dimer (ASPA-ASPB) in the absence of substrate at the allosteric site (apo-hAsp), hAsp dimer with the substrate molecule in the active site (hAsp ES), and hAsp dimer with the substrate attached at the allosteric site shown in Figure 3 (holo-hAsp). We pay primary attention to dynamics of the protein loops containing amino acid residues, which control entrance for NAA to the Zn-containing active site (Figure 4). The salt-bridges between positively charged Arg and Lys and negatively charged Glu and Asp can be formed and broken in protein dynamics; correspondingly, the gates to the active site can be either closed or open. Analysis of MD trajectories (5 runs by 100 ns) reveals different dynamical behavior of the gate-forming loops for the monomer and dimer protein molecules. We present in the Supporting Information graphs (Figure S4) showing a residue− residue distance fluctuation analysis of the key salt-bridges, Arg71-Glu293, Glu24-Arg63, Arg63-Glu285, and their H-bond occupancies, which suggest that the gates are likely to be closed in monomer species, whereas dynamics of the same pairs assumes that the gates are open in the dimer. Here, we point out that although in the crystals both monomers of the dimer are nearly identical (RMSD for the heavy atoms is 1.039 Å according to PDB ID: 2O536), their dynamical properties differ considerably. Only the first 20 ns of MD trajectories describing the dimer solvated in a water box demonstrate fairly similar dynamics in the subunits. Then, a shift in relative position of the subunits is observed, and arrangement of the residues at the interface is not symmetrical. This behavior explains why the allosteric site is found at the interface only for one subunit, here, ASPB. In our case, fluctuations of gate forming loops in the ASPB subunit resemble those in a separated monomer molecule, while in ASPA tendencies to gate opening are observed. MD trajectories evaluated for the complex (holo-hAsp) revealed that upon NAA binding to the allosteric pocket a mobility of the loops 62−74 and 282−294, which control access to the active site in ASPA, decreases substantially. This is illustrated by average RMSD values computed for five parallel 50 ns trajectories as shown in the left panels in Figure 5. The RMSD values for the 62−74 and 282−294 loops in ASPA in the holo-hAsp dimer are about twice as low as the values obtained for the apoprotein. Mobility of the loops is slowly restored after approximately 50 ns which is an averaged time when NAA remains bound to the allosteric site. The latter value was obtained from 5 parallel 100 ns runs. Upon NAA binding to the allosteric site, a substantial reorganization of crucial H-bonds and salt-bridges that define the dimer hAsp structure occurs. Also, the number of stable nonbonded interactions between hAsp subunits increases. The reorganization of hydrogen bond network confirmed by Hbond occupancy (Figure 6) and salt-bridges lengths analysis of

Figure 2. Three allosteric pockets (distinguished by different colors) predicted by Allosite. Note that the cumulative picture in the right panel shows another view on the dimer than in Figure 1 and in the left panels in this figure.

Table 1. Results of Allosite Scoring Analysis for Allosteric Binding Pockets volume (Å3) SASA (Å2) feature score perturbation score Allosite score

AS1

AS2

AS3

1443 853 0.76 0.85 0.78

1182 515 0.69 1.00 0.76

443 297 0.88 0.35 0.77

obtained by Allosite. Table 2 summarizes properties of the pockets with the highest score. Table 2. Results of SiteMap Allosteric Site Search: Scoring Functions and Allosteric Binding Properties for Detected Pockets site score size (site points) D score volume (Å3) exposure enclosure contact hydrophobicity hydrophilicity

SM1

SM2

SM3

1.015 407 1.037 1291 0.569 0.713 0.874 0.605 1.016

0.986 80 0.930 203 0.503 0.793 1.117 0.910 1.186

0.702 49 0.674 113 0.743 0.566 0.652 0.206 0.992

The pocket SM3 with the lowest score was excluded from consideration. In agreement with the results obtained by using Allosite, the SM1 pocket with the best values of scoring functions occupies a large fraction of the hAsp dimer interface partly overlapping with SM2 and AS1. A view on SM1 area is shown in Figure 3. Application of the docking procedure with Autodock4.2 allowed us to find a binding site of the substrate and to estimate the corresponding binding energy. The NAA molecule was bound to the allosteric binding site located on the top of the area detected by SiteMap. The overall distribution of this NAA conformation resulted in 80% of the samples obtained within the same grid. Orientation of NAA was selected as the position with the highest affinity score. The computed binding free energy is −4.7 kcal/mol. Inset in Figure 3 shows that in this environment, NAA interacts with the Lys292, Arg233 side chains and Glu290, Gly237 main chains of ASPB monomer (as 2002

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Figure 3. SiteMap area corresponding to the SM1 site (see Table 2) and the allosteric binding pocket of NAA found by Autodock at the ASPB side. We use another orientation than that in Figure 2 for better visibility. Inset shows a closer view on the NAA binding pocket at the dimer interface.

Figure 4. Views from different perspectives on the residues controlling entrance for a substrate molecule to the active site in each monomer subunit.

formation of the ES complex, thus closing the gate behind the substrate. The gates are supposed to be open after release of products: crystal structures6−8 of hAsp in the apo form clearly demonstrate that the side chain of Arg71 is oriented in the solute opening the gates to the active site. Also, consistently with kinetics studies, the results of simulation of reaction mechanism at the active site10 show that the rate limiting step is assigned to product release, but not to cleavage of the chemical bond in the substrate. Therefore, short presence of NAA at the active sites should not interfere with the effects of closing and opening gates to the active sites. Dynamical Network Analysis. Further studies of allosteric effect on intersubunit interactions in the hAsp dimer were carried out with the dynamical network analysis. We selected the last 20 ns of equilibrated dimer trajectories out of the total 100 ns length to search for communities in the hAsp dimer structure both in the apo- and holo-forms. We found that the residues responsible for interface interactions were distributed among three dynamical communities as distinguished by different colors in Figure 7. We note that the collective motions within communities should have an impact on the hydrolysis reaction; in particular, communities I and II share dynamics with the active site residues Glu24, Glu178, His116, which directly participate in catalytic activity of hAsp.10 Figure 7 illustrates that substrate binding to the allosteric site (indicated by asterisk in the right panel) affects the protein dynamical properties at the dimer interface as reflected by

the gate residues, caused by substrate binding to the allosteric site, can result in poor access of the enzyme active site. From these MD simulations, we conclude that NAA binding to the allosteric binding pocket shown in Figure 3 stimulates changes in interface interactions between two hAsp subunits, which causes increased rigidity of the loops with the gateforming residues; therefore, access of the substrate molecules to the active sites for the chemical reaction of NAA hydrolysis is obstructed. The results of accelerated molecular dynamics simulations for hAsp dimer with several substrate molecules distributed around the protein in solution show that any of the NAA molecules is periodically bound to the allosteric site. Therefore, we expect that at high NAA concentrations, the allosteric site at the ASPB subunit is well populated. A possible effect of presence of NAA bound at the active site on dynamical properties of the model systems was studied by comparing MD trajectories for the enzyme−substrate complex of hAsp dimer (hAsp ES) with and without NAA at the remote site. The initial coordinates of the ES complex for these simulations were taken from the previous results of quantum mechanics/molecular mechanics calculations.10 Figure S5 in the Supporting Information shows that the gates to the active site remain closed when NAA is trapped in the ES complex, independent of presence or absence of NAA at the remote site. This observation is in accord with the atomic-scale mechanism of NAA hydrolysis,10 which shows involvement of Arg71 in the 2003

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Figure 5. RMSD values for the loops with the gate-controlling residues for apo-hAsp (gray) and holo-hAsp (pink). Left panels, subunit ASPA; right panels, subunit ASPB.

Figure 6. H-bond occupancies for intersubunit interaction pairs in 50 ns MD simulation.

gate-forming residues Arg71, Glu293 observed in the crystal structure of the dimer (see Figure S1 in the Supporting Information). In holo-Asp (right side in Figure 7), analysis of the dynamical network reveals substantially different community fragmentation as compared to the apoprotein form. This time, only two communities completely within ASPA monomer (I and II on the right side in Figure 7) contain the gate-forming residues. We present in the Supporting Information (Figure S8) the graph illustrating the Kolmogorov−Smirnov test performed to compute p-values of RMSF for the C-alpha

network communities. In apo-hAsp (left side in Figure 7), the gate-forming residues shown in Figure 4 (Arg63, Arg71, Lys228, Lys291, Glu293, Glu285, Asp68) are spread among different dynamical communities (I, II, III on the left side in Figure 7) extending over both monomers. Such property is an indication of small dynamical correlation in the motion of the gate-forming species. This conclusion from the dynamical network analysis is consistent, first, with somewhat disordered motion of the gate forming loops of apo-hAsp dimer noticed in MD trajectories (Figure 5), second, with high flexibility of the 2004

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Figure 7. Dynamical communities containing the intersubunit interaction pairs in the hAsp dimer in the apo-form (left) and in the holo-form (right). Red asterisk in the holo-form indicates an area of the allosteric binding pocket.

Figure 8. Critical nodes with the gate-forming Arg63 in apo-hAsp (left) and hopo-hAsp (right).

both apo- and hopo-forms). The community analysis indicates that only two critical nodes (Arg196(A)−Lys200(A), Asn23(B)−Gln184(B)) are preserved upon NAA binding to the allosteric site. Therefore, location of critical nodes connecting ASPA and ASPB units in holo-hAsp indicates a rearranged dynamical network formed upon NAA binding. The latter observation is consistent with the distribution of inter subunit H-bond occupancies for apo- and holo-structures shown in Figure 6. Therefore, all modeling results support conclusion of an enhanced effect of interface noncovalent bonds on mobility of the transport channel for the substrate in ASPA subunit. Importantly, NAA shares one small community (community III in Figure 7) at the interface area. All allosteric site forming residues belong to the same community, but the edge with the highest weight is formed between NAA and Tyr289 (Bsubunit) nodes. Therefore, it is instructive to analyze signaling pathways connecting the Tyr289(B) node with the critical residues in the A-subunit. An example of such a signaling pathway to Arg71(A) is shown in panel a in Figure 9. Panels b and c illustrate that lengths of these signaling pathways decrease, but the total number of pathways increases upon NAA binding to the allosteric site. This means, in turn, that the coupled residue−residue signaling motions involved in formation of the enzyme−substrate complex are stronger when NAA is bound to the allosteric site. Namely, substrate binding to the allosteric site at the surface of ASPB restricts the motion of the gate residues in monomer ASPA, which blocks

atoms in apo-hAsp and holo-Asp along the 30 ns trajectory. These results support a conclusion that motion of the protein main chain is considerably influenced by allosteric interactions. It is instructive to characterize intercommunity interactions in dynamical networks by critical nodes; the latter indicate high signal degeneracy between community aggregates, i.e., the protein regions with high regional correlations. Positions of critical nodes point out on possible sites for signal disruption. Importantly, critical nodes responsible for intersubunit interactions between ASPA and ASPB in holo-hAsp differ from those in apo-hAsp. Figures S6, S7 in the Supporting Information illustrate differences in community fragmentation as well as in distribution of critical nodes in apo-hAsp and holoAsp. Comparison of fragments of the dimer interface with an assignment of the gate-forming residues in the ASPA subunit in apo-hAsp (Figure S6) and holo-hAsp (Figure.S7) shows an increased regularity in interface dynamics caused by allosteric binding of NAA. Figure 8 illustrates a tendency in loop motion taking the nodes with the gate-forming residue Arg63 as an example. The loops 62−74 and 282−294 are pulling apart opening the active site in the apo-form (left panel in Figure 8), whereas motion of these loops in the holo-form (right panel in Figure 8) tends to block an access to the active site. The total amount (12) of critical nodes located in the ASPA unit in holo-Asp increased as compared to that (6) in apohAsp, whereas in the ASPB unit there were no changes (7 in 2005

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Figure 9. (a) Specific signaling pathway from Tyr289 (ASPB-subunit) to Arg71 (ASPA-subunit). (b) Pathway lengths from Tyr289 (ASPB) to the gate-forming residues in ASPA for apo-hAsp (gray) and holo-hAsp (magenta). (c) Total number of signal pathways leading from Tyr289 (ASPB) to the gate-forming and active site residue Glu24 in ASPA for apo-hAsp (gray) and holo-hAsp (magenta).

approaches, Allosite and SiteMap, leads to a conclusion that the allosteric binding sites with the high affinity to NAA are located at the hAsp dimer interface. Molecular docking through the specified areas points out a specific binding pocket located at the dimeric interface. The Autodock algorithm based on a validated semiempirical force field predicts the binding free energy at this pocket −4.7 kcal/mol. The results of further molecular dynamics simulations show that binding of NAA at the specified allosteric site leads to reduced flexibility of the protein loops with the amino acid residues forming the gates to the enzyme active site. Consequently, access of the substrate to the catalytic sites is blocked upon allosteric binding of NAA at the dimeric interface. According to our results, substrate binding to the active site of the ASPA monomer is influenced by binding to the identified allosteric site at ASPB. Dynamical network analysis allows us to uncover how the residues forming gates to the enzyme active sites recognize perturbation caused by allosteric substrate binding; we can find the corresponding signaling pathways through the protein. This analysis also shows that NAA binding to the allosteric site shortens signal protein pathways and modifies the signaling route. Significance of studies of the hAsp-catalyzed hydrolysis of NAA is explained by a request to treat severe diseases caused by enzyme dysfunction due to genetic disorder. In particular, the mutation Glu285Ala is one of the most common in Canavan disease patients. According to our simulations with the dynamical network analysis, the residue at position 285 on the loop 282−294 occurs in the immediate vicinity of the

substrate binding to the catalytic site; however, perturbation at the allosteric site is transmitted via increasing number of signaling paths. Finally, we note occurrence of the known clinical mutations in hAsp at positions 231, 285, 288, and 295 on the revealed signaling pathways connecting the specified allosteric site and the gate-forming residues in the ASPA monomer. One of the most significant mutation in Canavan disease patients, Glu285Ala, is connected with a replacement of residue 285, which occurs in approximately 10% of signaling pathways. The residues Tyr231 and Tyr288, also important for Canavan disease patients, occur in approximately 20% of signaling pathways.

■

CONCLUSIONS All modeling approaches applied in this work consistently predict a molecular mechanism of allosteric inhibition in the NAA hydrolysis by hAsp. First, we note an important contribution from molecular dynamics simulations showing that consideration of the protein dimer is a necessary step in modeling. Dynamics of the monomeric species demonstrate that conformations with the closed gates to the enzyme active site are dominating, whereas dynamics of the dimer is consistent with availability of the entrance gate for substrate molecules. High flexibility of the most important gate forming salt bridge Arg71-Glu293 agrees with the results of crystallography studies (PDB ID: 2Q517) showing variety of conformations of the corresponding loops in the dimer. Next, the role of dimeric interface in allosteric modulation of the reaction is revealed. Application of two modeling 2006

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling

Integrative View of Protein, Gene and 3D Structural Information. Nucleic Acids Res. 2016, 45, 271−281. (10) Kots, E. D.; Khrenova, M. G.; Lushchekina, S. V.; Varfolomeev, S. D.; Grigorenko, B. L.; Nemukhin, A. V. Modeling the Complete Catalytic Cycle of Aspartoacylase. J. Phys. Chem. B 2016, 120, 4221− 4231. (11) Wang, W.; Jiang, C.; Zhang, J.; Ye, W.; Luo, R.; Chen, H.-F. Dynamics Correlation Network for Allosteric Switching of PreQ1 Riboswitch. Sci. Rep. 2016, 6, 31005. (12) Rahman, M.; Liu, H.; Wadood, A.; Chen, H.-F. Allosteric Mechanism of Cyclopropylindolobenzazepine inhibitors for HCV NS5B RdRp via Dynamic Correlation Network Analysis. Mol. BioSyst. 2016, 12, 3280−3293. (13) Xu, L.; Ye, W.; Jiang, C.; Yang, J.; Zhang, J.; Feng, Y.; Luo, R.; Chen, H.-F. Recognition Mechanism between Lac Repressor and DNA with Correlation Network Analysis. J. Phys. Chem. B 2015, 119, 2844−2856. (14) Yang, J.; Liu, H.; Liu, X.; Gu, C.; Luo, R.; Chen, H.-F. Synergistic Allosteric Mechanism of Fructose-1,6-bisphosphate and Serine for Pyruvate Kinase M2 via Dynamics Fluctuation Network Analysis. J. Chem. Inf. Model. 2016, 56, 1184−1192. (15) Zhang, J.; Luo, H.; Liu, H.; Ye, W.; Luo, R.; Chen, H.-F. Synergistic Modification Induced Specific Recognition between Histone and TRIM24 via Fluctuation Correlation Network Analysis. Sci. Rep. 2016, 6, 24587. (16) Tse, A.; Verkhivker, G. M. Small-World Networks of Residue Interactions in the Abl Kinase Complexes with Cancer Drugs: Topology of Allosteric Communication Pathways Can Determine Drug Resistance Effects. Mol. BioSyst. 2015, 11, 2082−2095. (17) Del Sol, A.; Tsai, C. J.; Nussinov, R.; Ma, B. The Origin of Allosteric Functional Modulation: Multiple Pre-Existing Pathways. Structure 2009, 17, 1042−1050. (18) Nussinov, R.; Tsai, C. J. Allostery in Disease and in Drug Discovery. Cell 2013, 153, 293−305. (19) Sgrignani, J.; Bon, M.; Colombo, G.; Magistrato, A. Computational Approaches Elucidate the Allosteric Mechanism of Human Aromatase Inhibition: A Novel Possible Route to Small-Molecule Regulation of CYP450s Activities? J. Chem. Inf. Model. 2014, 54, 2856−2868. (20) Wagner, J. R.; Lee, C. T.; Durrant, J. D.; Malmstrom, R. D.; Feher, V. A.; Amaro, R. E. Emerging Computational Methods for the Rational Discovery of Allosteric Drugs. Chem. Rev. 2016, 116, 6370− 6390. (21) Huang, W.; Lu, S.; Huang, Z.; Liu, X.; Mou, L.; Luo, Y.; Zhao, Y.; Liu, Y.; Chen, Z.; Hou, T.; Zhang, J. Allosite: a Method for Predicting Allosteric Sites. Bioinformatics 2013, 29, 2357−2359. (22) Shen, Q.; Wang, G.; Li, S.; Liu, X.; Chen, Z.; Song, K.; Yan, J.; Geng, L.; Huang, Z.; Huang, W.; Chen, G.; Zhang, J.; Lu, S. ASD v3.0: Unraveling Allosteric Regulation with Structural Mechanisms and Biological Networks. Nucleic Acids Res. 2016, 44, 527−535. (23) Le Guilloux, V.; Schmidtke, P.; Tuffery, P. Fpocket: An Open Source Platform for Ligand Pocket Detection. BMC Bioinf. 2009, 10, 168. (24) Halgren, T. A. Identifying and Characterizing Binding Sites and Assessing Druggability. J. Chem. Inf. Model. 2009, 49, 377−389. (25) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. Autodock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 2009, 30, 2785−2791. (26) Solis, F. J.; Wets, R. J-B. Minimization by Random Search Techniques. Math. Oper. Res. 1981, 6, 19−30. (27) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An Open Chemical Toolbox. J. Cheminf. 2011, 3, 33. (28) Word, J. M.; Lovell, S. C.; Richardson, J. S.; Richardson, D. C. Asparagine and Glutamine: Using Hydrogen Atom Contacts in the Choice of Side-Chain Amide Orientation. J. Mol. Biol. 1999, 285, 1735−1747.

critical nodes important for signal pathways; therefore, replacement of Glu285 affects allosteric regulation of the NAA hydrolysis.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00133. Additional details on crystal structure analysis, docking analysis, molecular dynamics simulations, and dynamical network analysis (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Prof. Alexander V. Nemukhin. Phone: +7-495-939-10-96. Email: [email protected]; [email protected]. ORCID

Alexander V. Nemukhin: 0000-0002-4992-6029 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. M. Khrenova for advice and helpful discussion. This study was supported by the Russian Science Foundation (project # 14-13-00124). We acknowledge the use of supercomputer resources of the Lomonosov Moscow State University41 and of the Joint Supercomputer Center of the Russian Academy of Sciences.

■

REFERENCES

(1) Bjartmar, C.; Kidd, G.; Mork, S.; Rudick, R.; Trapp, B. D. Neurological Disability Correlates with Spinal Cord Axonal Loss and Reduced N-Acetyl Aspartate in Chronic Multiple Sclerosis Patients. Ann. Neurol. 2000, 48, 893−901. (2) Tsai, G.; Coyle, J. T. N-Acetylaspartate in Neuropsychiatric Disorders. Prog. Neurobiol. 1995, 46, 531−540. (3) Vermathen, P.; Laxer, K. D.; Schuff, N.; Matson, G. B.; Weiner, M. W. Evidence of Neuronal Injury Outside the Medial Temporal Lobe in Temporal Lobe Epilepsy: N-Acetylaspartate Concentration Reductions Detected with Multisection Proton MR Spectroscopic Imaging − Initial Experience. Radiology 2003, 226, 195−202. (4) Matalon, R.; Michals, K.; Sebesta, M.; Deanching, M.; Gashkoff, P.; Casanova, J.; Optiz, J. M.; Reynolds, J. F. Aspartoacylase Deficiency and N-acetylaspartic Acid Uria in Patients with Canavan Disease. Am. J. Med. Genet. 1988, 29, 463−471. (5) Le Coq, J.; An, H.-J.; Lebrilla, C.; Viola, R. E. Characterization of Human Aspartoacylase: the Brain Enzyme Responsible for Canavan Disease. Biochemistry 2006, 45, 5878−5884. (6) Le Coq, J.; Pavlovsky, A.; Malik, R.; Sanishvili, R.; Xu, C.; Viola, R. E. Examination of the Mechanism of Human Brain Aspartoacylase through the Binding of an Intermediate Analogue. Biochemistry 2008, 47, 3484−3492. (7) Levin, E. J.; Kondrashov, D. A.; Wesenberg, G. E.; Phillips, G. N., Jr. Ensemble Refinement of Protein Crystal Structures: Validation and Application. Structure 2007, 15, 1040−1052. (8) Bitto, E.; Bingman, C. A.; Wesenberg, G. E.; McCoy, J. G.; Phillips, G. N. Structure of Aspartoacylase, the Brain Enzyme Impaired in Canavan Disease. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 456−461. (9) Rose, P. W.; Prli, A.; Altunkaya, A.; Bi, C.; Bradley, A. R.; Christie, C. H.; Di Costanzo, L.; Duarte, J. M.; Dutta, S.; Feng, Z.; Green, R. C.; Goodshell, D. S.; Hudson, B.; Karlo, T.; Lowe, R.; Peisach, E.; Randle, C.; Rose, A. S.; Shao, C.; Tao, Y. P.; Valasatava, Y.; Voigt, M.; Westbrook, J. D.; Woo, J.; Yang, H.; Young, J. Y.; Zardecki, C.; Berman, H. M.; Burley, S. K. The RCSB Protein Data Bank: 2007

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008

Article

Journal of Chemical Information and Modeling (29) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (30) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D., Jr. CHARMM General Force Field (CGenFF): A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671−690. (31) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell, A. D., Jr. Extension of the CHARMM General Force Field to SulfonylContaining Compounds and Its Utility in Biomolecular Simulations. J. Comput. Chem. 2012, 33, 2451−2468. (32) Hamelberg, D.; Mongan, J.; McCammon, A. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919−11929. (33) Wang, Y.; Harrison, C. B.; Schulten, K.; McCammon, J. A. Implementation of Accelerated Molecular Dynamics in NAMD. Comput. Sci. Discovery 2011, 4, 015002. (34) De Oliveira, C. A. F.; Hamelberg, D.; McCammon, J. A. Estimating Kinetics Rates from Accelerated Molecular Dynamics Simulations: Alanine Dipeptide in Explicit Solvent as a Case Study. J. Chem. Phys. 2007, 127, 175105. (35) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (36) Sethi, A.; Eargle, J.; Black, A. A.; Luthey-Schulten, Z. Dynamical Networks in tRNA: Protein Complexes. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 6620−6625. (37) Eargle, J.; Li, L.; Luthey-Schulten, Z. Dynamical Network Analysis, 2012, University of Illinois at Urbana−Champaign, LutheySchulten Group, NIH Resource for Macromolecular Modeling and Bioinformatics, http://www.scs.illinois.edu/schulten/tutorials/ network/ (last accessed May17, 2017). (38) Girvan, M.; Newman, M. Community Structure in Social and Biological Networks. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 7821− 7826. (39) Nussinov, R.; Tsai, C.-J.; Liu, J. Principles of Allosteric Interactions in Cell Signaling. J. Am. Chem. Soc. 2014, 136, 17692− 17701. (40) Nussinov, R.; Tsai, C.-J. Unraveling Structural Mechanisms of Allosteric Drug Action. Trends Pharmacol. Sci. 2014, 35, 256−264. (41) Voevodin, Vl.V.; Zhumatiy, S. A.; Sobolev, S. I.; Antonov, A. S.; Bryzgalov, P. A.; Nikitenko, D. A.; Stefanov, K. S.; Voevodin, Vad. V. Practice of ″Lomonosov″ Supercomputer. Open Systems J. (Mosc.) 2012, 7, 36−39.

2008

DOI: 10.1021/acs.jcim.7b00133 J. Chem. Inf. Model. 2017, 57, 1999−2008