A Computational Modeling and Simulation Approach to Investigate Mechanisms of Subcellular cAMP Compartmentation

Subcellular compartmentation of the ubiquitous second messenger cAMP has been widely proposed as a mechanism to explain unique receptor-dependent functional responses. How exactly compartmentation is achieved, however, has remained a mystery for more than 40 years. In this study, we developed computational and mathematical models to represent a subcellular sarcomeric space in a cardiac myocyte with varying detail. We then used these models to predict the contributions of various mechanisms that establish subcellular cAMP microdomains. We used the models to test the hypothesis that phosphodiesterases act as functional barriers to diffusion, creating discrete cAMP signaling domains. We also used the models to predict the effect of a range of experimentally measured diffusion rates on cAMP compartmentation. Finally, we modeled the anatomical structures in a cardiac myocyte diad, to predict the effects of anatomical diffusion barriers on cAMP compartmentation. When we incorporated experimentally informed model parameters to reconstruct an in silico subcellular sarcomeric space with spatially distinct cAMP production sites linked to caveloar domains, the models predict that under realistic conditions phosphodiesterases alone were insufficient to generate significant cAMP gradients. This prediction persisted even when combined with slow cAMP diffusion. When we additionally considered the effects of anatomic barriers to diffusion that are expected in the cardiac myocyte dyadic space, cAMP compartmentation did occur, but only when diffusion was slow. Our model simulations suggest that additional mechanisms likely contribute to cAMP gradients occurring in submicroscopic domains. The difference between the physiological and pathological effects resulting from the production of cAMP may be a function of appropriate compartmentation of cAMP signaling. Therefore, understanding the contribution of factors that are responsible for coordinating the spatial and temporal distribution of cAMP at the subcellular level could be important for developing new strategies for the prevention or treatment of unfavorable responses associated with different disease states.


