Allostery through the computational microscope: cAMP activation of a canonical signalling domain

Abstract

Ligand-induced protein allostery plays a central role in modulating cellular signalling pathways. Here using the conserved cyclic nucleotide-binding domain of protein kinase A's (PKA) regulatory subunit as a prototype signalling unit, we combine long-timescale, all-atom molecular dynamics simulations with Markov state models to elucidate the conformational ensembles of PKA's cyclic nucleotide-binding domain A for the cAMP-free (apo) and cAMP-bound states. We find that both systems exhibit shallow free-energy landscapes that link functional states through multiple transition pathways. This observation suggests conformational selection as the general mechanism of allostery in this canonical signalling domain. Further, we expose the propagation of the allosteric signal through key structural motifs in the cyclic nucleotide-binding domain and explore the role of kinetics in its function. Our approach integrates disparate lines of experimental data into one cohesive framework to understand structure, dynamics and function in complex biological systems.

ARTICLE

Received 24 Aug 2014 | Accepted 22 May 2015 | Published 6 Jul 2015

Ligand-induced protein allostery plays a central role in modulating cellular signalling pathways. Here using the conserved cyclic nucleotide-binding domain of protein kinase As (PKA) regulatory subunit as a prototype signalling unit, we combine long-timescale, all-atom molecular dynamics simulations with Markov state models to elucidate the conformational ensembles of PKAs cyclic nucleotide-binding domain A for the cAMP-free (apo) and cAMP-bound states. We nd that both systems exhibit shallow free-energy landscapes that link functional states through multiple transition pathways. This observation suggests conformational selection as the general mechanism of allostery in this canonical signalling domain. Further, we expose the propagation of the allosteric signal through key structural motifs in the cyclic nucleotide-binding domain and explore the role of kinetics in its function. Our approach integrates disparate lines of experimental data into one cohesive framework to understand structure, dynamics and function in complex biological systems.

DOI: 10.1038/ncomms8588

Allostery through the computational microscope: cAMP activation of a canonical signalling domain

Robert D. Malmstrom1,2, Alexandr P. Kornev3, Susan S. Taylor1,3 & Rommie E. Amaro1,2

1 Department of Chemistry and Biochemistry, University of California, San Diego, 9500 Gilman Drive, La Jolla, California 92093-0340, USA. 2 National Biomedical Computation Resource, University of California, San Diego, 9500 Gilman Drive, La Jolla, California 92093-0340, USA. 3 Department of Pharmacology, University of California, San Diego, 9500 Gilman Drive, La Jolla, California 92093-0340, USA. Correspondence and requests for materials should be addressed to R.E.A. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588

Since the introduction of the allosteric effect in L-threonine deaminase1,2, researchers have sought to understand the mechanisms of allostery. The classical phenomenological

models proposed by MonodWymanChangeux3 and (Pauling)KoshlandNemethyFilmer4,5 have been extended over time to emphasize population shifts or modulations in the conformational ensemble of the proteins6,7. Fundamentally, it is the character of the underlying free-energy landscape that determines the mechanism of allostery611.

While experimental and computational approaches provide insight into allosteric mechanisms12,13, a robust atomic-level description of the free-energy landscape remains a challenge due to a variety of methodological limitations. Without an atomic-level description of the conformational ensemble and its corresponding free-energy landscape, the details of allostery can remain hidden within ensemble averages or constrained perspectives of protein conformation and dynamics.

The objective of our work is to apply cutting-edge computational approaches to explore the conformational ensemble of a protein with and without a bound ligand, thereby gaining insight into how the ensemble gives rise to protein function. This work focuses on the mechanism of ligand-driven allostery in the cyclic nucleotide-binding domain of protein kinase As (PKA) regulatory subunit. Despite extensive research conducted on this system14, questions still remain.

The intercellular activation of PKA by cAMP is a prototypical example of ligandprotein allostery, which occurs through the cooperative binding of cAMP to tandem cyclic nucleotide-binding domains, designated A and B, in each of PKAs regulatory subunits14. Crystallographic data show that the regulatory subunit, isoform Ia, undergoes signicant conformational changes on binding cAMP. It transitions from an extended holoenzyme conformation that inhibits PKAs catalytic subunits to a compact cAMP-bound conformation that releases and activates PKAs catalytic subunits15,16. Within the cyclic nucleotide-binding domain, the inactive holoenzyme and the active cAMP-bound conformations of the regulatory subunit are characterized by the conformational changes in three structural motifsthe N3A motif, the phosphate-binding cassette (PBC) and the B/C helix17,18, leaving the core b-barrel subdomain largely unchanged (Fig. 1). These structural motifs represent the fundamental signal transduction component and form the binding interface between the cyclic nucleotide-binding domain and PKAs catalytic subunit.

The existing crystallographic data provide only two conformational states for the cyclic nucleotide-binding domains, the holoenzyme (a.k.a. H or inactive)16,19 and the cAMP-bound conformation (a.k.a. B or active)15. Nuclear magnetic resonance (NMR) dynamics data of the minimal functional component of PKAs regulatory subunit, the A cyclic nucleotide-binding

domain (a.a. 119244), indicate the presence of two dominant conformational states2022 and suggests conformational selection as the mechanism of allostery23. On cAMP binding, changes in NMR chemical shift data implicate the PBC as a major dynamic element, while few changes are observed in the N3A motif or the B/C helix.

However, such experiments are limited in their ability to elucidate details of the full structural ensemble. Instead, they provide averaged metrics that correlate to dynamic motions. Achieving an atomistic description of the underlying conformational ensemble and resolving the corresponding free-energy landscape will clarify the mechanisms by which cAMP modulates the function of the cyclic nucleotide-binding domain.

Because the PKA holoenzyme contains four cyclic nucleotide-binding domains, cAMP-induced activation is a complicated process. A recent study by Boras et al.24 provided a general blueprint of the events during activation at the subcellular level, but the molecular mechanism of activation remains to be elucidated.

To build a foundation on which to address these questions, we simulated fully solvated atomistic models of the conformational ensembles of the A cyclic nucleotide-binding domain (hereafter referred to as the cyclic nucleotide-binding domain or CBD), with and without cAMP bound. Extensive molecular dynamics (MD) simulations were integrated using Markov state model theory to explore PKAs conformational space and long-timescale dynamics. Markov state models depict the interaction dynamics of discrete interconnected states as a transition probability matrix, at a xed lag time, assuming that the transitions between the states are independent of previous transitions (that is, Markovian)2528. By assigning individual frames extracted from MD trajectories to discrete conformational states, sampling from many separate trajectories can be integrated into one coherent framework that captures the kinetics and thermodynamics of the conformational ensemble at atomic resolution.

Similar approaches have been used to study protein folding29, biased transition pathways between functional states30,31, identication of cryptic allosteric sites32 and ligand regulation of G-protein-coupled receptors33. Our work is unique in that it directly assesses the atomic conformational landscape using initial unbiased long-timescale MD simulations augmented with adaptive sampling of one protein in two functional states, as opposed to only combining many short-timescale simulations or sampling along a predetermined reaction coordinate.

