Using Markov State Models to Develop a Mechanistic Understanding of Protein Kinase A Regulatory Subunit RIα Activation in Response to cAMP Binding*

Background: Binding four cAMP molecules activates protein kinase A (PKA). Results: Mechanistic Markov models (MM) show that although one cAMP activates one catalytic subunit four are necessary for full activation. Conclusion: Conformational selection is the predominant mechanism in PKA activation, but heterodimer interactions are also required. Significance: This MM provides a mechanistic foundation for systems models of PKA signaling networks. Protein kinase A (PKA) holoenzyme consists of two catalytic (C) subunits and a regulatory (R) subunit dimer (R2C2). The kinase is activated by the binding of cAMPs to the two cyclic nucleotide binding domains (CBDs), A and B, on each R-subunit. Despite extensive study, details of the allosteric mechanisms underlying the cooperativity of holoenzyme activation remain unclear. Several Markov state models of PKA-RIα were developed to test competing theories of activation for the R2C2 complex. We found that CBD-B plays an essential role in R-C interaction and promotes the release of the first C-subunit prior to the binding to CBD-A. This favors a conformational selection mechanism for release of the first C-subunit of PKA. However, the release of the second C-subunit requires all four cAMP sites to be occupied. These analyses elucidate R-C heterodimer interactions in the cooperative activation of PKA and cAMP binding and represent a new mechanistic model of R2C2 PKA-RIα activation.

cAMP-dependent protein kinase A (PKA) is a ubiquitous cAMP receptor in eukaryotes (1). The holoenzyme consists of a regulatory (R) 2 subunit dimer and two catalytic (C) subunits (2,3). In the absence of cAMP, the R 2 dimer and two C-subunits interact to form a holoenzyme structure that inhibits kinase activity by binding and blocking the active sites of the C-subunit. Each R-subunit contains two tandem cAMP binding domains (CBDs): CBD-A and CBD-B (2). Activation of PKA is a well known example of allosteric regulation with cAMP binding to the four CBDs, weakening the interactions of the R 2 dimer with the C-subunits and leading to the unleashing of the catalytic activity of the C-subunit (4).
Allostery is a regulatory mechanism in which an effector molecule (ligand) binds at a site other than the active site of the enzyme (5). Binding of the ligand is most commonly associated with conformational changes in the enzyme that lead to the change of its activity. There are two main alternative allosteric activation mechanisms: "induced fit," proposed by Koshland (6), suggests that binding of a ligand precedes conformational changes in the protein structure, whereas "conformational selection," also known as the Monod, Wyman, Changeux model (7), holds that conformational changes in the binding protein precede ligand binding.
Recent studies of a single CBD showed that conformational selection is likely to be the dominant mechanism in the case of cAMP binding to PKA (8). This finding, however, does not fully explain the mechanisms of PKA activation for several reasons. First, the R 2 C 2 holoenzyme contains four binding sites for cAMP that are capable of communicating between each other, creating a cooperative effect of cAMP binding (9). Second, each CBD undergoes major structural changes during activation (10), and two of them (both CBD-As) have extended interfaces with the C-subunits that can interfere with the conformational changes (11,12). Thus, the unleashing of the C-subunits in the R 2 C 2 complex must also be a factor in the activation process. Because of the complexity of the system, it has been difficult to expand the NMR and crystal structure studies beyond a single CBD or an R-C heterodimer to the whole R 2 C 2 complex (12,13).
Here Markov state models (MMs) of PKA-RI␣ were used to test competing allosteric mechanisms of activation and cooperativity against currently available experimental data. These models show that R-C heterodimer interactions are essential to explain observations of the activity of PKA. Neither CBD-A nor CBD-B binding alone promotes the dissociation of the R 2 dimer and C-subunits, but instead, both increase the K D of the R-C complex. However, the models suggest that the first C-subunit activation is regulated by CBD-B binding and that the second activation is regulated by CBD-A binding. These results are important for understanding the roles of conformational selection and induced fit in the mechanisms of PKA activation and may be useful in drug design for targeting PKA and other cAMP-binding proteins.

Model Formation
Five MMs 3 were developed to probe different mechanistic theories of PKA-RI␣ activation. Briefly, each combination of binding partners (R, C, and cAMP) was considered to be a Markov state (Fig. 1). These states represent various molecular conformations that affect the transitions to and from each state. State transitions represent binding events, each defined by simple mass action kinetics. Each transition was considered to be reversible and therefore subject to thermodynamic constraints around a closed loop. For these models, we assumed based on previous studies that the CBD-B site must be bound first in any R-C conformation (9,14). After the first cAMP binding, to the best of our knowledge, the second cAMP binds the CBD-A site of the bound R-subunit, making the order of binding B-A-B-A for all R-C conformations. However, for R-subunit not bound to C-subunit, both the other CBD-B site and the adjacent CBD-A site are available and allow B-B-A-A binding. There-fore, in the R 2 C 2 models, both B-A-B-A binding and B-B-A-A binding were considered for free R 2 dimer without C-subunit, but only B-A-B-A was included for R 2 C 2 and R 2 C 1 states. The addition of kemptide as a buffer of free C-subunit was considered as an extension of the models.