Introduction
For nearly 40 years, subcellular compartmentation has been offered as an explanation for how cAMP, the ubiquitous and diffusible second messenger, can both regulate a multitude of cellular functions and elicit specific and selective responses. Despite widespread recognition of the importance of cAMP compartmentation in tightly controlling local signaling, exactly how compartmentation occurs is still poorly understood. The general definition of compartmentation in this context is when a gradient exists in the concentration of cAMP between two locations. As it relates to cell signaling, the concentration gradient is relevant when it affects the potential for cAMP to activate an effector, such as protein kinase A (PKA), in one location but not another. A number of processes have been suggested to contribute to this phenomenon, but studies have offered conflicting data that differ in their interpretation and assessment of key players.
Localized degradation by phosphodiesterases (PDEs) has been a prime focus of many studies attempting to understand the basis of cAMP compartmentation [1][2][3][4][5]. Phosphodiesterases are thought to contribute to the generation of cytosolic cAMP gradients either by acting as functional barriers to diffusion that result in lower levels of cAMP distal to its site of production or as sinks that deplete cAMP in localized areas. Evidence clearly demonstrates that PDE activity is an essential factor in cAMP compartmentation. This has been illustrated by employing a number of different experimental approaches, including Jurevcius and Fischmeister who used patch clamp electrophysiology to demonstrate that in frog ventricular myocytes, inhibition of PDE activity allows local stimulation of cAMP by β-adrenergic receptors to enhance distal Ca 2+ channel activity [6]. On the other hand, Zaccolo et al. used a genetically encoded FRET-based biosensor to demonstrate that β adrenergic stimulation elicits a localized pattern of cAMP production in neonatal cardiac myocytes that is disrupted by inhibition of PDE activity [7]. However, the question of whether or not PDE activity alone is sufficient to explain the compartmentalized behavior of cAMP signaling is still debated.
Computational modeling has proven to be a useful tool in investigating the relative contribution of PDEs to cAMP compartmentation [8]. All modeling studies support the idea that PDE activity is necessary for cAMP compartmentation. At least one study has predicted that it is theoretically possible for artificially high levels of PDE activity alone to explain compartmentation [9]. Other models using more realistic levels of PDE activity suggest that factors such as the shape of the cell and the rate of cAMP diffusion play critical roles in explaining the existence of cAMP gradients within a cell [10][11][12][13][14][15][16][17][18][19][20][21].
One way that cell shape may be a factor in compartmentalizing cAMP signaling is by affecting the surface-to-volume ratio. Studies using FRET-based biosensors in neurons have found cAMP levels to be higher in dendrites than cell bodies [22]. It was suggested that this could be due to the higher surface-to-volume ratio found in dendrites, resulting in greater membrane bound adenylyl cyclase activity and reduced cytosolic PDE activity. Subsequent modeling supported the feasibility of this hypothesis, without having to assume the involvement of other factors [18]. Feinstein et al. found that the surface-to-volume ratio of a cell can contribute to generation of cAMP gradients, but it was also necessary to assume that the movement of cAMP is slower than the rate of free diffusion [20]. Other models have been able to explain compartmentation independent of cell morphology, as long as it was assumed that cAMP diffusion is somehow restricted [10, 11, 13-15, 17, 23].
Although the potential effect that the surface-to-volume ratio of a cell has on cAMP compartmentation has been examined, the influence of the actual size and shape of subcellular compartments is less well understood. A major reason is that the physical nature of these microdomains is not well described. Previous modeling studies often circumvented this issue by using loosely defined membrane and cytosolic domains and treating the movement of cAMP between them as fluxes that do not require knowledge of the number, size, or location of these compartments.
Previous experimental studies have shown that receptors associated with cholesterol rich lipid rafts, which include caveolae, can elicit cAMP responses that are distinctly different from those produced by extracaveolar receptors found outside of lipid rafts [16,24,25]. Lipid rafts are liquid-ordered domains of the membrane rich in cholesterol and sphingolipids. Caveolae are a specific subset of lipid rafts that contain caveolins, proteins involved in the formation of signaling complexes that include β 1 and β 2 adrenergic receptors (βARs) as well as adenylyl cyclase isoforms 5 and 6 (AC5/6) [26][27][28][29]. In cardiac myocytes, activation of receptors associated with caveolar lipid rafts are involved in local cAMP production and PKA-dependent regulation of L-type Ca 2+ channel function [24,28,29]. There is evidence that these types of compartmentalized cAMP responses also occur in the transverse tubules (t-tubules) of cardiac myocytes [30]. T-tubules are invaginations of the plasma membrane that come in close proximity to the junctional sarcoplasmic reticulum (SR) forming dyadic junctions [31]. Therefore, it is possible that the size, shape, and distribution of caveolae, especially those associated with the restricted space at cardiac dyadic junctions, may contribute to compartmentation of cAMP signaling in cardiac myocytes.
The purpose of the present study was to apply novel computational approaches to predict whether PDE activity alone or in conjunction with restricted diffusion is sufficient to produce cAMP gradients in submicroscopic signaling domains.