We nd that the conformational ensembles of cAMP-free (CBDApo) and cAMP-bound (CBDcAMP) functional states exist within a shallow free-energy landscape that allows access to both experimentally determined functional conformations, indicating conformational selection as the principal mechanism of allostery with the isolated cyclic nucleotide-binding domain. We nd that the addition of cAMP modies the principal motions of the cyclic nucleotide-binding domain, which correspond to transitions between active and inactive states. Furthermore, our approach exposes the propagation of the allosteric signal through the cyclic nucleotide-binding domains signalling motifs.

Our results complement existing structural15,16 and dynamical2023,34 experimental data, and extend our understanding of the mechanism of ligand-induced allostery within the cyclic nucleotide-binding domain. Further, they provide a general framework for understanding the function of the ancient and ubiquitous cyclic nucleotide-binding domain35,36 and allosteric mechanisms more generally.

ResultsModelling the conformational landscape. To explore and map the conformational landscapes of the cyclic nucleotide-binding

Inactive conformation Active conformation

N3A motif

B/C helix

cAMP

N3A motif

PBC

PBC

B/C helix

Figure 1 | Experimentally determined conformational changes in the cyclic nucleotide-binding domain. (RIa a.a. 119 to 244) in cAMP bound active, right and holoenzyme inactive, left, with key signalling structural subdomains highlighted: PBC, red; B/C helix, ice blue; and N3A motif, ochre.

2 NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588 ARTICLE

domain, we integrated multiple MD simulations with Markov state models. Explicit solvent MD simulations of the cyclic nucleotide-binding domain were performed with and without cAMP-bound systems (CBDcAMP and CBDApo, respectively).

The MD simulations of both systems sampled a totalled 74 ms (Supplementary Table 1). The simulations were composed of four long-timescale, B13-ms simulations on the Anton supercomputer37 (at the Pittsburgh Supercomputing Center) and multiple parallel short-timescale, 0.51 ms simulations using

GPU-accelerated AMBER12 (refs 3840). For both systems, sampling simulations were initiated from experimentally determined atomic coordinates of the cyclic nucleotide-binding domain15,16 in both the inactive and active conformations. Initial sampling was followed by multiple rounds of directed sampling (1015 ns simulations using GPU-accelerated AMBER12 (refs 3840)) to rene the Markov state models41. Directed sampling simulations were initiated from intermediate conformations selected near poorly sampled transitions within the Markov state model. For both the CBDcAMP and CBDApo systems, no single simulation sampled a complete transition between the experimentally determined inactive and active structures. However, simulations started from either experimental conformation overlapped in conformational space, based on alpha carbon root mean squared distance (r.m.s.d.), indicating that the MD trajectories could be integrated meaningfully into a single Markov state model (Supplementary Fig. 1).

To build the Markov state models, we needed to dene the conformational space, and then divide the space into discrete conformational states. As in protein-folding models29, the cyclic nucleotide-binding domain conformations are dened by the atomic coordinates of each residues alpha carbon. We divided the sampled conformational space into discrete conformational states through r.m.s.d. clustering using MSMBuilder2 (ref. 42). We prealigned the structures to a common reference frame to calculate r.m.s.d. to capture important translational motions lost using the standard minimum r.m.s.d. approach, which aligns each pair of compared structures to each other before determining r.m.s.d.. Inspection of the individual MD trajectories indicated that the b-barrel of the cyclic nucleotide-binding domain was stable relative to the motions of the N3A motif, the PBC and the B/C helix (Supplementary Fig. 2). Therefore, we used the b-barrel as our common reference frame.

All conformations sampled from the MD simulations of both the CBDApo and CBDcAMP systems were clustered together to create a comprehensive framework of conformational states. Once the conformation states were identied, Markov state models were built for each system. This means that the models were built using the same division of conformational space, but each model included only states sampled in its respective simulations. (A comparison between this approach and clustering conformations from CBDApo and CBDcAMP separately is discussed below.)

The Markov state models were built by specifying the maximum distance or cutoff for the clustering algorithm, which determined how the conformational space was divided, and a lag time for the model. Model parameters were selected because they were Markovian, as indicated by implied timescale plots, and could maximize the number of conformational states. In addition, consistency between the Markov state models and the MD simulations were determined by the ChapmanKolmogorov test26 (that is, were the models generated trajectories similar to the results from the original MD simulations?; Supplementary Figs 310). An evaluation of several r.m.s.d. cutoffs and lag times indicated that the best Markov state model utilized a clustering cutoff of 3.0 and a lag time of 9.6 ns (Supplementary Figs 3 and 4).

Graphs of the Markov state models for CBDApo and CBDcAMP are visualized as a network (Fig. 2). Each node corresponds to a cluster of similar conformations or a single conformational state. The diameters of the nodes are proportional to the log of equilibrium population of the conformational state. Thus, smaller nodes indicate conformations that are more rarely sampled. The location of each node is determined by the r.m.s.d. of the structure of the clusters generator, approximately the centroid of the cluster, to the inactive and active crystallographic structures. We note that the location of each node gives a general indication of the conformation of each state relative to the experimental structures, but node locations are a projection from a higher-dimensional space. (Separation of the experimentally determined active node, indicated in Fig. 2, from the bulk of the other nodes is an artifact of the clustering algorithm in MSMBuilder2 (ref. 42), which selects the rst MD frame (that is, the equilibrated starting active conformation) as a generator for the rst cluster. Notice that the separation of the active conformation node from its neighbours is B3 or the r.m.s.d.

cutoff of the clustering algorithm.)

Ligand-meditated changes in the conformational ensemble. A comparison of the models of the conformational state ensembles for the CBDApo and CBDcAMP systems allowed us to assess how the binding of cAMP modies the overall free-energy landscape. We compared the distribution of conformational states between the two systems, the overall character of the free-energy landscape and the kinetics of traversing the free-energy landscape between functional conformations states. For this comparison, and throughout this work, we assumed that the crystallographic structures represent the functional conformations (active or inactive states) of the cyclic nucleotide-binding domain. The nodes containing the functional conformations are indicated in Fig. 2.

Active conformation

10

2

RMSD from inactive conformation ()

8

6

4

Inactive conformation

0 2 4 6 8 10 12

RMSD from active conformation ()

Figure 2 | Overlay of Markov State Models State Graphs for the cyclic nucleotide-binding domain with (CBDcAMP, magenta) and without (CBDApo, cyan) cAMP bound. The position of each conformational state node is the r.m.s.d. of the generator of each cluster relative to the reference conformations. The diameter of each conformational state node (cyan, magenta) is proportional to the relative the log of the equilibrium population. (that is, the lager the node the more probable the state at equilibrium.) Cyan and magenta nodes are transparent; therefore, dark purple nodes represent nodes that have occupancy in both CBDApo and CBDcAMP Markov state models (that is, dark purple nodes overlap).

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588

When comparing the conformational landscapes of CBDApo and CBDcAMP, we noticed two striking featuresthe large number of shared states (Fig. 2) and the character of the unique states. The models indicate that the systems share B70% of the total number of sampled conformational states (Supplementary Table 2). Notably, the shared conformational states of CBDApo and CBDcAMP comprise 66 and 38%, respectively, of the total population for each ensemble.