Data Sets Used for Fitting
Five published data sets from two published studies (15,16) were chosen to test the various mechanisms in the model. The data were obtained under physiologically relevant conditions at high Mg 2ϩ and high ATP concentrations to limit additional effects of Mg 2ϩ binding and ATP/ADP interactions. Furthermore, all experiments were conducted for long enough periods to reach thermodynamic equilibrium (17). Data sets I-III are PKA catalytic activity measurements made in the presence of progressively higher concentrations of R 2 dimer with a kemptide substrate. Briefly, 0.25 nM C-subunit was combined with 3, 30, or 300 nM RI␣ subunit. Increasing amounts of cAMP were added and incubated at 37°C for 40 min (15).
Data set IV tested the R-cAMP exclusive interactions in non-cAMP-saturated conditions in the absence of C-subunit. Briefly, 0.25 nM R-subunit was incubated with [ 3 H]cAMP under strict equilibrium conditions. Size exclusion chromatography was used to separate the free cAMP from the bound (16).
Data set V tested cAMP binding to the R 2 C 2 holoenzyme in the presence of an excess of C-subunit (5 M C-subunit and 0.25 nM RI␣). Again size exclusion chromatography was used to separate the free from bound cAMP (16). . The Markov state model of PKA-RI␣ R 2 C 2 holoenzyme. A representation of MM states for the activation of PKA-RI␣ R 2 C 2 holoenzyme by cAMP is shown. Arrows represent the dominant pathways during activation for each of the three tested R 2 C 2 regulatory mechanistic models. Grayed out states are those that are either not populated or only populated when R 2 dimer is in excess compared with C-subunit. The two R-and C-subunits are identical, but for simplicity of naming, the first R-subunit to bind C-subunit is named R, whereas the first R-subunit to bind cAMP is RЈ. The F factors are the thermodynamic constraints that regulate each loop. These factors are described in more detail in Table 1.

Computational Modeling and Fitting Algorithm
All MMs were implemented in the Virtual Cell computational environment (18) and translated into MATLAB code for use with an existing optimization algorithm (19). Computations were performed in parallel on 12 cores of a Linux cluster. Free model parameters (Table 1) were fit using a derivative-free optimization technique with a cost function based on a weighted sum of squares (WSS) residual between the experimental data (x i ; mean of replicates) and the modelcomputed values (m i ), weighted by the inverse of the standard error (S.E. i ); and accumulated for each experimental data point (i) over the number of data points (n j ) in each data set (j) (Equation 1).
The error for each data set was divided by the number of data points to give equal weight to each data set. All of the data were normalized so that the WSS could be pooled across experimental data sets. The best fit parameters were determined to be those giving the lowest total normalized error for each mechanistic model (Tables 2 and 3). For each mechanistic theory, the order of cAMP binding with the lower WSS error when compared with published experiments (15,16) is described in detail (all data are shown in Tables 2-4). These parameters (Para.) represent all potential degrees of freedom that could be fit in each model. The parameter names are consistent across all models. The regulatory mechanism theories were tested by constraining various factors. All biochemical interpretations are written as the effect of C-subunits on cAMP binding, but each could be written as the reverse as well. Unconst, unconstrained.
An R 2 value was calculated for the result of each model to determine how well each mechanistic model individually fits the given data set (Equation 2).
x is the mean of all the experimental data in a given data set. A reduced 2 statistic was also calculated ( 2 red ) as a measure of the goodness of fit of the model that also takes into account the number of degrees of freedom (DF) in the model and thereby favors simpler models with fewer adjustable parameters (Equation 3).
Finally, an F-test was performed to compare the fit of nested models to determine whether one gives a significantly better fit to the data than the other (20). For this we calculated an F statistic (Equation 4) where P is the broader model with more degrees of freedom and Q is the reduced model with fewer degrees of freedom. For this statistic, v1 ϭ DF Q and v2 ϭ N Ϫ (DF P ϩ 1). The significance level was set at 0.01.
In data sets IV and V, each measurement was reported as the amount of cAMP bound versus free cAMP. We accounted for the fact that both of these quantities were dependent variables in the model in calculating the error sum of squares. For data set V, because Ͼ97% of the total cAMP was unbound, the free cAMP acts like an independent variable. For data set IV, most of the cAMP was bound. Therefore, a hill equation was fit to data set IV, and the error between the model and the Hill equation at a given free cAMP concentration was used as the WSS error so that the error across all models could be calculated consistently.