Results
Experimental studies in cardiac myocytes have shown that activation of βARs associated with caveolar regions of the plasma membrane produce unique compartmentalized cAMP responses [24,25,28,29,32]. Other studies have used computational approaches to investigate the effect that cell morphology has on the generation of cytosolic cAMP gradients [18,20,23]. However, the importance that the organization and structure of submicroscopic signaling domains has on creating compartmentalized cAMP responses has not been addressed. To investigate how cAMP-mediated responses are localized and prevent initiation of global responses, we developed an idealized, partial differential equation (partialDE)-based 2D continuum model of a cardiac myocyte subspace with spatially distinct cAMP microdomains to allow for simulation of cAMP compartmentation and diffusion.
A schematic diagram of the longitudinal cross-section of an adult ventricular myocyte with dimensions of 100 μm x 20 μm [14] is shown in Fig 1A. This illustrates the repeating pattern of the sarcomeres, which are spaced 2 μm apart [33]. Fig 1B illustrates the 2D continuum model that we constructed to represent the subcellular sarcomeric space used in the simulations shown in Fig 2. The model represents the intracellular space between adjacent t-tubules. The ttubules are lined by caveolar domains spaced 100 nm apart [34,35]. Each unit is half the width of the cell (10 μm). Fig 1C shows a magnified section of the subsarcolemmal space consisting of a single caveolar domain (blue), which is 0.1 μm x 0.01 μm [34,36,37], and an adjacent 0.1 μm x 0.01 μm extracaveolar space (green).
β-adrenergic receptors (βARs) and adenylyl cyclases (ACs) (the site of cAMP production) were distributed equally among the 50 caveolar domains on each side of the sarcomeric space. Using this model we would expect to be able to readily track compartmentation as cAMP gradients. We define a "significant" gradient as one in which the concentration of cAMP drops by more than 15% of its value relative to the site of production. Furthermore, we identify gradients as being relevant to compartmentalized signaling if cAMP concentrations in one compartment reach levels likely to produce PKA activation (>1 μM).
Predictions from simulations using this 2D continuum partialDE model to simulate cAMP diffusion are shown as snapshots of cAMP concentration at different points in time across the microdomain space in Fig 2. The time course of spatial changes in cAMP concentration resulting from activation of βARs is shown when cAMP was allowed to move at a rate approximating free-diffusion (300 μm 2 /s) under conditions where no phosphodiesterases (PDEs) were present (Fig 2A), where PDEs were localized to the caveolar domain at concentrations consistent with those reported experimentally [14,15] (Fig 2B), and where PDE concentrations in the caveolar domain were increased 10-fold ( Fig 2C). Only miniscule gradients (sub-nanomolar) were observed during the 2.0 second simulations, even when the concentration of PDE was increased 10-fold. The prediction of the model led to no indicators of significant compartmentation, which we would have expected to observe as gradients of cAMP concentration within the subcellular sarcomeric space. Rather, the monochromatic color maps in Several studies [10-17, 20, 23, 38] have suggested that for compartmentation to occur, diffusion of cAMP must be substantially slower than the reported value of free diffusion in a dilute aqueous environment, which is 300 to 400 μm 2 /s [39,40]. Assuming that cAMP movement is affected by factors such as cytoplasmic viscosity and molecular crowding, the diffusion coefficient of molecules the size of cAMP has been estimated at 60 μm 2 /s [41]. However, slowing diffusion in our simulation was still insufficient to generate a spatial gradient, even when the amount of PDE activity was increased 10-fold above the levels believed to exist in cardiac myocytes ( Fig 2D). It has also been suggested that buffering of cAMP through its interactions with PKA can decrease the effective diffusion coefficient for cAMP even further, to values closer to 10 μm 2 /s [41]. Yet, even this marked reduction of cAMP diffusion rate in the simulation was insufficient to generate spatial gradients of cAMP ( Fig 2E).
In the 2D continuum model, PDEs were contained within the thin caveloar domain (i.e. all PDE is effectively along the plasma membrane). In order to more specifically address the question of whether PDEs can form a "functional barrier" to cAMP diffusion, we developed a 3D stochastic model of cAMP diffusion in a subcellular microdomain (depicted in Fig 1D) and implemented the model using MCell [42]. This approach allows for investigation of the  51 where βARs and AC5/6 are localized. Caveolar domains (blue) are 100 nm x 10 nm, spaced 100 nm apart. In this example, extracaveolar domains (green) are associated with the subsarcolemmal space of the t-tubules (between caveolar domains) as well as the peripheral sarcolemma. (D) A single caveolar domain and half of each adjacent extracaveolar flanking region. All of the PDE modelcules are located at a distance (L*) from the plasma membrane (the site of cAMP production). L indicates the most distal site from the plasma membrane in the compartment. The Mcell simulations were carried out in a subspace compartment (from L* to L) on this microdomain. (E) A schematic of the steady state distribution of cAMP along the microdomain in D as derived using the 1-dimensional continuum model. The concentration is A+B at the cAMP production site (z = 0). Beyond the PDE barrier (from L* to L), the concentration of cAMP is B.
contribution of spatial localization of microdomain specific signaling components (e.g., PDEs) to compartmentation. Evidence exists that the localization of signaling complexes is important in producing compartmentalized responses [28,43]  We then evaluated the effect of increasing the number of PDE molecules in the functional barrier by several orders of magnitude (Fig 3B-3E). Only when the number of PDE molecules was increased above 10,000 (Fig 3D and 3E) did a cAMP gradient become visible. This is illustrated most clearly by the accumulated concentration map at the bottom of each panel.
The results described above indicate that PDE activity alone is unlikely to produce significant cAMP gradients by acting as a functional barrier when cAMP was allowed to diffuse freely. We next tested if this was also the case when the rate of cAMP diffusion was decreased to 60 μm 2 /s, as shown in S3 Fig. This condition reflects the experimentally measured diffusion coefficient of cAMP like molecules that was determined by using fluorescein and the ϕ450 fluorophore, fluorescent molecules about the same size as cAMP that do not bind to PKA. In water, these molecules exhibit rates of free diffusion of~300 μm 2 /s, but inside cardiac myocytes the diffusion coefficient decreases to~60 μm 2 /s, attributable to collision with other macromolecules in the intracellular environment due to molecular crowding [44]. Despite the slower rate of diffusion, a functional barrier consisting of 10 PDE molecules was still not sufficient to produce a cAMP gradient (S3 Fig). In the setting of slower diffusion, it was necessary to increase the number of PDE molecules to at least 1000 (S3C-S3E Fig) before a small gradient was visible. Slowing the rate of cAMP diffusion also increased the concentration of cAMP observed at all levels of PDE activity.
We then repeated the simulations using a diffusion coefficient of 10 μm 2 /s (Fig 4). This reflects the further slowing of cAMP diffusion due to the effects PKA buffering as suggested experimentally [41]. Interestingly, under these conditions, there is evidence for cAMP compartmentation when the number of PDE molecules in the barrier is at least 100, which corresponds to a concentration of 41.5 μM (Fig 4B-4E).
Even with a diffusion coefficient of 10 μm 2 /s, the diffusion length (2 ffiffiffiffiffi Dt p ) of cAMP on relevant time scales (1-10 seconds) is much larger than the length scale of the caveolar domain (0.01-0.2μm). This leads to a nearly uniform concentration of cAMP in planes parallel to the plasma membrane, and therefore cAMP dynamics can be well-predicted by a 1D continuum model that can be solved analytically to obtain an expression for the steady state concentration of cAMP along the microdomain (see Methods and S1 Appendix). (Fig 1E shows  Thus far, our modeling results suggest that under physiologically relevant conditions, cAMP diffusion is not sufficiently restricted by the presence of PDE molecules to explain compartmentation. To further test the effects of the model parameter values on cAMP compartmentation, we used the 1D continuum model to perform a parameter sensitivity analysis. The black curves in Fig 5A shows the dependence of the cAMP compartmentation ratio R on PDE concentration for default parameters and a diffusion rate of 300 μm 2 /s as used in Fig 3. The red, blue, green, and cyan lines illustrate the sensitivity of R to changes in D/(k f L Ã ) (a parameter accounting for diffusion, the location of the PDE boundary, and the rate of cAMP association with PDE), k b (rate of cAMP dissociation from PDE), k cat (PDE catalysis rate), and J B (cAMP production rate), respectively. Each parameter was adjusted by ±20%, and the results are plotted as a pair of colored lines for each parameter change. Fig 5B and 5C show the sensitivity to the parameters for cAMP diffusion constants of 60 and 10 μm 2 /s, respectively, as in S3  Fig and Fig 4. For all cases, the cyan lines are almost entirely obscured by the black lines, indicating that the cAMP compartmentation is insensitive to changes in the cAMP production rate, J B . The compartmentation ratio R was most sensitive to D/(k f L Ã ) and k cat , however no change in R greater than 0.06 was observed. It is notable that the highest sensitivity of R (in terms of change of the absolute magnitude of the ratio) occurred in the ranges of PDE concentrations that are well above the physiologically relevant range. For PDE concentrations between 1 and 100 μM, R was insensitive to perturbations to all other parameters.
The simulations conducted thus far used models incorporating an idealized view of the 3D space between t-tubules in a cardiac myocyte. None of them contained realistic subcellular structures that might, in an actual cell, act as physical barriers to diffusion of cAMP. The cytosolic compartment of a cardiac myocyte is structurally complex and the site of cAMP production in t-tubules likely occurs in close proximity to the junctional SR, forming dyadic clefts. Movement of cAMP out of this space is also likely to be affected by the presence of mitochondria, which make up approximately 30% of the cardiac myocytes volume and are tightly packed around these structures. To examine the possibility that cAMP compartmentation might be observed in this type of restricted space, we created a 3D continuum anatomical barrier model using cryo-TEM images of adult mouse cardiac myocytes, as described  previously [45] (Fig 6A). The dimensions of the resulting dyadic cleft were approximately 1040 x 765 x 415 nm. Production of cAMP was generated by 15 AC molecules situated in the center of the dyadic cleft. These were surrounded by a hollow sphere of PDEs 25 nm thick and 200 nm in diameter. The t-tubules, SR, and mitochondria were assumed to be impenetrable barriers to direct diffusion of cAMP throughout the cytosol. We then examined the effects of varying PDE activity, as well as the diffusion coefficient, on cAMP gradients (Fig 7).   When cAMP was produced at the same rate as in the previous simulations (120 molecules/ s), and it was allowed to diffuse at a rate of 200 μm 2 /s, which is consistent with previous estimates of the diffusion coefficient in intact cells [22,46,47], no evidence of a gradient was observed when the number of PDE molecules surrounding the site of production was varied between 100 and 1000 (Fig 7). This behavior did not change when the diffusion coefficient for cAMP was reduced to 60 μm 2 /s. However, when the diffusion coefficient was further reduced to 10 μm 2 /s, a significant gradient was observed at all PDE concentrations.

Discussion
In this study, we developed several modeling approaches that included discrete cAMP microdomains to carry out simulations to investigate the importance of parameters proposed to be critical to compartmentation, including location and concentration of PDE activity, rates of cAMP diffusion, and anatomical barriers to diffusion.
There is a large body of literature demonstrating that the non-uniform distribution of PDE activity in different subcellular compartments plays a critical role in compartmentation of cAMP dependent responses [5]. This is achieved by various mechanisms, including interactions with A kinase anchoring proteins (AKAPs), which create signaling complexes that include not only PKA, but also adenylyl cyclase [48]. In the present study, we used a 2-dimensional continuum model to evaluate the effect of localizing PDE activity near the site of cAMP production. The results demonstrate that concentration of PDE activity in this location did not prevent cAMP diffusion throughout the space representing the area between adjacent ttubules.
Furthermore, a 10-fold increase in PDE activity was still insufficient to restrict cAMP diffusion, rather it reduced cAMP levels throughout the entire compartment. In previous modeling studies using similar levels of PDE activity, cAMP gradients could be generated if was assumed that the movement of cAMP was restricted by some mechanism [10,14,15,20,23,49]. However, we found that reducing the cAMP free diffusion coefficient to levels similar to those recently described in intact cardiac myocytes [41], did not promote compartmentation. PDEs were also not enough to prevent the experimentally observed higher concentrations of cytoplasmic cAMP from immediately diffusing into the caveloar spaces, where biosensors have suggested 10-fold lower cAMP concentrations under basal conditions (S2 Fig) [14,15].
An alternative explanation for how PDE activity creates compartmentalized cAMP responses is based on the common supposition that PDEs act as functional barriers preventing cAMP diffusion [1,5,50]. To test this hypothesis, we developed a 3D stochastic model of cAMP diffusion. When it was assumed that cAMP moved at rates equal to free diffusion (300 μm 2 /s), cAMP gradients could only be observed when the number of PDE molecules in the barrier was unrealistically high (Fig 3). Furthermore, levels of PDE activity sufficient to produce a gradient resulted in overall cAMP levels that are well below those required for activating cAMP. Decreasing the diffusion coefficient for cAMP reduced the level of PDE activity necessary to produce gradients. However, even with a diffusion coefficient of 10 μm 2 /s, it was still necessary to use artificially high levels of PDE activity, and the overall level of cAMP was still well below that necessary to activate any downstream signaling.
We also implemented a one-dimensional functional barrier model, which was possible because the diffusion length for relevant reaction (cAMP degradation by PDE) timescales is much larger than the length scale of the periodic structure of the caveolar-extracaveloar compartments. The one-dimensional functional barrier model allowed for derivation of an analytical expression for the steady-state spatial distribution of cAMP, which we used to identify the dependence of cAMP distribution on model parameters including diffusion coefficients, position of the functional barrier, concentration of PDEs and reaction kinetics. Beyond showing that there is very little compartmentation of cAMP under physiological conditions, the analytical solution indicates that the model results show very little sensitivity to variations in model parameters, an indication of robustness of the models.
The final set of simulations examined the role that physical barriers to diffusion might play in the generation of cAMP gradients. A 3D continuum model was implemented, which included subcellular structures, which acted to impede the diffusion of cAMP. It has been postulated that different "compartments" of cAMP can be carved out by creating an area of high AC surrounded by PDEs to form a barrier to prevent the cAMP produced by the AC from affecting the rest of the cell. Although the TEM images of cardiac myocytes do not show anything that would allow such a barrier to exist [51,52], this hypothesis is prominent in the literature [1,5,50], and so it is valuable to test its conceptual validity. It is possible that AKAPs might bind a greater proportion of PDEs at the edges of the cleft surrounding the AC molecules, but to date no study has suggested this type of cellular localization. If the PDE activity near the site of production exists in restricted anatomically bounded clefts, this may explain how cAMP levels near the site of production are kept low under basal conditions, preventing activation of PKA by much higher cAMP levels found throughout the rest of the cell [14,15].
The present study demonstrates that gradients consistent with those expected to result in compartmentalized responses can be produced, but only in anatomically restricted spaces, and the dimensions of those spaces are below the resolution limit of light microscopy. Therefore, one would not expect to be able to directly visualize these compartmentalized responses with techniques currently available. However, results obtained using FRET based biosensors together with the targeted application of agonists using scanning ion conductance microscopy have shown that activation of beta2-receptors produces evidence of cAMP responses localized specifically to t-tubules in adult ventricular myocytes, and that these cAMP responses do not propagate throughout the cell [30]. This is consistent with our modeling results demonstrating that cAMP production occurring in dyadic clefts along the tubules are compartmentalized.
If we assume that PDE concentration of an average cardiac myocyte is~0.1μM [53,54] and the volume is 31,400 μm 3 , this means that there are~1.3 x 10 6 PDE molecules per cell. If we further assume that there are approximately 13,000 dyadic clefts ( [55,56], 10,000-50,000) per myocyte, and all PDE activity in the cell is concentrated in these clefts, this would mean that there are 100 PDE molecules per cleft. At this concentration, we found no evidence of a cAMP gradient across the PDE barrier in the anatomical model when diffusion was set at 200 μm 2 /s. Even if we assumed that the number of PDE per cleft was 10 fold higher, this did not affect our ability to detect a gradient. However, reducing the cAMP diffusion made a significant difference. With a diffusion coefficient of 10 μm 2 /s, a cAMP gradient was observed at all PDE concentrations tested. The results of these simulations support the conclusion that PDE activity alone is not sufficient to explain compartmentation, but if diffusion of cAMP is limited by factors such as molecular crowding, PKA buffering, and anatomical barriers combined, then compartmentation may occur.
The diffusion coefficient of cAMP was determined by using fluorescein and the ϕ450 fluorophore, fluorescent molecules about the same size as cAMP that do not bind to PKA. In water, these molecules exhibit rates of free diffusion of~300 μm 2 /s, but inside cardiac myocytes the diffusion coefficient decreases to~60 μm 2 /s. This is consistent with the 4 to 5 fold decrease in mobility typically seen with molecules this size, and it has been attributed primarily to collision with other macromolecules in the intracellular environment due to molecular crowding. It turns out that cytoplasmic viscosity is only believed to be a minor player [44]. PKA can be found in both membrane and soluble cellular fractions of most cells. Our recent data suggest that PKA is targeted specifically to the mitochondrial outer membrane by A kinase anchoring proteins (AKAPs) [41]. Our future studies will be aimed at determining the quantitative effects of this anchoring on limiting the diffusion of cAMP.
It is worth noting that this work does not preclude gradients of cAMP across entire cells or non-steady state gradients. Several studies of cell motility in non-cardiac cells have shown that a cAMP gradient can exist across the cell [57]. Also, several neuronal studies have pointed to cAMP gradients as a major feature in the turning behavior of neuronal growth cones [58,59]. In these cases, the distances under consideration are significantly larger. Also, these systems have different organization of relevant enzymes; for example, one study showed by TEM imaging that significant clusters of AC molecules localize to the synapse in rat neurons [60]. However, these computational experiments show that having 10,000 steep gradients around each cleft or each caveolae is infeasible and suggest that another explanation for the observed compartmentalized nature of PKA activity must be considered.
In this study, we focused on the example of the dyadic cleft as a restricted space that would be expected to affect the generation of cAMP gradients and compartmentalized responses. Restricted spaces created by other means would be expected to have the same effect. For example, cultured neonatal cardiac myocytes may not have dyadic clefts, but they are flatter, which together with the tight packing of mitochondria beneath the plasma membrane may be another way of creating restricted spaces that contribute to compartmentation. It is also likely that factors yet to be identified contribute to compartmentalized responses in cardiac myocytes as well as other cell types.

2-dimensional continuum model
We constructed a 2 μm by 10 μm two-dimensional finite difference model representing the sarcomeric space between adjacent t-tubules of an adult ventricular myocyte (see Fig 1). Cytosolic domains associated with caveolae found in the plasma membrane of the t-tubules were modeled as 0.1 μm x 0.01 μm spaces. These caveolar domains were flanked on each side by 0.1 μm extracaveolar spaces. βARs and AC5/6 were placed in the plasma membrane associated with caveolar domains.
All the simulations were encoded in C and run on 48-Core AMD Opteron Processors. The implicit numerical method was used to integrate Eqs 1 & 2. All parameters used in the model can be found in Iancu-Harvey model [61,62]. The time step (Δt) was set to 0.001 s. Numerical results were visualized using MATLAB R2014a by The Math Works, Inc.
The concentration of cAMP in each compartment was calculated using the following equations:

Caveolar domain
G-protein activation module: cAMP produced by AC5/6: cAMP degraded by PDEs: The general formulation used for each PDE isoform (PDEx). cAMP dynamics: @cAMPðx; z; tÞ @t ¼ @cAMP AC5=6 @t À @cAMP PDE2 @t þ @cAMP PDE3 @t þ @cAMP PDE4 @t þD @ 2 cAMPðx; z; tÞ @x 2 þ D @ 2 cAMPðx; z; tÞ @z 2 Bulk domain cAMP dynamics: @cAMPðx; z; tÞ @t ¼ D @ 2 cAMPðx; z; tÞ @x 2 þ D @ 2 cAMPðx; z; tÞ @z 2 ð13Þ @cAMP @x where D is diffusion coefficient 300 μm 2 /s (Fig 2A-2C), 60 μm 2 /s ( Fig 2D) and 10 μm 2 /s (Fig 2E), and WL = 2 μm and LL = 10 μm. Definitions and initial values for model parameters were based on experimental data as described in [61,62] and shown in Tables 1 and 2.  nm. We define the z direction as orthogonal to the membrane. The boundary conditions were assumed to be no flux boundaries. The caveolar domain contained 15 β 1 ARs, which generated cAMP at 120 molecules/s. βARs and AC5/6 were placed in the plasma membrane associated with caveolar domains. cAMP freely diffused in space. MCell tracks diffusion in radial coordinates by dividing an octant of a sphere into 16384 directions (dphi % 0.005 degrees). PDE molecules were placed in the plane z = L Ã = 100 nm as functional barriers. The Interaction radius is set to default [63][64][65]. The time step (Δt) was set to 5.0 x 10 −9 s so that p b ¼ kÁs where k is binding rate, σ is surface grid density, N a is the Avogadro constant, and D is diffusion [64]. σ is set according to the following table (Table 3): The cAMP-PDE reaction was set to where K f = 1.2 x 10 7 M -1 s -1 , K b = 58.82 s -1 , and K cat = 14.70 s -1 [66]. Even with the smallest cAMP diffusion constant used in this study, the diffusion length for relevant time scales is much larger than the length scale of the caveolar/extracaveolar anatomical microstructure. Therefore, the distribution of cAMP in planes with fixed z are approximately uniform, and the distribution (effective concentration) of cAMP as a function of z can be approximated by a 1-dimensional continuum model. (See S1 Appendix for details). This model can be solved to obtain an expression for the steady state distribution of cAMP : where Unless otherwise specified, the rate constants K f , K b , and K cat are as defined above, the flux of cAMP into the domain is J B = 4.982 μm μM s -1 , and the location of the functional barrier is L Ã = 100 nm. PDE tot is the total concentration of bound and unbound PDE and is taken to be 4.1514 x 10 n μM for n = 0,1,2,3, and 4, where the concentration is averaged over the region from the plasma membrane at z = 0 to the PDE barrier at z = L Ã . (Note that the product of PDE tot and the cross-sectional area of the microdomain and L Ã is the total number of PDE molecules.)

3-dimensional continuum anatomical barrier model of the dyadic junction
For the models that included a physiological geometry, this geometry was developed from cryo-TEM images of adult mouse cardiac myocytes, as described previously [45] (see Fig 7). Briefly, a tetrahedral surface mesh was imported into Blender for finite element simulations and to smooth the sharp edges from segmentation of the TEM images (using Blamer) that would impede numerical modeling and lead to artifacts. BLAMer is a plug-in for the Open Source Blender visualization environment (http://www.blender.org) that provides an interactive interface to the GAMer (Geometry-preserving Adaptive Mesher) tool (http://nbcr.ucsd. edu/?page_id=1131) from the FEtk (Finite Element ToolKit) software package (http://nbcr. ucsd.edu/?page_id=495) maintained and distributed by the NIH-supported National Biomedical Computation Resource. GAMer produces high-quality simplex meshes of surfaces and volumes and was used via BLAMer by Hake et al. [45] to mesh the myocyte dyadic cleft anatomy from 3D electron tomographic data. This mesh was resegmented and imported into Virtual Cell. This mesh included two t-tubules surrounded by SR as well as two mitochondria.
Supporting Information