Many of the unique CBDApo conformational states contain structures in which the cAMP-binding site is formed poorly and, thus, sterically prevented from occurring when cAMP is bound (Supplementary Fig. 11). Many of the most populated unique CBDcAMP conformational states exhibit structures in which the C terminus is interacting with cAMP (Supplementary Fig. 12). These states include the dynamical capping of cAMPs adenosine ring by Y244, part of the C-helix, which substitutes the native capping residue W260 from the B CBD in the full regulatory subunit15,43. These capping interactions with Y244 have been observed experimentally in the truncated RIa (a.a. 1259)43 but not previously observed in simulations. In addition, uorescent experiments studying the dynamics of R239 (located on the C-helix) indicate that the C-helix has a higher propensity to interact with the core cyclic nucleotide-binding domain in the presence of cAMP34. The unique CBDcAMP conformational states are not sterically excluded from the CBDApo conformation ensemble, and interactions between the C terminus and the cAMP-binding sites were observed in the CBDApo ensemble. Therefore, the absence of unique states in the CBDcAMP ensemble is likely due to the absence of stabilizing interactions between cAMP and the C terminus. Importantly, both models sampled the active and inactive states; yet, neither state is the most populated conformational state in either system.

Using the most probable conformational state as the reference state, at physiological conditions, the depth of free-energy landscape is 6.6 kcal mol 1 for CBDApo and 7.6 kcal mol 1 for CBDcAMP. The free energies for the individual conformational states follow a Gaussian distribution with a mean depth of the free-energy landscape at 3.11.4 kcal mol 1 for CBDApo and 4.21.6 kcal mol 1 for CBDcAMP. The transition free energy out of a conformational state to any of its neighbouring conformational states (that is, only states that it is connected to) ranges from 5.8 to 8.6 kcals mol 1 for both systems.

(5.8 kcal mol 1 corresponds to a transition made at approximately the lag time of the Markov state model and, thus, is the lowest bound transition free energy calculated from the Markov state model.) The absolute difference in free energy between the active and inactive states is 0.1 kcal mol 1 for CBDApo and0.9 kcal mol 1 for CBDcAMP. Overall, the free-energy landscapes of both systems provide access to both functional states at physiological conditions, with cAMP binding deepening the free-energy landscape by B1 kcal mol 1.

Transition pathway theory44,45 identied the 10 highest ux paths between the functional conformational states for each ensemble (Supplementary Fig. 13). The proportion of the highest ux path relative to the total ux for each functional state ensemble was striking: For CBDApo, the highest ux path was only 9% of the total ux, and for CBDcAMP, it was only 4% of the total ux (Supplementary Table 3). The lack of a dominant path and the variety of paths identied in the analysis indicate that there are multiple pathways between functional conformational states. Inspection of the pathways indicated no common order of conformational changes or shared features at the bottlenecks of the pathways.

Therefore, to understand the transition between active and inactive conformation states across the free-energy landscape,

we needed to look at the transition pathways collectively. We measured the kinetics of the transition between active and inactive states using mean rst passage times (that is, the fastest average time it takes for the system to move across the free-energy landscape from active-to-inactive state, or vice versa; Table 1). Interestingly, the mean rst passage times for the transition from inactive to active state are similar with or without cAMP bound: B30 ms. However, the addition of cAMP slows the transition from active-to-inactive state by a factor of ve, from 16 to 84 ms.

These results indicate that CBDApo favours the inactive state, which inhibits PKAs catalytic subunit, whereas CBDcAMP favours the state that activates PKAs catalytic subunit. Thus, the dynamics of the cyclic nucleotide-binding domain derived from the Markov state models are consistent with the function of the regulatory subunit in PKA.

Taken together, our results indicate conformational selection as the governing mechanism of allostery in the cyclic nucleotide-binding domain. The weakest evidence for this mechanism is the number of shared states between both models because it is largely a function of the clustering algorithm. However, when the conformations for each system are clustered separately, the results are the same. In fact, the joint clustering approach better captures the principal motions of the CBD at a 3.0- clustering cutoff, as determined by implied timescale plots (Supplementary Fig. 14).

The stronger evidence supporting conformational selection is the change in the free-energy landscape and the corresponding transition kinetics. Weinkam et al.8 classied allosteric mechanisms based on the difference between the change in free energy of two functional states and the binding free energy of the ligand. Although the current Markov state models do not allow determination of the free energy of cAMP binding, and there has been no direct experimental measurements of cAMP binding in our construct, we can estimate the binding free energy of cAMP through previous work on the isolated A cyclic nucleotide-binding domain46. The free energy of binding of cAMP to the isolated A cyclic nucleotide-binding domain is B 3 kcal mol 1

(ref. 46). This value is greater than the difference in free energy between the active and inactive states, B 1 kcal mol 1,

supporting conformational selection as the general mechanism of allostery. In addition, the on-rate for cAMP to the A cyclic nucleotide-binding domain of R, B4.8 104 s 1 (ref. 46), is

faster than the transition times between active and inactive states in the CBDcAMP system. And the fact that the transition rates between inactive and active conformations are independent of ligand binding indicate conformational selection as described by Zhou9. Experimental results23,47 support conformational selection as the mechanism of ligand-induced allostery, thus complementing the interpretations of our Markov state model presented here.

One limitation of our model is the potential bias towards conformational selection because the simulations were initiated at both experimentally determined active and inactive states, which

Table 1 | Mean rst passage time between active and inactive conformations.

System Without cAMP With cAMPA-I (ls) AI (ls) A-I (ls) AI (ls)

CBD 16.1 34.0 83.7 30.3

PBC 2.9 0.9 4.9 0.2 B/C Helix 8.2 13.4 24.3 10.6 N3A 1.3 2.9 2.2 1.6

cAMP, cyclic AMP; CBD, cyclic nucleotide-binding domain; PBC, phosphate-binding cassette.

4 NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588 ARTICLE

CBD-Apo CBD-cAMP

were included in and connected through the Markov state models. However, nothing sterically hinders cAMP from binding either functional conformational state in the absence of the catalytic subunit, and no restraints were required to keep cAMP bound during the simulations. Thus, all four systems used to initiate the simulations are potential members of the conformational ensemble. Starting the simulations in each experimental conformation does not predispose them to becoming members of highly populated or networked states or bias the transition kinetics between active and inactive states.

Details of the allosteric mechanism. While our results indicate conformational selection as the general mechanism of allostery for the CBD, this nding is not wholly novel due to previous experimental work23,47. However, our work provides a foundation to delve further into the mechanisms of allostery at an atomic scale in ways inaccessible to current experimental methods: we can explore structural characteristics of the equilibrium conformational space, and the principal motions of the CBD and its functional macrostates.

To understand the structural characteristics of the conformational landscapes, we built stationary probability density volumes of the positions of the alpha carbons. In other words, we sought to determine the probability of nding the mass of the alpha carbons in a given volume of space at equilibrium. The volumes were determined by dividing the space in 0.1 voxels and calculating the probability of nding an alpha carbon within 3 (the r.m.s.d. cutoff for the Markov state model) of the voxel for each conformational state. By summing over all states, we obtained a map of the stationary probability density, which provides a visual representation of the conformational ensembles.

Fig. 3 shows the probability of density surfaces at 66, 95 and99.9%. The volume enclosed by each surface has a probability greater than the cutoff of being occupied by an alpha carbon. This gure shows the wide range of motion sampled by the a-helical structural subdomain at equilibrium; this range of motion is similar for the CBDApo and CBDcAMP systems. Similar to the graphical representation of the Markov state model, the stationary probability density volumes indicate the similarity of the conformations sampled by CBDApo and CBDcAMP.