Mechanistic Models of PKA
Dimer Model-The simplest model, the dimer model ( Fig.  2A), assumes that the R 2 C 2 holoenzyme can be treated as two independent R-C heterodimers with no cross-talk between them. This model assumes that both cAMP binding events can affect the activation of C-subunit. After being constrained thermodynamically, this model reduces to six states with five degrees of freedom: K D B , K D A , K D C , and two thermodynamic factors, F 1 and F 3 ( Table 1).
CBD-A-regulated Model-The CBD-A-regulated model ( Fig.  2B) assumes that although CBD-B must be bound before CBD-A only the binding of cAMP to CBD-A affects the activation of the C-subunit (9). The binding of cAMP to CBD-A is assumed to cause a conformational change in the R-subunit, activating the C-subunit. Moreover, because holoenzyme structures suggest that the activation of the first and second C-subunits might be dependent events, this possibility was permitted in the model (13). The CBD-Bs are considered independent of each other, but the CBD-As are considered dependent binders (13). The two R-subunits in the R 2 dimer were interchangeable, but for ease of naming, the first R-subunit to bind a C-subunit is R, whereas the first R-subunit to bind cAMP is RЈ (Fig. 1). Also, the C-subunit binding R is (C), whereas the

Model transitions and best fit values
The transition numbers correspond to the MM in Fig. 2 and denote the K D values of transition from one state to the next for the dimer, CBD-A-regulated, CBD-B-regulated, dually regulated, and unconstrained models. The functions relate the various parameters in Table 1 to the K D for a given transition. The K D values for each transition using the best fit parameters of each model are given. INF, transitions that are prevented to test a given biochemical mechanism; †, K D values that include a parameter that was not intentionally constrained, but the best fit set that parameter to be as arbitrarily high as was allowed by the model. C-subunit binding RЈ is CЈ. This model allows cAMP binding to CBD-AЈ on the RЈ-subunit to affect the activation of the C-subunit (21). In the R 2 C 1 transition state, the cAMP was assumed to bind to the free RЈ-subunit before binding to the R-subunit that is already bound by a C-subunit. After being constrained thermodynamically, the model reduced to 14 states with nine degrees of freedom: , and four thermodynamic factors relating the unleashing of the C-subunit and the binding of cAMP, F 3 , F 4 , F 7 , and F 9 .
CBD-B-regulated Model-Similar to the CBD-A-regulated model, this model also assumes that the CBD-Bs are independent of each other but that the CBD-As are dependent binders. However, unlike the CBD-A-regulated model, this model assumes that cAMP binding to CBD-B affects the unleashing of the bound C-subunit independently of the CBD-A being bound (Fig. 2C). This is consistent with the conformational selection model as CBD-A bound to the C-subunit cannot undergo the conformational transition that is considered to be a prerequisite for cAMP binding. After being constrained thermodynamically, this model had 12 states with nine degrees of freedom: , and four thermodynamic factors, F 2 , F 3 , F 5 , and F 9 .
Dually Regulated Model-The CBD-B-regulated model ( Fig.  2D) was expanded to allow both CBD-A and CBD-B binding to affect activation of the catalytic subunits. This added the possibility of forming a state with one R 2 dimer, one C-subunit, and four cAMPs (R 2 C 1 -cAMP 4 ) as well as a state with one R 2 dimer, two C-subunits, and two cAMPs (R 2 C 2 -cAMP 2 ). In this model, transitions 1-9, 11, 15-19, and 21 ( Fig. 2D) share the same properties as

Best fit parameter values and error
The best fit parameters for each model are shown. All K D value are in units of M, and all constants are unitless. The WSS error (Equation 1) and red 2 (Equation 3) are shown at the bottom as well as the S.D. and mean of the residuals for the five data sets. INF, transitions that are prevented to test a given biochemical mechanism; *, numbers were preconstrained to test various models; †, numbers were not intentionally constrained, but the best fit set them to be as arbitrarily high as was allowed by the model.

TABLE 4 K D and Hill coefficients
To compare cooperativity between experiments and models, the data sets and models after fitting were compared with a simple Hill equation. The resultant K D and Hill coefficients are shown for the experimental data (Exp., shown in boldface), dimer, CBD-B-regulated (CBD-B), CBD-A-regulated (CBD-A), dually regulated (Dually), and unconstrained (Uncon.) models. the CBD-B-regulated model. Transitions 12 and 20 describe the transfers to an R 2 C 2 -cAMP 2 state, and transitions 10 and 23 describe the transfers to an R 2 C 1 -cAMP 4 state. The model reduces to 14 states with 11 degrees of freedom:
Unconstrained Model-In an attempt to find the best model approximation to the data without assuming a specific regulatory mechanism a priori, an R 2 C 2 holoenzyme model was created to test the cooperativity of activation across the molecule (Fig. 2E). This model was developed without the constraints of either CBD-A-or CBD-B-regulated models. Each combination of R 2 dimer, C-subunit, and cAMP was available for binding with each reaction loop having its own thermodynamic reaction constant. The only constraint was that the CBD-B site must be bound before CBD-A. The model had 16 states with 14 degrees of freedom: , and nine thermodynamic factors F 1 -F 9 .