However, the differences between the stationary probability density volumes give us new structural insight into the mechanisms of allostery within the CBD. The rst notable difference is the unique density formed by the B/C helix proximal to cAMP in the CBDcAMP ensemble (Fig. 4a). This density corresponds to the interactions between the C terminus and the adenosine ring of cAMP discussed above. A second difference is the higher probability that the N3A motif extends over the B/C helix in the CBDApo ensemble (Fig. 4b). This density corresponds with the experimentally determined binding interface between PKAs catalytic and regulatory subunits16. Interestingly, the CBDApo system only increased density by about 10% in forming the interface. This small variation shows how seemingly small changes in a conformation ensemble give rise to inhibitory function of the CBDApo system by increasing the probability of formation of the binding interface between PKAs catalytic and regulatory subunits.

To determine the principal motions of the CBD, we built kinetic coarse-grained models from the Markov state models of each system to identify metastable states or the principal motions of CBD arising from the conformational ensemble. These models were generated by complementarily employing principal component and Markov clustering analysis25,28,32,42,48,49.

The CBDApo implied timescale plots (Supplementary Fig. 3) suggest at least one dominant slow motion followed by numerous

secondary faster motions. The dominant slow motion of the CBDApo ensemble is between conformations similar to the active and conformations similar to the inactive conformational states (Fig. 5). Of note, the second most dominant motion is integration of the unfolded N-helix (the N3A motif) into the core b-barrel (Supplementary Fig. 15). These conformations are only observed in the CBDApo ensemble as they deform the cAMP-binding site, sterically excluding cAMP. Integration of the unfolded N-helix requires partial unfolding of the core b-barrel, which was observed in previous NMR work21. We recognize that these conformations may be an artifact of the truncated CBD and, therefore, may not be part of the mechanisms of allostery for the regulatory subunit in vivo. Still, the consistency between the results of the simulations and the NMR work indicates that the models are able to replicate structural events seen experimentally.

Like the CBDApo system, the implied timescale plot for the CBDcAMP system indicates one dominant slow motion followed by a number of faster motions (Supplementary Fig. 4). Although the dominant motion is structurally similar to that of CBDApo (that is, it is also the transition between active and inactive conformations (Fig. 5)), it is slower than that of CBDApo (compare Supplementary Figs 3 and 4). Interestingly, the second slowest motion within the CBDcAMP system corresponds to those states where the C terminus is interacting with cAMP.

Reference 66.7% Surface 95.0% Surface 99.9% Surface

Figure 3 | Equilibrium probability density surfaces for CBDApo and CBDcAMP at 66, 95 and 99.9% cutoffs. The surfaces correspond to the probability density surface for each residues alpha carbon. Densities contributed from key structural subdomains are coloured: PBC, red; B/C helix, ice blue; and N3A motif, ochre. See Fig. 1 for the description of CBD representations.

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588

a

b

Figure 4 | Selected difference surfaces between equilibrium probability densities of CBDApo and CBDcAMP. (a) The difference surface for 10% or more density for the N3A Motif of CBDApo over CDBcAMP. Arrows indicate interface between PKAs catalytic and regulatory subunits (b) The different surface for 1% or more density for the B/C Helix of CBDApo over CDBcAMP.

To strengthen the relationship between motions of the CBD and its function, we generated a functional coarse-grained model of the conformational ensembles. The model identies which conformational stats have a higher probability of transitioning to either the active of inactive state, thus indicating the role of each conformational state in the function of the system. Using committer analysis from transition pathway theory44,45, the ensembles were divided into two functional macrostates. Each macrostate was composed of conformational states with a 450%

probability of transitioning to the active or inactive state. In other words, the functional coarse-grain model can be conceptualized as the continental divide of the free-energy landscape, indicating to which functional state (active or inactive) an individual conformational state is connected to most strongly.

The functional coarse-grain model of CBDApo (Fig. 4) contains similar populations for the active and inactive macro-states, with a preference for the inactive macrostate. CBDcAMP (Fig. 4), however, prefers the active macrostate (Table 2). These results agree with NMR studies that show similar inactive and active conformational populations at very-low-cAMP concentration and domination of the active macrostate at high-cAMP concentrations21. In addition, the relationship of the half-life of either functional macrostate is similar to the mean rst passage time between functional conformational states in that transitions from inactive to active states have similar kinetics in both systems, whereas the presence of cAMP slows the transition from active to inactive state (Table 2).

The functional coarse-grain model and the kinetic coarse-grain model closely correlate with each other as determined by

conformational state membership in the functional or kinetic macrostates. Correlation between the two coarse-grained models is expected as they both depend on the transition kinetics of the Markov state models; however, which motion corresponded to the change in functional states is not clear a priori. This correlation indicates that the functional macrostates determined by transition pathway theory correspond to the slowest motion of the cyclic nucleotide-binding domain (Fig. 4). It also supports the assumption that the crystallographic conformations represent the active and inactive states of the cyclic nucleotide-binding domain and that its slowest motion is the transition between active and inactive states.

The coarse-grained models combined with observations of the change in stationary probability densities for the two systems support a general explanation that the modication of the transition kinetics by cAMP gives rise to function. In CBDApo, transitions between states occur with a tendency towards the inactive conformation, which increases the likelihood of generating the interface that allows for the regulatory subunit to bind the catalytic subunit leading to inactivation of PKA. With cAMP bound, the transition between inactive and active states is slower, which is likely the result of interactions formed between the C terminus and cAMP. The preference for the active-like conformations in the CBDcAMP system decreases the probability of forming the binding interface and allowing activation of PKA.

Motions of the signalling motifs. To gather more detail on the allosteric mechanism, we examined the propagation of an allosteric signal through the helical signalling motifs of the cyclic nucleotide-binding domain by building new Markov state models for the three principal signalling motifs using the same MD trajectories (apo and cAMP-bound). While the motions of the signalling motifs are captured in the Markov state models of the full CBD, as discussed above, individual models elucidate the conformational ensembles of focused structural elements and provide an unobstructed view of their modulation by cAMP. The Markov state models for the individual motifs were developed and analysed using the same methods for the full system (Supplementary Figs 510). They indicate conformational selection as the general mechanism of allostery, due to well-connected active and inactive conformational states (Fig. 6) with multiple pathways across the free-energy landscape for each of the three key motifs (Supplementary Table 3).

When we compared the conformational state distribution between the apo and cAMP-bound systems, we observed that the N3A motif and the B/C helix share a majority of their conformational states (Fig. 6 and Supplementary Table 2), similar to the full cyclic nucleotide-binding domain. The B/C helix generally explores the conformational space between the experimentally determined active and inactive structures (Fig. 6a). However, the N3A motif explores conformations dissimilar to either experimental conformation (Fig. 6c). These models are consistent with the observed motions of these domains (Supplementary Fig. 2), particularly with the unfolding of the N-helix (Supplementary Fig. 11).

These results support and help explain why there is no observed change in chemical shifts in NMR experiments21 for the N3A motif or the B/C helix despite the changes in the crystallographic structures (Fig. 1). Simply put, in solution the conformational space explored by the N3A motif and the B/C helix is largely independent of cAMP, though cAMP does modulate the relative populations of the conformational states and their kinetics.

In contrast, the PBC showed signicant changes in chemical shifts on binding cAMP21. While these changes appear to be

6 NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588 ARTICLE

Kinetic coarse-grain model Functional coarse-grain model

CBD-Apo

RMSD from inactive conformation ()

10

8

6

4

RMSD from inactive conformation ()

RMSD from inactive conformation ()

10

8

6

4

10

8

6

4

2 0 4 6 8 10 12 2

2 0 4 6 8 10 12

2

2

RMSD from active conformation () RMSD from active conformation ()

RMSD from active conformation () RMSD from active conformation ()

RMSD from inactive conformation ()

CBD-cAMP

10

8

6

4

2 0 4 6 8 10 2

2 0 4 6 8 10

Figure 5 | Comparison of conformational state assignment for kinetic and functional coarse-grain models. Conformational state metastable states for the rst/slowed motions in the kinetic coarse-grain Markov state model graph are coloured below according to macrostate membership (red, white). Functional coarse-grain models are coloured either white or maroon, depending on whether the state would go to the active or inactive conformations, respectively.

Table 2 | Population and half-life for CBDApo and CBDcAMP coarse-grain functional states.

Macrostate CBDApo population (%) t1/2 mean (s.d.; ls) Macrostate CBDcAMP population (%) t1/2 mean (s.d.; ls) Active 36 2.29 (0.09) Active 77 36.59 (1.12)

Inactive 64 8.62 (0.41) Inactive 23 7.41 (0.45)

cAMP, cyclic AMP; CBD, cyclic nucleotide-binding domain.

driven by interaction with cAMP, the Markov state model of the PBC illuminates another aspect of these changes: that, on binding of cAMP, the PBC loses about half of its conformational states, representing a signicant loss in conformational entropy (Fig. 6b, Supplementary Table 2). In addition, the kinetic and functional coarse-grain models of the PBC show signicant changes between apo and cAMP systems (Supplementary Fig. 17 and Supplementary Table 5) unlike the other signalling motifs (Supplementary Figs 16 and 18 and Supplementary Tables 4,6). Without cAMP, the PBC has a shallow free-energy landscape with quick transition between metastable macrostates (Table 1, Supplementary Fig. 17, and Supplementary Table 5). With cAMP bound, there is one dominant (slow) transition between two metastable macrostates (Table 1, Supplementary Fig. 17 and Supplementary Table 5; also compare implied timescale plots in Supplementary Figs 7 and 8). Unlike the other models (Fig. 5 and Supplementary Figs 16 and 18), the conformational states that comprise these metastable macrostates do not correspond to the division of microstates seen in the functional coarse-grain model (Supplementary Fig. 17). Instead the functional coarse-grain model indicates a strong propensity of the PBC to assume a cAMP-binding conformation or active conformation even in the absence of cAMP (Supplementary Fig. 17 and Supplementary Table 5). Overall, the model shows that, without cAMP bound,

the PBC explores a variety of conformations, but the binding site has a tendency to assume the active cAMP-bound conformation. cAMP binding selects the active conformation while allowing transition between active and inactive state but signicantly slows the transition between these states (Table 1, Supplementary Table 5).

The spatial arrangement of the helical signalling motifs suggests that the propagation of the allosteric signal would be from the PBC to the B/C helix to the N3A motif (Fig. 1). Even the natures of the conformation ensembles, described by the Markov state models, suggest propagation of the allosteric signal from a switch like the PBC out to the conformationally diverse N3A motif.

However, the kinetics of the signalling motifs suggest another mechanism. The mean rst passage times between active and inactive states suggest that the movement of the B/C helix is the rate-limiting step for the functional transition because it has the slowest transitions between active and inactive states (Table 1). Like the full CBD, the B/C helix has similar transition times between active and inactive states whether or not cAMP is bound. However, cAMP signicantly decreases the rate of the active-to-inactive transition (Table 1). This rate reduction is due to the stabilization of the active conformation by cAMP through interactions between the C helix and cAMP (Fig. 4), as observed and discussed above.

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588

a

B/C helix

Phosphate-binding cassette (PBC)

N3A motif

25

conformation () RMSD from inactive conformation ()RMSD from inactive conformation ()

20

15

10

5

0 5 10 15 20 25

RMSD from active conformation ()

b

6

5

4

3

2

1

1 2 3 4 5

RMSD from active conformation ()

RMSD from active conformation ()

c

12

10

RMSD from inactive

8

6

4

0 5 10 20

Figure 6 | Overlay of Markov State Model Conformational State Graphs for the SubDomains of cyclic nucleotide-binding Domain with (CBDcAMP, magenta) and without (CBDApo, cyan) cAMP Bound. Markov state models conformational state graphs use the same conventions dened in Fig. 2.

In PKA and other proteins, the experimental free energy of binding of cAMP to the cyclic nucleotide-binding domain is dominated by enthalpic contributions5052, suggesting that the contacts between cAMP and the PBC and B/H helix are drivers for this phenomenon. Our model indicates that the binding of cAMP near the PBC dampens the dynamics of that signalling motif, greatly increasing its propensity to adopt an active-like conformation, and that the motions of the B/C helix are the rate-limiting step in this transition.

DiscussionOverall, we nd that the prototype cyclic nucleotide-binding domain is a highly dynamic system and that cAMP modulates its function by modifying the dynamics of the signalling motifs. By comparing CBDApo and CBDcAMP ensembles, we see that the free-energy landscape allows access to both active and inactive states through an ensemble of transition pathways (that is, there is no one dominant pathway between states). Collectively, our analysis supports conformational selection as the general mechanism of allostery in this system. The addition of cAMP modulates the transition rates between the functional states but only between active-to-inactive states and not vice versa. This modulation of rates occurs by slowing motion in the CBDcAMP ensemble, which is caused by interactions with cAMP. Finally, our work indicates that the change in dynamics of the B/C helix is the rate-limiting step in adopting a fully active conformation. Furthermore, the motion of the B/C helix is essential for signal propagation to the N3A motif and formation of the interface with the catalytic subunit. This mechanism could be tested experimentally though NMR by observing shifts in the relative populations of active and inactive ensembles caused by mutations in the C terminus, specically, mutations of residue 244, which should modify interactions with cAMP.

While the general mechanism of allostery for the cyclic nucleotide-binding domain is conformational selection, some elements of the mechanism appear to behave more like induced t. Particularly notable in this regard are the interactions between the C-helix and cAMP, which generate multiple conformational states that are not sampled in the apo state, as compared with the motions of the individual signalling motifs that appear more similar to entropic mechanisms of allostery. Furthermore, in the context of activation of the full kinase, these mechanisms may change, as interaction with the catalytic subunit and the remaining portions of the regulatory subunit will almost certainly modulate the conformational ensemble of the cyclic nucleotide-binding domain. This work suggests that a clearer picture of the free-energy landscape, achieved through an integrated experimental and computational approach, is required to fully understand allosteric behaviour.

Although a number of our conclusions were already observed experimentally, the power of our computational Markov-state-model-based approach is that it unies multiple lines of existing experimental data into a single, cohesive framework and, thereby, achieves a more complete understanding of allostery. This approach builds on an emerging paradigm for understanding protein function and dynamics based on computational exploration and visualization of protein conformational ensembles and their underlying functional free-energy landscapes in atomic detail.

Methods