The Addition of Kemptide
Several studies have shown that the presence of kemptide affects PKA activation either through C-subunit buffering or  . Removing states to test cAMP regulation. One of two key states was removed from the dually regulated model, and then the new models were fit to the data to probe which cAMP binding event regulates C-subunit unleashing. A, resulting fits to data sets I-V from the dually regulated model (solid line) with the R 2 C 2 -cAMP 2 state removed (long dotted line) or with the R 2 C 1 -cAMP 4 removed (short dotted line). B, WSS error for each of the models without refitting when the indicated state was removed. Error bars are the S.E. of the experimental data set.
through binding to the C-subunit catalytic cleft and promoting activation (22,23). A simple model of kemptide buffering was created (Fig. 4A) and fit to experimental data with varying concentrations of kemptide (8,70, and 140 M) in the presence of excess ATP and Mg 2ϩ (23) to determine the buffering capabilities of kemptide on free C-subunit concentration (Fig. 4B). The kemptide competition coefficient (K eq ) relates the amount of free C-subunit, [C], available to bind the R 2 dimer to the amount of C bound by kemptide, [C_kemp].
The equilibrium constant is a function of the rate of phosphorylation of kemptide as well as the K D of interaction between C and phosphorylated and non-phosphorylated kemptide. The dually regulated model was recreated to include kemptide buffering of C-subunit and fit to the five data sets with data sets I-III now having 100 M kemptide. The resulting parameter values were then compared with those found in the absence of kemptide buffering to determine whether it affected any of the conclusions from the models.

Validating Models with Mutant PKA
The predictive value of the models was tested using four additional independent data sets from variants of PKA harboring mutations resulting in either a silent CBD-A site (RIG201E) or a silent CBD-B site (RIG325D) (16). The mutants were tested under the same conditions as data sets IV and V to quantify the percent saturation of PKA-RI␣ with increasing concentrations of cAMP. It was assumed that maximum saturation of the R-subunit was achieved with 50 M cAMP. The dimer, CBD-A, CBD-B, and dually regulated models were altered to silence the CBD-A or CBD-B binding site and were compared against the corresponding mutant data both in the absence of C-subunit and as under saturating C-subunit conditions. Making no other changes to the parameters estimated from the five wild type data sets, the mutant model results were compared with experimental measurements using the statistics described above.