System preparation. Parameterization of cAMP. We parameterized cAMP for the MD simulations as a small organic molecule using the generalized amber force eld53 to allow for consistent treatment of ligands in future studies. The coordinates for cAMP were taken from the 1RGS15 crystal structure. Partial atomic charges for cAMP were determined by single-point energy calculations using Schrodingers QM module Jaguar (Suite 2012: Jaguar, version 7.9, Schrdinger, LLC, New York, NY, 2012) using the HartreeFock level of theory and 6-311g** bases. Generalized amber force eld atom types were assigned using antechamber54, and the parameters les were prepared with AMBERs xleap.

MD System preparation. Using Schrodingers Maestros (Suite 2012: Maestro, version 9.3, Schrdinger, LLC, New York, NY, 2012) PDB prep and the Desmond System Preparer, we prepared four systems for MD simulations, CBDApo and CBDcAMP in both the active and inactive conformations. The heavy-atom atomic coordinates for the active and inactive conformations of CBDApo were taken for crystal structures 1RGS15 and 2QCS16 resides 119 to 244, respectively. We added cAMP to the inactive conformation by aligning the cAMP-binding site of 1RGS and 2QCS and cutting and pasting the coordinates. For the CBDApo active conformation, the coordinates for cAMP were deleted from 1RGS. Within

8 NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588 ARTICLE

PBD prep, each system was capped to remove terminal charges, the systems were protonated at pH 7.0 with the pKa titratable residues determined with the

Maestro-integrated PROPKA, and, due to the low resolution of 1RGS, all crystallographic waters were removed. Using the Desmond System Preparer,we solvated each system in a cubic water box using TIP3 waters55, counter ions, and 120 mM NaCl. No restraints were placed on cAMP.

MD simulations. System parameterization. System coordinates from Maestro were converted into an xleap-readable format using in-house Python scripts. Each system was parameterized in xleap using the Amber99SB56 force eld and periodic boundary conditions were implemented.

GPU-enabled MD simulations. All systems were minimized and equilibrated using the GPU version of Amber12 (refs 3840). We minimized the system in four stages: (1) 500 steps of hydrogen-only minimization, (2) 500 steps of solvent minimization, (3) 500 steps with only the backbone constrained and (4) 5,000 steps of full minimization. We equilibrated the system using harmonic equilibration at 310 K over four sequential 500-ps runs, decreasing the restraint potential on the backbone each step. GPU-enabled AMBER12 production runs were carried out as an NTP ensemble at 310 K and 1 bar with a 2-fs time-step and partial mesh Ewald electrostatic approximation. MD input les are provided as part of the data sharing.

Anton MD simulations. MD simulations on Anton37 were performed on the same parameterized, minimized and equilibrated systems as GPU-enabled AMBER12 simulations. The Anton simulation was run in the NTP ensemble, using Antons Berendsen thermostat-barostat, at 310 K and 1 bar with a 2-fs time-step and partial mesh Ewald electrostatic approximation.

Initial and adaptive sampling. We sampled the conformational ensemble from the equilibrated conformation of each system, starting each run with a different set of initial velocities. Supplementary Table 1 contains the total sampling time for each system used in this analysis.

In conjunction with model building, we performed multiple adaptive sampling runs. These runs were initiated from preliminary Markov state model conformational states with o3 members. These conformational states corresponded to clustered regions of the conformational ensemble poorly sampled by the MD simulations and, therefore, the transitions to and from those states were poorly sampled. For each under-sampled conformational state, one conformation was selected. From that conformation, three new 1015 ns MD simulations were run, each with a new initial velocity using GPU-enabled AMBER12. The short MD simulations were sufcient to explore the local conformational landscape and return to a local energy minimum.

Building Markov state models. Trajectory preparation. MD trajectories were processed using CPPTRAJ57 and VMD58. All frames were aligned using Ca of residues 153 to 199 and 211 to 223, the stable b-barrel, to crystallographic active conformation. The alignment provided a common reference frame to compare all conformations while maximizing the difference between conformations. The frame rate for the trajectories was unied to 120 ps due to different conformational sampling rates between the GPU-enabled AMBER12 (0.5 ps) and Anton (120 ps) MD simulations. The trajectories were converted into NAMDs59 .dcd trajectory format for analysis with MSMBuilder2.

Building the Markov State Model. We used MSMBuilder2 (refs 32,42), release2.51, to preform cluster analysis, Markov state model building and Markov state model selection. To have a common set of microstates for both the CBDApo and CBDcAMP systems, all trajectories were clustered together using the hybrid k-centres k-medoids clustering algorithm using the custom distance metric option to calculate r.m.s.d. without rst performing an alignment on the atomic coordinates (an option not available in MSMBuilder2). MSMBuilder2 allows for renement of clustering through renement of frame membership within a cluster as well as global renement of cluster generators. For our systems, we found that any renement only made the space discretization worse, as determined by implied timescale plots. Therefore, we didnt use any renement options for our clustering. Markov state models were then built on the CBDApo and CBDcAMP trajectories separately using cluster generators to assign state membership to each frame of the trajectories. When the best conformational landscape partition and lag time were determined (see discussion below), a nal Markov state model was built using the maximal likelihood estimator option in MSMBuilder2.

Model selection. For the Markov state model of the whole systems, we selected a clustering cutoff distance of 3.0 and a lag time of 80 steps or 9.6 ns. For consistency in comparing models, a lag time of 9.6 ns was used for all.

Markov state model visualization. The Markov state models were visualized using in-house Python scripts leveraging the publically available NetworkX module60.

Markov state model analysis. Equilibrium population. The equilibrium population was determined from the rst eigenvalue of the Markov state models transition probability matrix that corresponds to the equilibrium distribution of the conformational states.

Transition pathway theory analysis. Transition pathway theory analysis was preformed using MSMBuilder2 transition pathway theory scripts42. Total ux was calculated as described in Prinz et al.45. Transition pathway theory-based clustering was preformed using in-house Python scripts leveraging MSMBuilder2 modules.

The script assigned transition pathway theory macrostates to each microstate of the original Markov state model based on the probability of transitioning to either the inactive or the active conformation (cutoff 50%). The conformational state assignments for each frame of the trajectories were reassigned using in-house Python scripts and a new Markov state model based on the macrostate assignments was built with MSMBuilder2.

Markov cluster analysis macrostate models. Markov cluster analysis clustered the states of the Markov state model based on their kinetic relationships within the Markov state model using a random walk algorithm to determine local sinks in the Markov state model48. We performed the Markov cluster analysis using in-house Python scripts that assigned new Markov cluster analysis state assignments to the original state assignments and built a Markov state model. In addition, we used MSMBuilders2s PCCA (ref. 42) module to validate the division identied

through Markov clustering.

Mean rst passage time. Absolute MFPT between the putative functional microstates was calculated as described in Singhal et al.61

Half-life. The half-life of coarse-grained microstates was calculated using a kinetic Monte Carlo scheme on the Markov state models transitions matrix as described in Shukla et al.31 Trajectories were generated by randomly selecting a starting state within a macrostate, then moving between states based on the transitions probability as determined by the Markov state model. The passage time of the trajectory was the number of steps in the trajectory multiplied by the lag time(9.6 ns). The half-life was calculated by taking the average of 10 median passage times from a collection of 1,000 trajectories.

Calculating kinetics and free energies. Assuming that transitioning out of one conformational state into its neighbouring conformational states can be modelled as parallel reactions, we can treat the kinetics for the conformational changes between connected conformational states as rst-order reactions. Therefore, we determined the rate constant using

k

lnp1=p0 t

where k in the rst-order rate constant, p0 is the probability of starting in a given state (that is, the starting concentration or 1), p1 is the probability of still being in that state at a given lag time (that is, the concentration after a give time) and t is the lag time (9.6 ns). The free energy of transition was calculated using Eyrings equation as follows:

DG RT ln

h kkBT k

where R is the gas constant, h is Planks constant, kB is Boltzmanns gas constant, k is the transmission coefcient assumed to be 1, T is temperature (310 K) and k is the rate constant calculated above.

The free energy of a single microstate was calculated using

DG RT ln

pi pr

where pi is the probability of a given microstate, pr is the probability of the reference microstate (the most probable microstate in the ensemble), R is the gas constant (0.0012 kcal mol 1 K 1) and T is temperature (310 K)62.

Markov state model signalling motifs. The Markov state models of the signalling motifs were built and analysed using the same methods described above except we only used a subset of Ca. The residues used for the analysis of each subdomain were: 119 to 151 for the N3A motif, 199 to 211 for the PBC and 226 to 244 for the B/C helix. Markov state models were built with a clustering cutoff of 4.5 for the N3A motif, 1.5 for the PBC and 4.5 for the B/C helix. The lag time used was9.6 ns.

Data sharing. MD trajectories, sample MD input scripts and the Markov state model analysis scripts are available for download at http://dx.doi.org/10.6075/J0Z60KZS

Web End =http://dx.doi.org/10.6075/J0Z60KZS .

References

1. Changeux, J.-P. The feedback control mechanisms of biosynthetic L-threonine disminase by L-isoleucine. Cold Spring Harb. Symp. Quant. Biol. 26, 313318 (1961).

2. Monod, J. & Jacob, F. Teleonomic mechanisms in cellular metabolism, gowth, and differentiation. Cold Spring Harb. Symp. Quant. Biol. 26, 389401 (1961).

3. Monod, J., Wyman, J. & Changeux, J. On the nature of allosteric transitions.J. Mol. Biol. 12, 88118 (1965).4. Koshland, D. E., Nmethy, G. & Filmer, D. Comparison of experimental binding data and theoretical models in proteins containing subunits. Biochemistry 5, 365385 (1966).

5. Pauling, L. The oxygen equilibrium of hemoglobin and its structural interpretation. Proc. Natl Acad. Sci. USA 21, 186191 (1935).

6. Cui, Q. & Karplus, M. Allostery and cooperativity revisited. Protein Sci. 17, 12951307 (2008).

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588

7. Motlagh, H. N., Wrabl, J. O., Li, J. & Hilser, V. J. The ensemble nature of allostery. Nature 508, 331339 (2014).

8. Weinkam, P., Pons, J. & Sali, A. Structure-based model of allostery predicts coupling between distant sites. Proc. Natl Acad. Sci. USA 109, 48754880 (2012).

9. Zhou, H.-X. From induced t to conformational selection: a continuum of binding mechanism controlled by the timescale of conformational transitions. Biophys. J. 98, L15L17 (2010).

10. Tsai, C. J., Kumar, S., Ma, B. & Nussinov, R. Folding funnels, binding funnels, and protein function. Protein Sci. 8, 11811190 (1999).

11. Boehr, D. D., Nussinov, R. & Wright, P. E. The role of dynamic conformational ensembles in biomolecular recognition. Nat. Chem. Biol. 5, 789796 (2009).

12. Feher, V. A., Durrant, J. D., Van Wart, A. T. & Amaro, R. E. Computational approaches to mapping allosteric pathways. Curr. Opin. Struct. Biol. 25, 98103 (2014).

13. Collier, G. & Ortiz, V. Emerging computational approaches for the study of protein allostery. Arch. Biochem. Biophys. 538, 615 (2013).

14. Taylor, S. S., Ilouz, R., Zhang, P. & Kornev, A. P. Assembly of allosteric macromolecular switches: lessons from PKA. Nat. Rev. Mol. Cell Biol. 13, 646658 (2012).

15. Su, Y. et al. Regulatory subunit of protein kinase A: structure of deletion mutant with cAMP binding domains. Science 269, 807813 (1995).

16. Kim, C., Cheng, C. Y., Saldanha, S. A. & Taylor, S. S. PKA-I holoenzyme structure reveals a mechanism for cAMP-dependent activation. Cell 130, 10321043 (2007).

17. Sjoberg, T. J., Kornev, A. P. & Taylor, S. S. Dissecting the cAMP-inducible allosteric switch in protein kinase A RIalpha. Protein Sci. 19, 12131221 (2010).

18. Kornev, A. P., Taylor, S. S. & Ten Eyck, L. F. A generalized allosteric mechanism for cis-regulated cyclic nucleotide binding domains. PLoS Comput. Biol. 4, e1000056 (2008).

19. Kim, C., Xuong, N.-H. H. & Taylor, S. S. Crystal structure of a complex between the catalytic and regulatory (RIalpha) subunits of PKA. Science 307, 690696 (2005).

20. Das, R. & Melacini, G. A model for agonism and antagonism in an ancient and ubiquitous cAMP-binding domain. J. Biol. Chem. 282, 581593 (2007).

21. Das, R. et al. cAMP activation of PKA denes an ancient signaling mechanism. Proc. Natl Acad. Sci. USA 104, 9398 (2007).

22. Das, R., Abu-Abed, M. & Melacini, G. Mapping allostery through equilibrium perturbation NMR spectroscopy. J. Am. Chem. Soc. 128, 84068407 (2006).

23. Akimoto, M. et al. Signaling through dynamic linkers as revealed by PKA. Proc. Natl Acad. Sci. USA 110, 1423114236 (2013).

24. Boras, B. W., Kornev, A., Taylor, S. S. & McCulloch, A. D. Using Markov state models to develop a mechanistic understanding of protein kinase A regulatory subunit RIa activation in response to cAMP binding. J. Biol. Chem. 289, 3004030051 (2014).

25. Pande, V. S., Beauchamp, K. & Bowman, G. R. Everything you wanted to know about Markov State Models but were afraid to ask. Methods 52, 99105 (2010).

26. Prinz, J. H. et al. Markov models of molecular kinetics: generation and validation. J. Chem. Phys. 134, 174105 (2011).

27. Chodera, J. D. J. & No, F. Markov state models of biomolecular conformational dynamics. Curr. Opin. Struct. Biol. 25C, 135144 (2014).

28. Malmstrom, R. D., Lee, C. T., Van Wart, A. T. & Amaro, R. E. Application of molecular-dynamics based markov state models to functional proteins. J. Chem. Theory Comput. 10, 26482657 (2014).

29. Lane, T. J., Bowman, G. R., Beauchamp, K., Voelz, V. A. & Pande, V. S. Markov state model reveals folding and functional dynamics in ultra-long MD trajectories. J. Am. Chem. Soc 133, 1841318419 (2011).

30. Yang, S., Banavali, N. K. & Roux, B. Mapping the conformational transition in Src activation by cumulating the information from multiple molecular dynamics trajectories. Proc. Natl Acad. Sci. USA 106, 37763781 (2009).