RESULTS
Dimer Model-The first question is whether the two R-C heterodimers talk to each other. The dimer model, using the best fit parameters shown in Table 2, was only able to match three of the five data sets (II, IV, and V) with R 2 values greater than 0.95 (Fig. 2). Under the assumptions of this model, cAMP had subnanomolar affinity for the CBD-B site and nanomolar affinity for the CBD-A site on free R-subunits. cAMP had micromolar affinity for both binding sites in the R-C dimer. This model also showed an increase in K D between R-and C-subunits by 7 orders of magnitude after two cAMPs were bound.
The dimer model was unable to match the K D from data sets I and V to within 20% of the experimentally observed value. Data sets I-III show decreasing cooperativities with increasing concentrations of R-subunit. In contrast, increasing R-subunit concentration in the dimer model resulted in increasing cooperativity, suggesting that most of its cooperativity is based on cAMP binding to free R-subunit (Table 4). This model could not match data sets I and III closely but did fit data set II with an R 2 of 0.99. It was also able to reproduce the free R-subunit cAMP binding with no predicted value further than two S.E. from the experimental data, suggesting that the R-C heterodimer interactions necessary for data sets I, II, III, and V are less significant for binding to free R-subunit. The 2 red for this model was 4.2.
CBD-A-regulated Model-Next, the condition where only CBD-A binding affects C-subunit release was considered. The best fit parameters for the CBD-A model (Tables 2 and 3) led to a lower WSS error than that of the dimer model (4.33 versus 5.09; Fig. 2B). CBD-Bs for both the free R-subunit and the R 2 C 2 holoenzyme were bound with 214 nM affinity due to thermodynamic constraints. This is an order of magnitude higher than has been reported for binding to R 2 dimer (22 nM (9)) and an order of magnitude lower than binding to holoenzyme (2.9 M (16)). The first cAMP binding the CBD-A showed subpicomolar affinity when binding to the free R-subunit that is 3 orders of magnitude lower than previously reported (1-20 nM (9, 16)) even when the CBD-B site is already bound. These K D values imply that this is not a reasonable model for PKA activation despite its low total WSS error ( Table 2). More interestingly, cAMP bound CBD-A in the R 2 C 2 holoenzyme with a K D of 19.3 M, which leads to the R 2 C 2 -cAMP 2 state being nearly unpopulated under physiological conditions. Therefore, even in a model intended to be regulated by CBD-A, only CBD-B is filled before the C-subunit is unleashed. In this model, the K D for R-C interactions increased by 4 orders of magnitude after cAMP bound to the first CBD-A and by an additional 5 orders of magnitude when the second CBD-A was bound. This increase in dissociation of R 2 dimer and C-subunit because of the presence of cAMP is consistent with published reports (9).
The CBD-A-regulated model was able to fit data sets I, II, and V with strong correlations (R 2 ϭ 0.98 for all). However, the model showed a distinct biphasic slope for the binding of cAMP to free R 2 dimer in data set IV because of the difference in K D between binding CBD-B and CBD-A. This is inconsistent with experimental findings and resulted in a slightly weaker correlation (R 2 ϭ 0.97). The high amount of CBD-B bound in the holoenzyme at low cAMP concentration was also somewhat inconsistent with data set V, but overall this model had the second lowest WSS error but the third lowest 2 red (4.36) due to the increase in degrees of freedom compared with the dimer model. CBD-B-regulated Model-The reverse was then considered where only CBD-B binding was allowed to affect C-subunit release. The CBD-B-regulated model was able to fit two of the five data sets (II and IV) as well as or better than the CBD-Aregulated model (0.98 and 1.00 versus 0.98 and 0.97, respectively; Fig. 2, B and C). However, data sets III and V had the worst fits of any of the models (R 2 of 0.88 and 0.94, respectively).
The best fit parameters (Tables 2 and 3) showed 1 nM affinity for the cAMP binding to CBD-B and CBD-A on free R-subunit; this is consistent with literature values (9). This model showed an increase of 3 orders of magnitude in K D for cAMP binding to the CBD-B in the holoenzyme compared with the R 2 dimer. This is consistent with a previous report (16) and allowed this model to fit both data sets II and IV very well (R 2 ϭ 0.98 and 1.00). To fit the activation curves for data sets I and II, this model showed a decrease in affinity between the R 2 dimer and C-subunits by 7 orders of magnitude due to cAMP binding. The binding to the first CBD-B site decreased the affinity between the R 2 dimer and C-subunits by 3 orders of magnitude and an additional 4 orders of magnitude due to the second CBD-B site. However, although this matches data set V at low concentrations of cAMP, at high concentrations, cAMP was prohibited from binding due to the excess C-subunit that keeps the R-subunit bound to the C-subunit. This effect reduced the strength of the agreement between the data and model (R 2 ϭ 0.94).
The one data set that could not be fit with R 2 Ͼ 0.9 was III. The experiments showed that under saturating cAMP concentrations with more R-subunit than C-subunit PKA was not fully activated. This result is consistent with previous work (22) and suggests an equilibrium between active and inactive PKA under saturating cAMP conditions. However, the CBD-B-regulated model predicted that there should be 100% PKA activation under saturating cAMP concentrations regardless of R-subunit concentration because the R 2 dimers cannot have both four cAMPs bound and a C-subunit. This is contrary to observation, resulting in a weaker correlation (R 2 ϭ 0.88) and a higher 2 red value of 10.41 between the model and data than either the CBD-A-regulated Model or the dimer model.

CBD-B and CBD-A Dually Regulated
Model-Next, a model was considered that allowed CBD-A and CBD-B to affect C-subunit release. This model had the lowest WSS error (3.02) and was the best fit model using the 2 red statistic (Tables 2 and 3). The nanomolar affinities for cAMP binding to the CBD-B of free R-subunit are consistent with a previous report (9). cAMP binding to the CBD-A of the R 2 dimer had 1 nM affinity, which is again consistent with literature (9, 16) and does not produce the same biphasic curve as the CBD-A-regulated model (R 2 ϭ 1.00). The model also showed an increase of 3 orders of magnitude in the K D for binding of cAMP to the holoenzyme compared with free R-subunit. The dissociation constant of cAMP binding to the holoenzyme is 0.762 and 2.91 M for the first and second cAMP binding, respectively. This is consistent in magnitude with previously reported values (16) and is in keeping with the activation curves of data sets I, II, and III with R 2 values of 0.98 for each.
The affinity between the R 2 dimer and the C-subunits in the model was decreased by 6 orders of magnitude due to cAMP binding. The existence of an R 2 C 1 -cAMP 4 state allowed the model to approach 80% activation asymptotically under satu-rating cAMP conditions with excess R 2 dimer compared with C-subunit, which is consistent with data set III. This addition to the CBD-B-regulated model also promoted a secondary pathway to reach an R 2 -cAMP 4 -bound state through the R 2 C 1 -cAMP 4 state. The added degree of freedom allows the R 2 C 1 -cAMP 4 -bound state that was necessary to accurately fit data sets III and V without creating the biphasic fit to data set IV. The model fit all five data sets well with R 2 values greater than 0.98 and a 2 red value of 3.02. The R 2 C 2 -cAMP 2 and R 2 C 1 -cAMP 4 states were individually removed after fitting the data to determine which CBD regulated the first and second C-subunit activations. The removal of the R 2 C 2 -cAMP 2 state had no effect on the fits for the first four data sets and a negligible effect on the fifth, suggesting that this state does not play a significant role in the steady-state activation of PKA (Fig. 3). This result also suggests that CBD-B regulates the unleashing of the first C-subunit, which has been consistent through all the models (Fig. 3). Similarly, the R 2 C 1 -cAMP 4 state was removed to determine which CBD regulates the second C-subunit activation. Fits for data sets I, II, and IV were mostly unchanged, but the model was unable to match the binding of cAMP at the high concentrations shown by data set V. Furthermore, this model showed 100% activation at high cAMP concentrations when compared with data set III. Together this suggests that the unleashing of the first C-subunit is dominated by cAMP binding to CBD-B, whereas the unleashing of the second C-subunit requires the fourth cAMP to bind CBD-A.
Unconstrained Model with Additional Degrees of Freedom-Finally, to test that the dually regulated model had the necessary and sufficient degrees of freedom, a last model was developed. This model included two additional states compared with the dually regulated model, R 2 C 2 -cAMP 3 and R 2 C 2 -cAMP 4 , as well as the possibility of CBD-BЈ binding affecting the dissociation of C-subunit (Fig. 2E). However, the best fit computed free constants F 6 and F 8 to be as high as the model constraints would allow. This resulted in negligible occupancy of these additional states. Moreover, CBD-BЈ did not affect the C-subunit (F 1 ϭ 1) even when it was not constrained. These results reduce the unconstrained model to the form of the dually regulated model, suggesting that these additional degrees of freedom are not needed to explain the current data (Tables 2 and 3). The best fit values describing state transitions for the unconstrained model were within the error of the corresponding transitions of the dually regulated model including all of the dissociation and thermodynamic constants, but the unconstrained model has a higher 2 red value of 3.63 due to the additional degrees of freedom.
Adding Kemptide-The best fit for the kemptide competition coefficient was a K eq value of 61.9 M for the interaction (Fig. 4). This addition was then added to the dually regulated model and fit to the same five data sets with data sets I-III also including the kemptide buffering. The resultant fit was then compared with the fit in the absence of kemptide buffering.
This model resulted in a comparable fit to the model in the absence of kemptide (WSS ϭ 3.07). The fit parameters were all changed by less than a factor of 2 (Table 5) except for K D C , which decreased by 60% to compensate for kemptide buffering. The mechanism of cAMP binding remained unchanged by the addition. This suggests that the conclusions drawn from the previ-ous models were not affected by the inclusion or exclusion of kemptide buffering, which predominately effects the association between R-and C-subunits.
Model Comparison-The dimer model fit two of the data sets well, but the other three all had R 2 values less than 0.98. The discrepancy was mainly due to a difference in the cooperativity between the dimer model and the experimental data. This led to a WSS error of 5.09 and a 2 red value of 4.2. The CBD-A-regulated model was able to fit three of the five data sets and showed an improvement in accuracy (WSS ϭ 4.33) but a worse goodness of fit ( 2 red ϭ 4.36) because of the greater number of degrees of freedom. Moreover, the best fit required the K D for cAMP binding CBD-B on free R-subunit to be over 100 nM, which is much higher than has been experimentally reported (9), suggesting that despite its comparatively low error sum of squares it is not a likely model for PKA activation.
The CBD-B-regulated model was able to fit two of the five data sets (Table 5) but did not match several trends including percent activation under saturating cAMP conditions in data set III. This led to the highest sum of squared error (12.55) and a worse goodness of fit ( 2 red ϭ 10.41). Therefore, this model was also unable to explain the observed activation of cAMP binding of PKA. The dually regulated model alone was able to fit all five of the data sets with an R 2 value greater than 0.98. This led to the lowest WSS error (3.02) and a 2 red value of 3.16. This model therefore achieved a much better fit than any of the alternative models without overfitting the data. To test significance, an F-test was run to compare nested models. This returned p values of 0.0038 and 2.6 ϫ 10 Ϫ10 for the CBD-A-and CBD-Bregulated models, respectively, when compared against the dually regulated model. Both of these fall under a 0.01 confidence that these are statistically different. To test for a necessary and sufficient model of PKA activation, a final unconstrained model was created with 14 degrees of freedom. However, even when given these extra degrees of freedom, the best fit reduced down to the dually regulated model (with only minor alterations) with a WSS error of 3.02 and a 2 red value of 3.63. The higher 2 red value is the result of the additional degrees of freedom in the unconstrained model with no improvement in fit. The F-test returned a p value of 0.999, which prevents the rejection of the null hypothesis that they are the same. Together this suggests that the dually regulated model is the simplest faithful representation of currently available data.
Validating Models with Mutant PKA Data-To test the predictive ability of the models, the dimer, CBD-A-regulated, CBD-B-regulated, and dually regulated models were modified only by removing specific states (with no changes to the previously estimated optimal parameters) to mimic the CBD-A-and CBD-B-inhibited mutants, RIG201E and RIG325D. The results of these mutant models were compared with published measurements using the RIG201E and RIG325D mutants of cAMP binding to free R 2 dimer and R 2 C 2 holoenzyme under conditions similar to those of data sets IV and V. The mutant MM models are shown in Fig. 5.
The dimer model matched the cAMP binding to the CBD-Binhibited mutant well (R 2 ϭ 0.99 for both free R-subunit and R 2 C 2 holoenzyme) but matched the CBD-A-inhibited mutant poorly (R 2 ϭ 0.88 for free R-subunit and 0.02 for R 2 C 2 holoenzyme; Fig. 6A). Combined, the WSS error for the mutant data sets was 15.95 and did not fit the experimental results.
The CBD-A-regulated model had the second lowest WSS error for fitting the original five data sets but the highest WSS error when compared with data from mutants. This model could not match any of the validation data sets with an R 2 value greater than 0.9, and two of the comparisons even resulted in negative R 2 values (Fig. 6B). The lowest R 2 values were for binding to the free R-subunits. The extraordinarily high K D for cAMP binding to CBD-B on free R-subunit that is necessary for the CBD-A-regulated model is inconsistent with previously reported values and with the data acquired with the mutants. Combined, the WSS error for the validation data sets was 339.4.
The CBD-B-regulated model matched the binding of cAMP to free R-subunit for both mutants with R 2 values greater than 0.9. However, the model for both of the R 2 C 2 holoenzyme data sets significantly underpredicted cAMP binding, resulting in R 2 values of 0.01 for the CBD-A-inhibited mutant and Ϫ0.27 for the CBD-B-inhibited mutant (Fig. 6C). In total, the WSS error for these four validation data sets was 46.2.
The dually regulated model, which was the best fit for the original five data sets, also matched the validation data the best. Three of the four validation data sets were matched with an R 2 value of 0.99, and all four had a combined WSS error of 5.87 (Fig. 6D). The model matched the general trend of cAMP binding to the R 2 C 2 holoenzyme with inhibited CBD-A but predicted lower binding, resulting in an R 2 value of 0.62. Because this model was predictive for cAMP binding to free R 2 dimer but not to R 2 C 2 holoenzyme for the CBD-A-inhibited mutant, there might be a CBD-B interaction in the holoenzyme that is not taken into account in the current model.