31. Shukla, D., Meng, Y., Roux, B. & Pande, V. S. Activation pathway of Src kinase reveals intermediate states as targets for drug design. Nat. Commun. 5, 3397 (2014).

32. Bowman, G. R., Huang, X. H. & Pande, V. S. Using generalized ensemble simulations and Markov state models to identify conformational states. Methods 49, 197201 (2009).

33. Kohlhoff, K. J. et al. Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways. Nat. Chem. 6, 1521 (2014).34. Anand, G. S. et al. Cyclic AMP- and (Rp)-cAMPS-induced conformational changes in a complex of the catalytic and regulatory (RI{alpha}) subunits of cyclic AMP-dependent protein kinase. Mol. Cell. Proteomics 9, 22252237 (2010).

35. Berman, H. M. et al. The cAMP binding domain: an ancient signaling module. Proc. Natl Acad. Sci. USA 102, 4550 (2005).

36. Kannan, N. et al. Evolution of allostery in the cyclic nucleotide binding module. Genome Biol. 8, R264 (2007).

37. Shaw, D. E. et al. Atomic-level characterization of the structural dynamics of proteins. Science 330, 341346 (2010).

38. Case, D. A. et al. Amber 12 (University of California, 2012).

39. Gotz, A. W. et al. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 8, 15421555 (2012).

40. Salomon-Ferrer, R., Goetz, A. W., Poole, D., Le Grand, S. & Walker, R. C. Routine micorsecond molecular dynamics simuations with AMBER - Part II: Particle Mesh Ewald. J. Chem. Theory Comput. 9, 38783888 (2012).

41. Bowman, G. R., Ensign, D. L. & Pande, V. S. Enhanced modeling via network theory: adaptive sampling of markov state models. J. Chem. Theory Comput. 6, 787794 (2010).

42. Beauchamp, K. A. et al. MSMBuilder2: modeling conformational dynamics at the picosecond to millisecond scale. J. Chem. Theory Comput. 7, 34123419 (2011).

43. Ringheim, G. E., Saraswat, L. D., Bubis, J. & Taylor, S. S. Deletion of cAMP-binding site B in the regulatory subunit of cAMP-dependent protein kinase alters the photoafnity labeling of site A. J. Biol. Chem. 263, 1824718252 (1988).

44. Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory for markov jump processes. Multiscale Model. Simul. 7, 11921219 (2009).

45. Prinz, J. H., Keller, B. & Noe, F. Probing molecular kinetics with Markov models: metastable states, transition pathways and spectroscopic observables. Phys. Chem. Chem. Phys. 13, 1691216927 (2011).

46. Herberg, F. W., Taylor, S. S. & Dostmann, W. R. Active site mutations dene the pathway for the cooperative activation of cAMP-dependent protein kinase. Biochemistry 35, 29342942 (1996).

47. Badireddy, S. et al. Cyclic AMP analog blocks kinase activation by stabilizing inactive conformation: conformational selection highlights a new concept in allosteric inhibitor design. Mol. Cell. Proteomics 10, M110.004390 (2011).48. Berezovska, G., Prada-Gracia, D., Mostarda, S. & Rao, F. Accounting for the kinetics in order parameter analysis: lessons from theoretical models and a disordered peptide. J. Chem. Phys. 137, 194101 (2012).

49. Chodera, J. D., Singhal, N., Pande, V. S., Dill, K. A. & Swope, W. C. Automatic discovery of metastable states for the construction of Markov models of macromolecular conformational dynamics. J. Chem. Phys. 126, 155101 (2007).

50. Moll, D., Schweinsberg, S., Hammann, C. & Herberg, F. W. Comparative thermodynamic analysis of cyclic nucleotide binding to protein kinase A. Biol. Chem. 388, 163172 (2007).

51. Mller, S. et al. Cyclic nucleotide mapping of hyperpolarization-activated cyclic nucleotide-gated (HCN) channels. ACS Chem. Biol. 9, 11281137 (2014).

52. Son, I., Selvaratnam, R., Dubins, D. N., Melacini, G. & Chalikian, T. V. Ultrasonic and densimetric characterization of the association of cyclic AMP with the cAMP-binding domain of the exchange protein EPAC1. J. Phys. Chem. B 117, 1077910784 (2013).

53. Wang, J. M., Wolf, R. M., Caldwell, J. W., Kollman, P. A. & Case, D. A. Development and testing of a general amber force eld. J. Comput. Chem. 25, 11571174 (2004).

54. Wang, J. M., Wang, W., Kollman, P. A. & Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 25, 247260 (2006).

55. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926935 (1983).

56. Hornak, V. et al. Comparison of multiple Amber force elds and development of improved protein backbone parameters. Proteins 65, 712725 (2006).

57. Roe, D. R. & Cheatham, T. E. PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 30843095 (2013).

58. Humphery, W., Dalke, A. & Schulten, K. VMD-Visual Molecular Dynamics.J. Molec. Graph 14, 3338 (1996).59. Phillips, J. C. et al. Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 17811802 (2005).

60. Hagberg, A. A., Schult, D. A. & Swart, P. J. In: Proc. 7th Python Sci. Conf. (eds Varoquaux, G., Vaught, T. & Millman, J.) 1115 (2008).

61. Singhal, N. & Pande, V. S. Error analysis and efcient sampling in Markovian state models for molecular dynamics. J. Chem. Phys. 123, 204909 (2005).

62. Altis, A., Otten, M., Nguyen, P. H., Hegger, R. & Stock, G. Construction of the free energy landscape of biomolecules via dihedral angle principal component analysis. J. Chem. Phys. 128, 245102 (2008).

Acknowledgements

We thank Dr Victoria Feher and Professor Giuseppe Melacini for their critical reviews of the manuscript and constructive feedback and Stephanie Sides for her editorial assistance. This work was funded in part by the National Institutes of Health (NIH) through the NIH Directors New Innovator Award Program (DP2-OD007237 to REA, R01-GM034921 to S.S.T.) and an NSF XSEDE supercomputer resources grant(RAC CHE060073N to R.E.A.). Funding from the National Biomedical Computation Resource (NBCR), NIH P41 GM103426 from the National Institutes of Health, is gratefully acknowledged. Anton computer time was provided by the Pittsburgh Supercomputing Center (PSC) and the National Center for Multiscale Modelingof Biological Systems (MM Bios) (grant P41GM103712-S1 from the National

10 NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8588 ARTICLE

Institutes of Health). D.E. Shaw Research generously made the Anton machine at PSC available.

Author contributions

This work is the production of collaboration between the labs of R.E.A. and S.S.T. R.E.A. devised and oversaw the computational components of the work; S.S.T. oversaw the biological aspects of the work. A.P.K. and S.S.T. provided key insight into P.K.A. and R.D.M. R.E.A. A.P.K. and S.S.T. interpreted the models. R.D.M. developed the methods, performed all the simulations and analysis, and was the principal author on the paper. R.E.A. also helped write the paper; A.P.K. and S.S.T. provided edits.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Malmstrom, R. D. et al. Allostery through the computational microscope: cAMP activation of a canonical signalling domain. Nat. Commun. 6:7588 doi: 10.1038/ncomms8588 (2015).

NATURE COMMUNICATIONS | 6:7588 | DOI: 10.1038/ncomms8588 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

& 2015 Macmillan Publishers Limited. All rights reserved.