DISCUSSION
Despite several decades of study, the mechanism by which PKA is activated by cAMP is still unclear. Recently, even the previous canon that C-subunits are released prior to activation has come into question, leaving competing theories about the events that lead to activation (24,25). If the conformational   OCTOBER 24, 2014 • VOLUME 289 • NUMBER 43 selection mechanism of activation is applicable to PKA, then to bind cAMP, both CBD-A and CBD-B have to be able to explore their conformational states. This is feasible in the case of CBD-Bs as these domains are exposed to solvent, and the interface between CBD-Bs and C-subunits is relatively small (26). In contrast, CBD-As have a large, mostly hydrophobic interface with the C-subunits that makes them very stable. Thus, if the conformational selection allosteric mechanism is relevant for PKA activation, then CBD-As can bind cAMP only after separation from the C-subunit (CBD-B-regulated model). Conversely, if a non-conformational allosteric mechanism drives PKA activation, then cAMP can bind CBD-As before separation from the C-subunit (CBD-A-regulated model).

Markov State Models in Protein Kinase A Activation
Our data show that the assumption that cAMP binding to CBD-B does not directly affect the activation of the C-subunit is inconsistent with the current empirical evidence. The CBD-Aregulated model was unable to replicate cAMP binding to the R 2 C 2 holoenzyme and required parameters that were outside the range of previously reported cAMP binding events. Furthermore, when any of the models were given the option of allowing CBD-A to be bound to the R 2 C 2 holoenzyme, the best fits set the K D of cAMP binding to CBD-A so high that this state would not be populated under physiological concentrations of cAMP until after the C-subunit is unleashed (Table 2). This strongly suggests that the binding of cAMP to the first CBD-B triggers the unleashing of the first C-subunit prior to the binding of cAMP to CBD-A, ruling out an induced fit mechanism. This result is consistent with the conformational selection mechanism that recently was reported in single CBD-A studies (8,27).
When a CBD-B-regulated model (conformational selection mechanism) is considered, the data show that CBD-B regulation alone cannot reproduce the equilibrium between bound and free C-subunits under saturating cAMP concentrations (data set III). Instead, a combination of CBD-B and CBD-A regulation most accurately represented the empirical data. The first C-subunit activation is regulated by binding the CBD-B according to a conformational selection mechanism; the second C-subunit activation is better described by a CBD-A-dominated activation similar to a non-conformational selection mechanism. It is possible that certain conformations are available with only one C-subunit present that are not favorable for the R 2 C 2 holoenzyme.
The dimer model representation of PKA-RI␣ was unable to reproduce the cooperativity of activation in data sets I-III and V, highlighting the role of R-C heterodimer interactions. Also, the interactions between the CBD-As both in the free R-subunit and in the R-C complex were necessary for the formation and fitting of the best model. This is consistent with recent publications (26,28). However, the diverse R 2 C 2 structures that various PKA isoforms exhibit suggest that although this model is consistent with current PKA-RI␣ experiments it would be expected that RI␤, RII␣, and RII␤ will exhibit significantly different mechanistic activations possibly with less interaction between the R-C heterodimers (13,29,30).
When the models were used to predict data that had not been used for parameter estimation, several of the models with comparatively low fitting errors were shown to have little to no predictive power. However, the dually regulated model had the lowest fitting error and very high predictive power. The mutant validation highlighted the unrealistic parameters in the CBD-A-regulated model in particular and to a lesser extent in the dimer model. Again, the dually regulated model was the only model able to predict the binding of cAMP to both the CBD-A and CBD-Bs in the free R-subunit and in the R 2 C 2 holoenzyme conformation. This suggests that not only does the model match broad PKA characteristics but also accurately predicts CBD interactions caused by mutations. Only the CBD-A-inhibited holoenzyme was underpredicted in the validation data. This suggests that there may be a CBD-B interaction that has not previously been considered.
The dually regulated model was also compared with other literature data (data not shown). This model accurately predicted PKA dissociation under saturating cAMP presented in Vigil et al. (22). Additionally, Herberg et al. (9) reported that cAMP activated PKA-RI␣ with a K D of 166 nM when there were equal concentrations of R-and C-subunits, whereas the dually regulated model predicted a slightly higher K D of 250 nM under the same conditions. This moderate difference could be explained by variations between experimental designs between laboratories. For example, the activity assays of Døskeland and co-workers were run at 37°C for 40 min (15), whereas the assays of Herberg et al. (9) were run at room temperature for 2 min.
Limitations-One surprising result from the dually regulated model was that the best fit for F 9 , which relates the order in which cAMP binds to free R-subunit, was ϳ2.0. When F 9 was fixed at any other value, the fit to data set IV and the validation of the CBD-A-inhibited mutant was significantly diminished (data not shown). This demonstrates that the effective K D for cAMP binding to CBD-A when one CBD-B is bound is the same as the effective K D for cAMP binding CBD-A when both CBD-Bs are bound to the free R 2 dimer. This suggests that for the free R 2 dimer, the first CBD-A might be independent of CBD-B binding, which experimental results suggest is possible (14). In these models, the assumption was that CBD-B always binds first. However, with this result for free R 2 dimer, it is also possible that CBD-A binding first to free R-subunit should be included in future work.
Furthermore, the dually regulated model could not differentiate between release of the C-or CЈ-subunit. Therefore, the model cannot rule out the possibility that the dominant effect of binding CBD-B was release of the CЈ-subunit. However, this could be tested by the creation of an R 2 heterodimer mutant with one R-subunit having both CBD-A and CBD-B silent and one WT R-subunit.
One can speculate based on experimental evidence (8, 27) that even though the release of the second C-subunit is after the binding of CBD-A by cAMP it is possible that it is a non-standard form of conformational selection. Future experiments to determine the interaction interface between the R 2 dimer and the remaining C-subunit when four cAMPs are bound could help explain this mechanism.
Finally, recent work has highlighted the role direct binding between protein kinase A-anchoring proteins and phosphodiesterases plays in PKA activation in vivo (32,33). Therefore, the current model, which is based on in vitro data, cannot fully capture in vivo activity.
Conclusion-When targeting inactive states of PKA, one has to consider that the CBD-B is in a dynamic equilibrium between two conformations, whereas CBD-A is locked in one conformation and is inaccessible to cAMP until the C-subunit is unleashed. Previous results imply that a single mutation inside CBD-B can cause a significant shift in the equilibrium that leads to increased stability of the R 2 C 2 holoenzyme (31). This work illustrates how Markov model simulations can be used for understanding the fine molecular mechanism of cAMP binding in other CBD-containing proteins.