Computational Modeling Reveals Frequency Modulation of Calcium-cAMP/PKA Pathway in Dendritic Spines.

Dendritic spines are the primary excitatory postsynaptic sites that act as subcompartments of signaling. Ca2+ is often the first and most rapid signal in spines. Downstream of calcium, the cyclic adenosine monophosphate (cAMP)/protein kinase A (PKA) pathway plays a critical role in the regulation of spine formation, morphological modifications, and ultimately, learning and memory. Although the dynamics of calcium are reasonably well-studied, calcium-induced cAMP/PKA dynamics, particularly with respect to frequency modulation, are not fully explored. In this study, we present a well-mixed model for the dynamics of calcium-induced cAMP/PKA dynamics in dendritic spines. The model is constrained using experimental observations in the literature. Further, we measured the calcium oscillation frequency in dendritic spines of cultured hippocampal CA1 neurons and used these dynamics as model inputs. Our model predicts that the various steps in this pathway act as frequency modulators for calcium, and the high frequency of calcium input is filtered by adenylyl cyclase 1 and phosphodiesterases in this pathway such that cAMP/PKA only responds to lower frequencies. This prediction has important implications for noise filtering and long-timescale signal transduction in dendritic spines. A companion manuscript presents a three-dimensional spatial model for the same pathway.

Significance

Cyclic adenosine monophosphate (cAMP)/protein kinase A (PKA) activity triggered by calcium is an essential biochemical pathway for synaptic plasticity, regulating spine structure, and long-term potentiation. In this study, we predicted that for a given calcium input, adenylyl cyclase 1 and phosphodiesterase 1 kinetics reflect both the high and the low frequencies with different amplitudes and cAMP/PKA acts as a leaky integrator of calcium because of frequency attenuation by the intermediary steps. These findings have implications for cAMP/PKA signaling in dendritic spines in particular and neuronal signal transduction in general.

Introduction

Influence of Ca2+ in response to neurotransmitter-mediated ion channel opening is one of the first steps in synaptic plasticity (1,2). Transient high calcium elevation in postsynaptic domains results in the activation of multiple protein kinases, including calcium-calmodulin-dependent kinase II (CaMKII) and cyclic adenosine monophosphate (cAMP)-dependent protein kinase A (PKA) among others (3, 4, 5, 6, 7, 8). These protein kinases are associated with long-term potentiation (LTP), memory formation, consolidation, and retrieval and are known to have distinct roles in various aspects of these processes (see Table 1 in Giese and Mizuno (9)).

The dynamics of the protein kinases associated with LTP have been studied using systems biology approaches for more than two decades now (4,10, 11, 12). In 1999, a seminal work by Bhalla and Iyengar (13) highlighted the role of coupled signaling between PKA, CaMKII, and Mitogen-activated protein kinase (MAPK) and how they govern the emergent properties of signaling in neurons. Subsequently, Hayer and Bhalla (14) expanded this model to include trafficking of α-amino-3-hydroxy-5-methyl-4-isoxazoleproprionic acid (AMPA) receptors into the synapse and showed that the protein kinase activity and signaling cross talk can affect trafficking dynamics. These models are examples of comprehensive pathway maps that account for the extensive detailed interactions between various biochemical components in neurons. Separately, detailed models for specific pathway components have also been developed. For example, the dynamics of CaMKII as a function of calcium influx have been well-studied both experimentally (15, 16, 17) and theoretically (18,19). These studies have established the role of CaMKII as a critical molecule in synaptic plasticity and the role played by CaMKII bistability in determining LTP versus long-term depression (LTD) (20, 21, 22).

Tightly coupled with calcium dynamics through various signaling pathways are the dynamics of another important second messenger, cAMP (23,24). Two types of enzymes play important roles in driving parallel Ca2+ and cAMP oscillations: adenylyl cyclases (ACs) and phosphodiesterases (PDEs) (25). ACs, which synthesize cAMP from ATP, are either directly or indirectly regulated by Ca2+ signaling. Different AC isoforms and their spatial organization in cells are thought to play an important role in regulating the spatiotemporal dynamics of cAMP (26). AC1 is an isoform of AC that is highly expressed in the brain and is directly activated by Ca2+ through the calcium-calmodulin complex. PDEs hydrolyze cAMP to AMP and are regulated by cAMP-dependent PKA activity downstream of cAMP (27,28). The isoforms PDE1 and PDE4 are highly expressed in the hippocampus and are modulated by the calcium-calmodulin complex and phosphorylation by PKA, respectively (29).

The primary effector of cAMP, PKA, is a heterotetrameric protein kinase consisting of two catalytic and two regulatory subunits. Postsynaptically, calcium influx through N-methyl-D-aspartate receptors (NMDARs) triggers the synthesis of cAMP by ACs and induces PKA activity. PKA activation plays an important role in LTP and LTD (30) and neurite outgrowth (31), resulting in memory formation (32) and emotional expression, among other cognitive processes (33). Within dendritic spines, PKA is enriched after stimulation of upstream signaling pathways (34,35) and structural plasticity (36). Spinal PKA activity has also been associated with maintenance of synaptic strength (35). Several mechanistic models have been developed to model cAMP signaling pathway dynamics (37, 38, 39, 40, 41, 42, 43, 44). These models have investigated the mechanisms by which molecular components such as NMDARs and β-adrenergic receptors (βARs) are involved in LTP through the cAMP pathway (43,44). Some others have studied localized cAMP signaling by PKA-mediated phosphorylation of PDE4 (40) or localization and colocalization of PKA by A-kinase-anchoring proteins in different subcellular domains (spine versus dendrite) (41,42). However, despite the critical role of PKA in signal transduction and physiological response within dendritic spines, little is known about how calcium-cAMP/PKA cross talk occurs in dendritic spines, particularly in response to periodic oscillations of calcium.

One of the interesting dynamical features of calcium influx into postsynaptic spines is the frequency of oscillations (45, 46, 47). Intracellular calcium oscillations ranging from seconds to minutes are known to be responsible for distinct physiological roles in the nervous system (48, 49, 50, 51). Oscillatory gating dynamics of plasma-membrane-embedded channels or cyclic calcium release from intracellular stores can generate calcium spikes (23). Changes in the amplitude and frequency of these oscillating signals along with the intracellular concentrations of these second messengers can trigger diverse cellular responses in neurons (23,52, 53, 54). Studies have shown that calcium oscillates spontaneously in spines (55) and that CaMKII activity is sensitive to calcium frequency (56, 57, 58, 59, 60). Given these observations, we asked what, if any, are the relationships between Ca2+ frequency and cAMP/PKA dynamics? To answer this question, we built a mathematical model of calcium-induced cAMP/PKA dynamics that integrates key steps in the modulation of signaling components in the cAMP/PKA pathway. This model was constrained using experimental data available in the literature. Experimental measurements of calcium oscillations in dendritic spines show that calcium oscillates spontaneously in spines and these measurements are used to generate predictions of cAMP/PKA kinetics for different calcium inputs. Our model predicts that components in this pathway serve as frequency modulators of calcium dynamics such that cAMP/PKA only picks up the lower frequency of calcium influx; that is, cAMP/PKA is a leaky integrator of calcium dynamics. This finding sheds light on the kinetic mechanisms for robustness of cAMP/PKA responses in dendritic spines.

Materials and Methods

Experimental methods

Plasmids

GCaMP6f was purchased from Addgene (Cambridge, MA).

Ethical approval

All animal experiments in this study were conducted in accordance with the National Institutes of Health Guide for the care and use of laboratory animals, and all protocols were approved by the University of California San Diego Institutional Animal Care and Use Committee. All surgical procedures were terminal, and isoflurane was used as an anesthetic to prevent any pain or discomfort.

Rat primary hippocampal neuron culture

Hippocampi were dissected from E19 Sprague-Dawley male and female rat pups (Charles River, Wilmington, MA) in ice-cold HEPES-buffered saline solution (Gibco, Gaithersburg, MD), supplemented with penicillin-streptomycin (Lifetech, Guangdong, China), glucose, and additional HEPES. Hippocampi were dissociated using papain (Worthington Biochemical, Lakewood, NJ). Cells were washed and resuspended in phenol-red-free Neurobasal Medium (Lifetech) supplemented with L-glutamine (Stem Cell Technologies, Vancouver, Canada) and SM1 (NeuroCult; Stem Cell Technologies). Mixed neuron and astrocyte cell suspension was plated at a concentration of 200,000 cells/mL onto glass bottom 35 mm petri dishes (CellVis, Mountain View, CA) coated with poly-L-lysine (Sigma-Aldrich, St. Louis, MO). Cells were incubated in a humidified 37°C, CO2 incubator (HeraCell, Thermo Fisher Scientific, Waltham, MA). At DIV (days in vitro) 20, cells were transfected with GCaMP6f using Lipofectamine 2000 (Invitrogen, Carlsbad, CA) for 48 h before imaging.

Live-cell imaging

Neurons were imaged in standard medium using a resonant-scanning A1R confocal microscope (Nikon, Tokyo, Japan) equipped with a 60× 1.49 NA oil immersion objective, 489.1 nm laser, 525/50 nm emission filter, 488 dichroic, and DU4 GaASP detectors. Cells were kept at 37°C and 5% CO2 using a stage-top incubation system (Okolab, Ambridge, PA). Images of basal calcium dynamics were acquired using NIS Elements AR 5.10.01 at four frames per second for 4–5 min. In all experiments, images were acquired using identical parameters and settings (laser excitation power, acquisition time, and time-lapse interval). We replicated the findings using at least three independent culture preparations.

Identification of dendritic spines

Dendritic spines were selected from proximal and distal regions of dendrites. The images shown in Fig. 6 are portions of a dendrite image. For analysis, we selected only dendritic spines in which the head and neck, attached to the dendrite shaft, were visible in a single plane of focus, which was necessary to achieve sufficiently rapid image acquisition. Only experiments in which focus was perfectly maintained throughout the recording session were included in the time-lapse analyses. Dendritic spine selection was designed to minimize sampling bias by including approximately equal numbers of dendritic spines of different morphology, which was assessed on multiple planes of focus before starting the single-plane acquisition.

Figure 6

Model predictions for AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA for experimental calcium input measured from GCaMP6f expressed in DIV 22 primary rat hippocampal neurons. Fluorescence intensity measurements shown for calcium have an arbitrary starting time, and the model system at time zero (start of the recording) is already in equilibrium. The calcium concentration is estimated based on the relative increases and decreases in calcium to the minimum intensity measured. (A) The calcium input is at a low range, almost on the range of resting calcium concentration. This calcium amount slows down the rate of AC1 and PDE1 activation by the calcium-calmodulin complex, and as a result, only a low amount of cAMP is produced. AC1 ⋅ CaM and PDE1 ⋅ CaM seem to pick up only calcium peaks, whereas these calcium dynamics seem to barely affect the cAMP and PKA dynamics over the time course of recorded calcium (4 min). (B) The second set of measured calcium in the spine shows a higher amount of calcium, which accelerates the rate of AC1 and PDE1 activation and cAMP and PKA formation. Because of the higher range of calcium oscillations in this data set, AC1 ⋅ CaM, PDE1 ⋅ CaM, and cAMP seem to be affected more by calcium transients. (C) This set of measured calcium data shows the highest recorded calcium. Because of the high range of calcium oscillations (10 μM), the induced calcium transients are significant in this case, and even PKA shows a subtle effect of the most prominent calcium spikes. The data are representative of n = 10 neurons from three independent experiments. Representative images at the indicated time points are shown for each spine selected for analysis. Scale bars, 10 μm. To see this figure in color, go online.

Image analysis

The goal of the live-cell imaging data was to verify the existence of calcium oscillations in postsynaptic spines and to extract the frequency of calcium oscillations as input for the model to test model robustness. Initial image analysis was done using NIS Elements AR Analysis 5.10.01. To measure calcium in spines, images were deconvoluted in two dimensions using 20 iterations of blind deconvolution with medium noise. Individual spines were identified, and intensity was measured over time. Fluorescence intensities were background-corrected by subtracting intensity from a region with no transfected cells from spine intensities and represented as ΔF/F, where the minimal fluorescence intensity was subtracted from the instantaneous fluorescence intensity and divided by the minimal fluorescence intensity. Further processing was done using Microsoft Excel and MATLAB (The MathWorks, Natick, MA). Figure preparation was done using OMERO.

Mathematical model development

To simulate the temporal dynamics of calcium-stimulated cAMP/PKA pathway in dendritic spines, we developed a well-mixed model that accounts for the dynamics of the different biochemical species shown in Fig. 1. Here, we list the model assumptions, identification of the key components in the pathway, and the governing equations.

Figure 1

Schematic of the neuronal cAMP-PKA pathway. The main biochemical components in the pathway modeled are shown. Each arrow depicts the corresponding reactions between the different biochemical components of the model. Dashed arrows represent enzymatic reactions modeled using Michaelis-Menten kinetics, and solid arrows represent reactions modeled using mass action kinetics. The reaction numbers on each arrow correspond to the reactions listed in Table 1. To see this figure in color, go online.

Table 1

cAMP-PKA Pathway Reactions, Reaction Types, and Reaction Rates Used in the Model

	List of Reactions	Reaction Type	Reaction Rate
Module 1: Ca2+-CaM Complex Formation
1	2Ca2+ + CaM ⇌kb_1kf_1 Ca2 ⋅ CaM	MA	R = kf1[Ca2+]2[CaM] − kb1[Ca2CaM]
2	2Ca2+ + Ca2 ⋅ CaM ⇌kb_2kf_2 Ca4 ⋅ CaM	MA	R2 = kf2[Ca2+]2[Ca2CaM] − kb2[Ca4CaM]
Module 2: AC1 Activation by Ca2+-CaM Complex
3	AC1 + Ca2 ⋅ CaM⇌kb_3kf_3 AC1 ⋅ Ca2 ⋅ CaM	MA	R3 = kf3[AC1][Ca2CaM] − kb3[AC1Ca2CaM]
4	2Ca2+ + AC1 ⋅ Ca2 ⋅ CaM ⇌kb_4kf_4 AC1 ⋅ Ca4 ⋅ CaM	CB	R_4=(k_f4[AC1Ca_2CaM][Ca^2+]/K_m4+[Ca^2+])−k_b4[AC1Ca_4CaM]
Module 3: PDE1 Activation by Ca2+-CaM Complex
5	PDE1 + Ca2 ⋅ CaM ⇌kb_5kf_5 PDE1 ⋅ Ca2 ⋅ CaM	MA	R5 = kf5[PDE1][Ca2CaM] − kb5[PDE1Ca2CaM]
6	2Ca2+ + PDE1 ⋅ Ca2 ⋅ CaM ⇌kb_6kf_6 PDE1 ⋅ Ca4 ⋅ CaM	CB	R_6=(k_f6[PDE1Ca_2CaM][Ca^2+]/K_m6+[Ca^2+])−k_b6[PDE1Ca_4CaM]
Module 4: cAMP Production
7	ATP + AC1 ⋅ Ca4⋅ CaM→Kcat_7 cAMP	MM	R_7=(K_cat7[AC1Ca_4CaM][ATP]/K_m7+[ATP])−K_r7[cAMP]
Module 5: PKA Activation
8	R2C2 + 2cAMP ⇌kb_8kf_8 R2C2 ⋅ cAMP2	MA	R8 = kf8[R2C2][cAMP]2 − kb8[R2C2cAMP2]
9	R2C2 ⋅ cAMP2 + 2cAMP⇌kb_9kf_9 R2C2 ⋅ cAMP4	MA	R9 = kf9[R2C2cAMP][cAMP]2 − kb9[R2C2cAMP4]
10	R2C2 ⋅ cAMP4⇌kb_10kf_10 2PKAc + R2 ⋅ cAMP4	MA	R10 = kf10[R2C2cAMP4] − kb10[PKAc]2[R2cAMP4]
Module 6: PDE Phosphorylation by PKA
11	PDE1 + PKAc→Kcat_11 PDE1P	MM	R_11=(K_cat11[PKAc][PDE1]/K_m11+[PDE1])−K_r11[PDE1P]
12	PDE4 + PKAc→Kcat_12 PDE4P	MM	R_12=(K_cat12[PKAc][PDE4]/K_m12+[PDE4])−K_r12[PDE4P]
Module 7: cAMP Degradation
13	cAMP + PDE1 ⋅ Ca4⋅ CaM→Kcat_13 AMP	MM	R_13=(K_cat13[PDE1Ca_4CaM][cAMP]/K_m13+[cAMP])−K_r13[AMP]
14	cAMP + PDE4→Kcat_11 AMP	MM	R_14=(K_cat14[PDE4][cAMP]/K_m14+[cAMP])−K_r14[AMP]
15	cAMP + PDE4P→Kcat_15 AMP	MM	R_15=(K_cat15[PDE4P][cAMP]/K_m15+[cAMP])−K_r15[AMP]
Module 8: AC1 Activation and cAMP Production by Gs
16	ISO + β2AR ⇌kb_16kf_16 ISO ⋅ β2AR	MA	R16 = kf16[ISO][β2AR] − kb16[ISO ⋅ β2AR]
17	ISO ⋅ β2AR + Gsαβγ⇌kb_17kf_17 ISO ⋅ β2AR ⋅ Gsαβγ	MA	R17 = kf17[ISO ⋅ β2AR][ Gsαβγ] − kb17[ISO ⋅ β2AR ⋅ Gsαβγ]
18	ISO ⋅ β2AR ⋅ Gsαβγ→kf_18 ISO ⋅ β2AR ⋅ Gβγ + GsαGTP	MA	R18 = kf18[ISO ⋅ β2AR ⋅ Gsαβγ]
19	β2AR + Gsαβγ⇌kb_19kf_19β2AR ⋅ Gsαβγ	MA	R19 = kf19[β2AR][Gsαβγ] − kb19[β2AR ⋅ Gsαβγ]
20	ISO + β2AR ⋅ Gsαβγ⇌kb_20kf_20 ISO ⋅ β2AR ⋅ Gsαβγ	MA	R20 = kf20[ISO][β2AR ⋅ Gsαβγ] − kb20[ISO ⋅ β2AR ⋅ Gsαβγ]
21	ISO ⋅ β2AR ⋅ Gsαβγ→kf_21 ISO ⋅ β2AR ⋅ Gβγ + GsαGTP	MA	R21 = kf21[ISO ⋅ β2AR ⋅ Gsαβγ]
22	ISO ⋅ β2AR ⋅ Gβγ→kf_22 ISO ⋅ β2AR + Gβγ	MA	R22 = kf22[ISO ⋅ β2AR ⋅ Gβγ]
23	GsαGTP →kf_23 GsαGDP	MA	R23 = kf23[GsαGTP]
24	GsαGDP + Gβγ→kf_24 Gsαβγ	MA	R24 = kf24[GsαGDP][ Gβγ]
25	AC1 + GsαGTP ⇌kb_25kf_25 AC1 ⋅ GsαGTP	MA	R25 = kf25[AC1][GsαGTP] − kb25[AC1 ⋅ GsαGTP]
26	AC1 ⋅ GsαGTP + Ca4CaM ⇌kb_26kf_26 AC1 ⋅ GsαGTP ⋅ Ca4CaM	MA	R26 = kf26[AC1 ⋅ GsαGTP][Ca4CaM] − kb26[AC1 ⋅ GsαGTP ⋅ Ca4CaM]
27	ATP + AC1 ⋅ GsαGTP ⋅ Ca4CaM→Kcat_27 cAMP	MM	R_27=(K_cat27[AC1⋅G_sαGTP⋅Ca_4CaM][ATP]/K_m27+[ATP])

Enzymes are shown in bold type. CB, cooperative binding; MA, mass action; MM, Michaelis-Menten.

Model assumptions

Compartment size

We model the spine as a single compartment with a volume of 0.1 fL (10−16 L) based on the average volume of a single dendritic spine (61). All the biochemical species, including the membrane-bound ones, are well-mixed across this compartment. A separate spatial model is presented in (62).

Timescales

We focus on cAMP/PKA dynamics at the minute timescale. Calcium input was prescribed as sinusoidal functions with exponential decays based on experimental observations in the literature (63) and our experimental measurements. Different calcium inputs are described in detail below.

Reaction types

In general, the catalytic reactions such as enzymatic reactions, phosphorylation by protein kinases, and dephosphorylation by phosphatases are modeled using Michaelis-Menten kinetics. All other reactions are modeled using mass-action kinetics, assuming that the concentrations are present in large quantities. The binding reactions are modeled using either mass-action or cooperative kinetics. Catalytic reactions representing enzymatic activation or inhibition are modeled using Michaelis-Menten kinetics. For binding/unbinding reactions, we started with mass-action kinetics, and then if experimental data showed high nonlinearity or cooperativity, we used kinetic rates with higher nonlinearities (Fig. S1). For enzyme activations, especially the final step of activations, we used Michaelis-Menten kinetics. The reactions and corresponding reaction rates are given in Table 1.

Enzyme regulations by calcium-calmodulin complex

We assumed that the enzymes regulated by calcium-calmodulin, AC1 and PDE1, are fully activated only when all four binding sites of calmodulin are bound by calcium ions (64,65). These fully activated enzymes are shown as AC1 ⋅ CaM (AC1 ⋅ Ca4 ⋅ CaM) and PDE1 ⋅ CaM (PDE1 ⋅ Ca4 ⋅ CaM).

Kinetic parameters

Kinetic parameters for these reactions were estimated using the COPASI simulator (66) and are given in Table 2. The model was constrained such that dose-response curves from the simulations matched the experimentally reported dose-response curves. To estimate the model parameters, we primarily used data from brain tissues. For AC1 activation by Ca2+, we used mouse hippocampus data (67). For cAMP-ATP dynamics, data from mouse neuroblastoma (N18TG2 cells) were used (68). For PDE1 activation and PDE1/PKA phosphorylation, we used data from bovine brain (60 kDa) (69,70). In cases in which data from brain tissues were not available, we used other cell types. For example, for PDE4/PKA phosphorylation, we used data from SF9 cells (Spodoptera frugiperda Sf21 insect cells) (71). The vast majority of experimental data available in the literature are based on dose-response curves and steady-state data. We mostly found steady-state data to fit the parameters because time course data were not often available in the literature. Although kinetic parameters have not been fitted to time course data, the oscillation timescales and concentration levels have been considered by constraining the lower bound and upper bound of calculated parameters.

Table 2

Reaction Parameters Calculated for the Model

Reaction Rate	kf	Unit	kb	Unit	Kcat	Unit	Km	Unit	Kr	Unit
R1	0.1	μM−2 ⋅ s−1	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R2	0.5	μM−2 ⋅ s−1	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R3	45.46	μM−1 ⋅ s−1	1	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R4	6.89	s−1	0.39	s−1	N/A	N/A	93.02	μM	N/A	N/A
R5	87.65	μM−1 ⋅ s−1	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R6	49.80	s−1	1.09	s−1	N/A	N/A	6.95	μM	N/A	N/A
R7	N/A	N/A	N/A	N/A	7.41	s−1	193.73	μM	0.44	s−1
R8	0.19	μM−2 ⋅ s−1	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R9	27.89	μM−2 ⋅ s−1	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R10	0.21	s−1	1	μM−2 ⋅ s−1	N/A	N/A	N/A	N/A	N/A	N/A
R11	N/A	N/A	N/A	N/A	0.10	s−1	0.42	μM	0.12	s−1
R12	N/A	N/A	N/A	N/A	7.67	s−1	1.34	μM	0.02	s−1
R13	N/A	N/A	N/A	N/A	1a	s−1	35a	μM	0	s−1
R14	N/A	N/A	N/A	N/A	8.66	s−1	1.21	μM	0.10	s−1
R15	N/A	N/A	N/A	N/A	1.59	s−1	0.84	μM	0	s−1
R16	5.556	μM−1 ⋅ s−1	5	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R17	0.6	μM−1 ⋅ s−1	0.001	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R18	20	s−1	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R19	0.04	μM−1 ⋅ s−1	0.0003	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R20	2.5	μM−1 ⋅ s−1	0.5	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R21	20	s−1	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R22	80	s−1	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R23	1	s−1	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R24	100	μM−1 ⋅ s−1	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A
R25	38.5	μM−1 ⋅ s−1	10	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R26	6	μM−1 ⋅ s−1	0.9	s−1	N/A	N/A	N/A	N/A	N/A	N/A
R27	10	μM−1 ⋅ s−1	2273	s−1	56.84	s−1	233b	μM	N/A	N/A

N/A, not applicable.

All the kinetic parameters are calculated by fitting the experimental data to the kinetic equations. The kinetic parameters of reaction 13 are directly extracted from (70) and the kinetic parameters of reactions 16–27 are directly extracted from (43).

.

.

Based on these assumptions, we constructed a well-mixed model of the calcium-induced cAMP/PKA pathway in dendritic spines.

Key components and governing equations of the model

The primary focus of this study is the regulation of cAMP by predominant Ca2+-sensitive isoforms of ACs and PDEs in the brain and, more specifically, in the hippocampus. The isoforms of AC expressed at high levels in rat brain are AC1, AC2, and AC5 (72). AC1 is predominantly found in the hippocampus and cerebellum, AC5 is restricted to the basal ganglia, and AC2 is more widely expressed but at much lower levels (72). The predominantly Ca2+/calmodulin-stimulated AC1 is selectively expressed in areas associated with learning and memory, including the cortex, hippocampus, and cerebellum. The Ca2+-inhibitable AC5 is exclusively expressed in striatal medium-sized neurons (72). Therefore, only AC1 has been considered in this model. Of all the 11 known members of the PDE family, only PDE1 and PDE4 are directly regulated by a rise in calcium (73). After identifying the key components of the cAMP pathway, using the existing experimental data in the literature, we designed the the model using the following eight modules:

Ca2+/CaM complex formation

Calcium can enter the spine through NMDAR-mediated influx at the plasma membrane or through the activation of ryanodine and IP3 receptors from internal stores (74). Resting calcium is in the range of 0.05–0.1 μM and can rise to 10–100-fold when neurons are activated (75). Calmodulin is a small protein involved in the regulation of many calcium-dependent events (76, 77, 78). Calmodulin has two calcium-binding sites at each N-terminus and C-terminus lobe with different Ca2+ affinities (79). Calcium has a higher affinity for C-lobes than the N-lobes (80) but the N-lobe binds calcium faster (81). Based on these observations, we modeled the formation of the calcium-calmodulin complex using two sequential binding reactions with mass-action kinetics. For simplicity, we did not consider the cooperative binding within each of these individual lobes (82). Reaction 1 in Table 1 accounts for binding of the first two calcium ions to the C-lobe, resulting in the formation of Ca2 ⋅ CaM, and reaction 2 in Table 1 accounts for the binding of third and fourth calcium ions to the N-lobe, resulting in the formation of Ca4 ⋅ CaM. Kinetic parameters of these reactions (kf1, kb1, kf2, kb2) were obtained by fitting experimental data (80) to the model in COPASI (Fig. 2 A; Table 2).

Figure 2

Steady-state parameter scans to evaluate the model based on available experimental data in the literature: (A) saturation curve of C- and N-lobes of calmodulin predicted by our model and a comparison with experimental data from (80). (B) Ca2+-stimulated AC activity predicted by our model and a comparison with experimental data from (67). (C) cAMP production at increasing concentrations of ATP predicted by our model and a comparison with experimental data from (68). (D) Activation of RIIβ holoenzymes by cAMP in A- and B-domains predicted by our model and a comparison with experimental data from (102). (E) Activation of PDE isoenzyme PDE1A by different concentrations of Ca2+ before and after the phosphorylation of the isoenzyme by PKA, predicted by our model, and a comparison with experimental data from (69). (F) PDE4A-catalyzed cAMP hydrolysis at increasing cAMP concentration before and after PKA phosphorylation, predicted by our model, and a comparison with experimental data from (71). It must be noted that none of the experimental data sets in the literature reported error bars or standard deviations for the measured data points. To see this figure in color, go online.

AC1 activation by Ca2+/CaM complex

Calcium-calmodulin-stimulated AC1 plays an important role in learning, memory formation, and LTP by coupling calcium activity to cAMP in neurons (83). Unlike other members of the AC family, AC1 is not stimulated by G-coupled receptors unless it is activated by intracellular calcium (84). AC1 responds to NMDA-mediated calcium elevations (and not calcium released from intracellular stores) and is mostly located close to the postsynaptic densities of dendritic spines (85). AC1 is the predominant and neurospecific calcium-stimulated AC that is mostly expressed in dendrites and axons and at synapses (86). At very high concentrations of calcium (>10 μM), AC1 is inhibited (87). Our model takes this effect into account by fitting experimental data from (67) to the AC1 activity for steady-state Ca2+ (Fig. 2 B). AC1 synthesizes cAMP from ATP and has been shown to be associated with hippocampal-dependent learning abilities (85). We modeled AC1 activation by calcium as a two-step process. First, AC1 binds to Ca2 ⋅ CaM by mass-action kinetics to form AC1 ⋅ Ca2 ⋅ CaM (reaction 3 in Table 1). Second, two other Ca2+ ions bind to AC1 ⋅ Ca2 ⋅ CaM in a cooperative binding reaction to form AC1 ⋅ Ca4 ⋅ CaM (reaction 4 in Table 1). Only fully activated AC1 ⋅ Ca4 ⋅ CaM can catalyze cAMP production. The first term in R4 accounts for cooperative kinetics, whereas the second term accounts for basal degradation. As previously mentioned, the calcium/calmodulin binding in this study was simplified to two binding reactions: the first two calcium ions binding calmodulin and producing Ca2 ⋅ CaM and then Ca2 ⋅ CaM binding to two more calcium ions and producing Ca4 ⋅ CaM. Based on Masada et al.’s (64) study showing AC1 is not activated by half-occupied calmodulin, we assumed that only the product of the second reaction (Ca4 ⋅ CaM) can activate AC1. The parameters used in this module are listed in Table 2, and a comparison against experimental data is shown in Fig. 2 B.

PDE1 activation by Ca2+/CaM complex

PDE1 is the other key enzyme that regulates calcium-cAMP interactions (65). PDE1A (the PDE1 subtype used in our model) is widely expressed at high levels in the hippocampus (88). PDE1 is the only type of PDE that is activated by calcium-calmodulin (89,90). PDE1 becomes activated by an increase in intracellular calcium level and has low cAMP affinity in comparison with other members of the PDE family (91). PDE1 isozymes show low cAMP affinities and low basal PDE1 activities. Based on the study by Sharma et al. (92), the isozymes function mainly during cell activations when both cAMP and Ca2+ concentrations are elevated. The phosphorylation of PDE1A (brain 60 kDa PDE1) by PKA increases the calcium concentration required for activation of PDE1A by decreasing the enzyme’s affinity for calmodulin (92). Unlike AC1, PDE1A does not discriminate between different sources of calcium (73). We modeled PDE1 activation by calcium as a two-step process: first, PDE1 binds to Ca2 ⋅ CaM by mass-action kinetics to form PDE1 ⋅ Ca2 ⋅ CaM (reaction 5 in Table 1). Second, two more Ca2+ ions bind to PDE1 ⋅ Ca2 ⋅ CaM in a cooperative binding reaction to form PDE1 ⋅ Ca4 ⋅ CaM (reaction 6 in Table 1). Only fully activated PDE1 ⋅ Ca4 ⋅ CaM can hydrolyze cAMP into AMP. The first term in R6 accounts for cooperative kinetics, whereas the second term accounts for basal degradation. The parameters used in this module are listed in Table 2, and a comparison against experimental data from (69) is shown in Fig. 2 E.

cAMP production

ACs catalyze the conversion of ATP into cAMP. cAMP is implicated in both LTP and LTD by activating PKA (93, 94, 95, 96). In our model, fully activated AC1, AC1 ⋅ Ca4 ⋅ CaM, catalyzes the formation of cAMP from ATP (reaction 7 in Table 1). This reaction is modeled by Michaelis-Menten kinetics with a basal degradation term. The basal activity of the type I AC is regulated by calcium-calmodulin (97). In the absence of Ca2+ influx, AC1 shows a low level of basal activity, producing modest levels of cAMP (98). Therefore, the AC1 basal cAMP production has not been explicitly modeled as an individual reaction. The steady-state amount of generated cAMP with ATP is evaluated using experimental data from (68) shown in Fig. 2 C.

PKA activation

PKA has two catalytic subunits and a regulatory dimer that binds to the catalytic subunits (99,100). RIIβ is the dominant PKA isoform expressed in the brain (101). cAMP binds to its effector, PKA, and regulates the activity of regulatory subunits. Then, PKA phosphorylates different protein kinases responsible for neuron growth and differentiation, such as LKB1 and GSK-3β in the hippocampus (91). Activation of PKA by cAMP is modeled in three steps, inspired by (101), with some modifications in the order of A- and B-domain binding of cAMP that was observed in experiments by (102). First, two cAMP molecules bind to the A-domain (reaction 8 in Table 1); then the other two cAMP molecules bind to the B-domain (reaction 9 in Table 1). Finally, when all four cAMP molecules are bound to A- and B-domains, the two catalytic subunits are released (reaction 10 in Table 1). These reaction schemas are evaluated by experimental data from (102) for steady-state cAMP concentration and PKA formation shown in Fig. 2 D, and kinetic parameters are given in Table 2.

PDE phosphorylation by PKA

Both PDE1A and PDE4A can be phosphorylated by PKA. PKA phosphorylation of PDE1A alters its affinity toward both calcium and calmodulin (69). PDE4 is a cAMP-specific PDE, and its activity is increased two- to sixfold upon PKA phosphorylation (91). Different PDE4 isoforms (PDE4A, PDE4B, PDE4D) are targeted to synapses to regulate the cAMP level and modulate synaptic plasticity in learning and memory (91). PDE4A (the PDE4 subtype used in our model) is highly expressed in the CA1 subregions of the rat hippocampus and is the major PDE4 subtype involved in mediating memory (88). The phosphorylation reactions of the two different PDE isoforms are modeled by Michaelis-Menten kinetics with a basal degradation, in which PKA acts as an enzyme in both of these reactions (reactions 11 and 12 in Table 1). These reactions are evaluated at steady state to show how PKA phosphorylation affects the activity of PDE1A and PDE4 in comparison with experiments from (69) and (71) in Fig. 2, E and F, respectively.

cAMP degradation

We assumed that cAMP can be degraded to AMP by PDE1 ⋅ Ca4 ⋅ CaM, PDE4, and PDE4P (PDE4 phosphorylated by PKA). All cAMP inhibition reactions are modeled by Michaelis-Menten kinetics with a basal degradation term, in which PDEs act as enzymes (reactions 13–15 in Table 1). cAMP hydrolysis by PDE4 and PDE4P evaluated against experimental data from (71) is shown in Fig. 2 F.

AC1 activation by Gs

There is some evidence indicating that the activity of NMDARs may be potentiated by stimulating G-protein-coupled receptors. Raman at al. (103) have shown that stimulation of βARs during excitatory synaptic transmission can increase charge transfer and Ca2+ influx through NMDARs. Gereau and Conn (104) have also tested the hypothesis that coactivation of mGluRs and a G-coupled receptor (the βAR) could lead to large increases in cAMP accumulation in hippocampus and thereby increase synaptic responses in area CA1. They report that coactivation of metabotropic glutamate receptors (mGluRs) and βARs is not accompanied by an increase in excitatory postsynaptic currents or by a decrease in synaptic inhibition in area CA1, suggesting that it is not mediated by a lasting change in excitatory or inhibitory synaptic transmission. However, coactivation of these receptors leads to a persistent depolarization of CA1 pyramidal cells with a concomitant increase in input resistance.

G can have an additional effect on AC1 activation. Stimulation of β-adrenergic receptors by an agonist such as isoproterenol (ISO) can activate the G subtype of GTP-binding protein, which stimulates AC isoforms (105). A study by Wayman et al. (84) suggests that AC1 enzyme is stimulated by G-coupled receptors in vivo when it is activated by intracellular calcium. Based on these experimental observations, reactions 15–27 (Table 1) has been reproduced from (43) exactly as is to take into account the effect of AC1 activation and cAMP production by G.

The list of the above reactions and estimated parameters are shown in Tables 1 and 2, respectively.

Ordinary differential equations

The temporal dynamics of each species, c, is defined aswhere c, i ∈ 1, 2, ..., 21, represents the concentration of the ith species as a function of time and J is the net reaction flux for the ith species. These ordinary differential equations (ODEs) and the corresponding initial concentrations of the different species are shown in Table 3. The sensitivity analysis with respect to kinetic parameters and initial conditions of the model with respect to transient cAMP concentration are shown in Fig. 3.

Table 3

Initial Concentrations and ODEs

Species	Concentration (nM)	ODE
Ca2+	100	(d[Ca^2+]/dt)=−2R_1−2R_2−2R_4−2R_6
CaM	1240	d[CaM]/dt=−R_1
Ca2 ⋅ CaM	0	(d[Ca_2CaM]/dt)=R_1−R_2−R_3−R_5
Ca4 ⋅ CaM	0	(d[Ca_4CaM]/dt)=R_2−R_26
AC1	5000	(d[AC1]/dt)=−R_3−R_25
AC1 ⋅ Ca2 ⋅ CaM	0	(d[AC1Ca_2CaM]/dt)=R_3−R_4
AC1 ⋅ Ca4 ⋅ CaM	0	(d[AC1Ca_4CaM]/dt)=R_4
PDE1	3800	(d[PDE1]/dt)=−R_5−R_11
PDE1 ⋅ Ca2 ⋅ CaM	0	(d[PDE1Ca_2CaM]/dt)=R_5−R_6
PDE1 ⋅ Ca4 ⋅ CaM	0	(d[PDE1Ca_4CaM]/dt)=R_6
ATP	2×10^6	(d[ATP]/dt)=−R_7−R_27
cAMP	0	(d[cAMP]/dt)=R_7−2R_8−2R_9−R_13−R_14−R_15+R_27
R2C2	220	(d[R_2C_2]/dt)=−R_8
R2C2 ⋅ cAMP2	0	(d[R_2C_2cAMP_2]/dt)=R_8−R_9
R2C2 ⋅ cAMP4	0	(d[R_2C_2cAMP_4]/dt)=R_9−R_10
R2 ⋅ cAMP4	0	(d[R_2cAMP_4]/dt)=R_10
PKAc	10	(d[PKAc]/dt)=2R_10
PDE1P	0	(d[PDE1P]/dt)=R_11
PDE4	1000	(d[PDE4]/dt)=−R_12
PDE4P	0	(d[PDE4P]/dt)=R_12
AMP	0	(d[AMP]/dt)=R_13+R_14+R_15
ISO	4a	(d[ISO]/dt)=−R_16−R_20
β2AR	19.2a	(d[β_2AR]/dt)=−R_16−R_19
ISO ⋅ β2AR	1.2a	(d[ISO⋅β_2AR]/dt)=R_16−R_17+R_22
Gsαβγ	3654a	(d[G_sαβγ]/dt)=−R_17−R_19+R_24
ISO ⋅ β2AR ⋅ Gsαβγ	1.2a	(d[ISO⋅β_2AR⋅G_sαβγ]/dt)=R_17−R_18+R_20−R_21
ISO ⋅ β2AR ⋅ Gβγ	0	(d[ISO⋅β_2AR⋅G_βγ]/dt)=R_18+R_21−R_22
β2AR ⋅ Gsαβγ	709.2a	(d[β_2AR⋅G_βγ]/dt)=R_19−R_20
GsαGTP	1.2a	(d[G_sαGTP]/dt)=R_18+R_21−R_23−R_25
GsαGDP	0	(d[G_sαGDP]/dt)=R_23−R_24
Gβγ	118.8a	(d[G_βγ]/dt)=R_22−R_24
AC1 ⋅ GsαGTP	52.8a	(d[AC1⋅G_sαGTP]/dt)=R_25−R_26
AC1 ⋅ GsαGTP ⋅ Ca4 ⋅ CaM	1.2a	(d[AC1⋅G_sαGTP⋅Ca_4CaM]/dt)=R_26

All the initial concentrations are estimated in COPASI using the experimental data fitted to the model parameters. Membrane-bound molecule concentrations are converted from to nM by volume to surface ratio of the compartment.

The concentrations are taken from (43).

Figure 3

Sensitivity analysis with respect to cAMP transient concentration for (A) parameters of reactions 1–12, (B) parameters of reactions 13–27, and (C) initial concentrations and compartment volume. Kinetic parameters of reactions 3–6 and initial concentrations of calmodulin, AC1, and PDE1 show the highest impact on cAMP transient concentration. To see this figure in color, go online.

Ca2+ input signals as the stimulus

The four different types of input calcium are modeled as

sinusoidal oscillations ( = 0.5 Hz) with exponential decay that is repeated every 5 min ( = 0.003 Hz):

calcium bursts with exponential decay that are being repeated every 5 min ( = 0.003 Hz):

nonoscillating calcium

in which [Ca2+] = 0.1 μM is the resting concentration of calcium in the cytoplasm;

and experimental measurements of calcium: spontaneous calcium oscillations measured in neurons are shown in Fig. 6. The measured calcium is in the range of 0.1–10 μM, with oscillations in the second timescale.

Numerical methods

The system of ODEs for this model is solved for a time course of 1800 s and interval size of 0.1 (s) through a deterministic simulation in COPASI (106). Using a deterministic algorithm (LSODA), the system of ODEs is solved with a dense or banded Jacobian when the problem is stiff and automatically selects between nonstiff (Adams) and stiff (BDF) methods by initially using the nonstiff method and dynamically monitoring data to decide which method to use (107). Parameter estimations were conducted in COPASI using the evolutionary programming method (108, 109, 110). The link to the model will be available on https://github.com/donya26/cAMP-PKA-temporal.git.

Sensitivity analysis

In the sensitivity analysis of kinetic parameters with respect to transient cAMP concentration (Fig. 3), kinetic parameters of reactions related to activation of AC1 and PDE1 enzymes by Ca2+/CaM complex (reactions 3–6) show the most significant effect on cAMP concentration (Fig. 3 A). However, kinetic parameters related to AC1 activation by G shows almost no effect on cAMP concentration (Fig. 3 B). Finally, initial concentration of these Ca2+-sensitive enzymes (AC1 and PDE1) and calmodulin shows the most significant effect on cAMP concentration (Fig. 3 C). It is important to note that the compartment size has a trivial effect on the model (Fig. 3 C). These results show the significance of regulation of cAMP/PKA pathway through Ca2+-sensitive enzyme activations. Sensitivity analyses were conducted in COPASI by numerical differentiation using finite differences (111).

Stochastic simulations

Spines are small compartments, and neuronal activity in the spines is mediated through changes in the probability of stochastic transitions between the open and closed states of ion channels (112). In this model, the focus has been on cAMP regulation, and for the sake of simplicity, we did not consider such stochastic calcium-related events. However, to make sure that the stochasticity will not change the obtained results, we also simulated the model stochastically (τ-leap method-using COPASI). Although to conduct the simulations stochastically, we had to make some changes to the model (e.g., we had to make all the reactions irreversible and change some of the reaction rates), we saw almost the same trend with deterministic and stochastic methods (Fig. S2). However, to better capture oscillation regulation of different components of the pathway, we modeled the system deterministically.

Results

Large-frequency modulation of immediate downstream effectors of calcium

The first module of the model focused on the formation of the calcium-calmodulin complex. Calcium input was modeled as a sinusoidal function with a 2 s oscillation period with an exponential decay every 5 min (Fig. 4 A; (63)). From here onwards, we refer to the larger frequency as for the oscillations in the second timescale and the smaller frequency as for the oscillations in the minute timescale. The dynamics of calmodulin consumption leading to the formation of Ca2 ⋅ CaM respond only to the low-frequency oscillations of calcium in the minute scale (Fig. 4 A). Ca2 ⋅ CaM demonstrates both low-frequency (minute-scale) and high-frequency (second-scale) oscillations (Fig. 4). In the next step, two more Ca2+ ions bind to Ca2 ⋅ CaM to form Ca4 ⋅ CaM. Ca4 ⋅ CaM also responds to second-scale oscillations of calcium; however, the oscillation amplitudes are almost 10-fold lower in comparison with Ca2 ⋅ CaM (Fig. 4 A). Another key difference between the observed oscillation patterns of Ca2 ⋅ CaM and Ca4 ⋅ CaM is that the Ca2 ⋅ CaM oscillation pattern seems to be mainly driven by CaM on the minute scale, whereas the Ca4 ⋅ CaM oscillation pattern looks almost exactly like calcium patterns with less significant oscillations on the second scale.

Figure 4

Ca2+-CaM complex formation and Ca2+-CaM-activated enzymes’ time course predicted by our model: (A) calcium input with sinusoidal oscillations (0.5 Hz) and exponential decay (0.003 Hz) as the stimulus in the model. Free calmodulin concentration over a 15 min and a 10 s time course showing only minute-scale oscillations. The first step of calcium-calmodulin binding of two calcium ions to the C-lobe of calmodulin over a 15 min and a 10 s time course. Ca2+ ⋅ CaM shows both minute-scale (0.003 Hz) and second-scale (0.5 Hz) oscillations. The second step of calcium-calmodulin binding of another two calcium ions to N-lobe of calmodulin over a 15 min and a 10 s time course. Ca4 ⋅ CaM also shows both minute-scale (0.003 Hz) and second-scale (0.5 Hz) oscillations. However, G as the other component responsible for AC1 activation shows no oscillations because it is not activated by Ca2+. (B) AC1 partial activation by Ca2 ⋅ CaM to form AC1 ⋅ Ca2 ⋅ CaM, followed by binding of two extra calcium ions to fully activate AC1 and form AC1 ⋅ Ca4 ⋅ CaM. Unlike AC1 activated by Ca2+, AC1 activated by G shows no oscillations, (C) PDE1 partial activation by Ca2 ⋅ CaM initially followed by further binding of two calcium ions to fully activate PDE1 and form PDE1 ⋅ Ca4 ⋅ CaM. To see this figure in color, go online.

AC1 and PDE1 are activated by the Ca2+/CaM complex. Ca2 ⋅ CaM binds to AC1 and forms AC1 ⋅ Ca2 ⋅ CaM. Both AC1 and AC1 ⋅ Ca2 ⋅ CaM show second- and minute-scale calcium oscillations (Fig. 4 B), and their oscillation patterns seem to be mostly driven by the Ca2 ⋅ CaM oscillation pattern shown in Fig. 4 C, with less significant second-scale oscillations. In the second step, AC1 ⋅ Ca2 ⋅ CaM fully activates AC1 by forming AC1 ⋅ Ca4 ⋅ CaM with the other two Ca2+ ions. In comparison to AC1 ⋅ Ca2 ⋅ CaM, AC1 ⋅ Ca4 ⋅ CaM shows less significant second-scale oscillations (almost 10-fold lower), and the minute-scale oscillation seems to be driven by the calcium oscillation pattern (Fig. 4 B). AC1 ⋅ Ca4 ⋅ CaM dynamics mostly follow Ca4 ⋅ CaM dynamics with less sensitivity to Ca2+. The dynamics of PDE1 activation by Ca2 ⋅ CaM is very similar to AC1. However, both PDE1 ⋅ Ca2 ⋅ CaM and PDE1 ⋅ Ca4 ⋅ CaM show much higher sensitivity to second-scale oscillations of Ca2+ (Fig. 4 C).

Small-frequency modulation of cAMP/PKA and PDE phosphorylation

After enzyme activation by the calcium-calmodulin complex, cAMP is synthesized by AC1 ⋅ Ca4 ⋅ CaM and AC1 ⋅ GGTP ⋅ Ca4 ⋅ CaM from ATP. The contribution of Ca2+-CaM-activated AC1 in cAMP production is much more substantial in comparison to G-activated AC1 and cAMP produced by AC1 ⋅ GGTP ⋅ Ca4 ⋅ CaM is almost negligible (Fig. S3). This is in agreement with experimental observations of Baker et al. (113) and Halls and Cooper (26) showing that the Gs-coupled receptor does not stimulate the Gs-insensitive AC1, suggesting calcium-calmodulin binding as the main activation mechanism for AC1s. Then PKA is activated by cAMP and activated PKA phosphorylates PDE1 and PDE4. Finally, cAMP is degraded by PDE1, PDE4, and PDE4P and shows minute-scale oscillations and leaky oscillations on the second scale (Fig. 5 A). PKA is activated by cAMP and as a result, it only shows minute-scale oscillations (Fig. 5 B). Although PDE1-P (phosphorylated PDE1) shows minute-scale oscillations (Fig. 5 C), PDE4P (phosphorylated PDE4) barely shows any oscillations (Fig. 5 D). Degradation of cAMP by PDE4 and PDE4P, which are indirectly affected by calcium oscillations, might filter the second-scale oscillations induced by AC1 ⋅ CaM in the cAMP synthesis and change the second-scale oscillations to leaky oscillations (Fig. 5 A). Thus, our model predicts that PDEs may filter the large-frequency () oscillations of the upstream components driven by calcium dynamics and capture only the smaller frequency ().

Figure 5

cAMP and PKA time course and PDE1 and PDE4 phosphorylation by PKA predicted by our model: (A) cAMP dynamics predicted by our model show mainly minute-scale oscillations and leaky oscillations on the second scale shown in the figure inset. (B) PKA dynamics predicted by our model shows only minute-scale oscillations and no second-scale oscillations (figure inset). (C) Phosphorylation of PDE1 by PKA and formation of phosphorylated PDE1 (PDE1P) that oscillates only on the minute scale. (D) PDE4 and phosphorylated PDE4 (PDE4P) formed by PKA that is indirectly affected by calcium show almost neither of the second- or minute-scale oscillations. To see this figure in color, go online.

cAMP/PKA pathway filtering of Ca2+ signals

The response of cAMP/PKA to different calcium patterns raises the following questions: do spines experience vastly varying calcium patterns? Can experimentally measured calcium patterns be used to generate predictions from our model? To answer these questions, we measured basal calcium fluctuations in dendritic spines of hippocampal neurons. Cultured neurons have been previously reported to form neuronal networks in vitro that exhibit spontaneous action potentials (114), which can be measured using reporters of neuronal activity such as calcium sensors. We expressed the genetically encoded calcium biosensor GCaMP6f (115) in primary rat hippocampal neurons. To quantify Ca2+ oscillations, we assigned calcium concentrations to correspond to changes in the measured fluorescence intensity (ΔF/F), based on previously reported intracellular calcium concentrations (75). Using the measured calcium in three spines from three neurons as the stimulus in the model, we found our model predicts different types of cAMP/PKA dynamics based on the calcium input (Fig. 6). We chose these three spines to reflect the different calcium patterns tested in Fig. 7. The first spine measured has minimal changes in fluorescence intensity and displays a low range of predicted calcium concentrations (0.1–0.2 μM) (Fig. 6 A). There seem to be no apparent calcium spikes, and the calcium concentration is barely enough to fully activate AC1 and PDE1 and show steady-state cAMP and PKA. AC1 ⋅ CaM and PDE1 ⋅ CaM integrate some calcium peaks; however, cAMP and PKA do not exhibit significant peaks. The second spine measured shows a higher fluorescence intensity that is almost 10-fold higher than the first data set (Fig. 6 B). In the second spine, AC1, PDE1, cAMP, and PKA reach a steady state. AC1 and PDE1 seem to reflect most of the calcium major peaks, whereas cAMP appears to filter out some of the less significant peaks, and PKA almost shows none of these peaks. The third spine measured shows many calcium spikes and a very high range of fluorescence intensity, which we correlated to a large range of calcium concentrations (Fig. 6 C). The slopes at the first 60 s of simulation show that the activation of AC1 and PDE1, and formation of cAMP and PKA occurs much faster in the third spine (Fig. 6 C) in comparison to the other two (Fig. 6, A and B). With high ranges of calcium transients, AC1, PDE1, and cAMP seem to integrate most of the calcium peaks, and even PKA has nominal peaks (Fig. 6 C). Thus, we find that for experimentally derived calcium transients, cAMP/PKA serve to filter out the larger frequency ().

Figure 7

Effect of calcium stimulation patterns on AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA dynamics: (A) nonoscillating calcium with an initial concentration of 1 μM induces only one peak in each of AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA. (B) Calcium bursts with exponential decays with 0.003 Hz frequency () induce minute-scale oscillations in each of AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA. (C) The combination of sinusoidal calcium oscillations of 0.5 Hz with exponential decays of 0.003 Hz induces both second-scale and minute-scale oscillations in both AC1 ⋅ CaM and PDE1 ⋅ CaM, as the plot insets show. Although cAMP barely shows a leaky oscillation on the second scale, PKA merely shows minute-scale oscillations. To see this figure in color, go online.

Ca2+ signal frequency modulates cAMP/PKA dynamics

If it is true that cAMP/PKA filters out dynamics to pick dynamics, then this must hold true for different values of these two frequencies. To further test our predictions, we designed three different calcium inputs (Fig. 7): first, no oscillations of calcium; second, is zero and only is present; and third, both and are present. In the first set of simulations, Ca2+ input has an initial constant concentration of 1 μM with no oscillations or calcium bursts during the 15 min of oscillations (Fig. 7 A). From AC1 ⋅ CaM and PDE1 ⋅ CaM to cAMP and PKA, it is evident that there is only one peak of the concentrations for all these species, with a 40–60 s delay between AC1 ⋅ CaM and PDE1 ⋅ CaM and cAMP/PKA. This behavior has also been observed by (36), in which PKA returns to basal level ∼10 min after stimulation by a pulse of calcium. For the second set of simulations, starting from the same initial calcium concentration (1 μM), we simulated a calcium burst every 5 min (Fig. 7 B). These calcium bursts at every 5 min are reflected in AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA dynamics. The range of concentrations for all of these molecules is the same as the concentration seen in the first simulation (Fig. 7 A). The third set of simulations has sinusoidal oscillations with 2 s periods in addition to the calcium bursts at every 5 min (Fig. 7 C). AC1 ⋅ CaM and PDE1 ⋅ CaM show both second- and minute-scale oscillations, whereas cAMP and PKA only show the minute-scale oscillations (Fig. 7 C). In comparison to the second simulation (Fig. 7 B), the concentration of AC1 ⋅ CaM, PDE1 ⋅ CaM, cAMP, and PKA are lower (Fig. 7 C). We also chose different values and confirmed that cAMP only picks up minute timescales (Fig. S4).

Discussion

cAMP/PKA activity triggered by Ca2+ is an essential biochemical pathway for synaptic plasticity, regulating spine structure, and LTP (36,116). In this study, we have developed a mathematical model, constrained by existing experimental data (Fig. 2), to simulate the dynamics of the cAMP/PKA pathway in response to calcium input. Analysis of this model allowed us to make the following predictions: first, for a given calcium input, AC1, and PDE1 kinetics reflect both the large and the small frequencies, albeit with different amplitudes; second, cAMP/PKA dynamics reflect only the small frequency of the calcium input; and finally, cAMP/PKA act as a leaky integrator of calcium because of frequency attenuation by the intermediary steps (Fig. 8). We also found that spines generated different spontaneous calcium frequencies in cultured neurons (Fig. 6), thus motivating the need to study multiple timescales. These findings have implications for cAMP/PKA signaling in dendritic spines in particular and neuronal signal transduction in general.

Figure 8

The frequency control in the modeled cAMP/PKA pathway from calcium stimulation through PKA. Calcium oscillates with a combination of frequencies ( = 0.5 Hz and  = 0.003 Hz). From the calcium-calmodulin complex activation of AC1 and PDE1 enzymes through cAMP production and degradation, the effect of becomes trivial. Ultimately, through the downstream pathway, remains the only significant oscillatory frequency in this pathway. To see this figure in color, go online.

The ability to modulate the output frequency in response to a signal input is an important feature in biochemical signal transduction; it ensures robustness and noise filtering. As a result, the output signal, cAMP/PKA in this case, will respond to calcium input for strong signals but not to weak signals (Fig. 7). In general, the combination of activation and inactivation kinetics along with feedback loops, many of which oscillate, may be responsible for noise attenuation in signal regulation (117,118). Noise attenuation has a significant role in regulating biological switches. These switches must sense changes in signal concentration while buffering against signal noise (119). Our model predicts that depending on the sensitivity of regulating enzymes to Ca2+, the signaling input can be modified to the desired signaling output with a different range of frequency. For example, starting with a stimulus Ca2+ signal with a combination of and oscillation frequencies, the cAMP/PKA pathway can modify the initial signal and filter it so that components downstream in the pathway such as PKA receive a signal with an oscillation frequency of (Fig. 8).

As a result, cAMP/PKA can be thought of as a leaky integrator of calcium signals (Fig. 5, A and B). This signal modification is important in the context of how cAMP/PKA directly impacts spine morphology (120), receptor trafficking (121), and transcription factor activation (e.g., cAMP-response element binding protein (CREB) and Exchange protein activated by cAMP (Epac)) (32,122). Although these transcription factors are not directly considered in our study, it is possible that the frequency modulation exhibited by the pathway we modeled plays a critical role in contributing to neuronal robustness to small fluctuations while retaining the ability to rapidly respond to large calcium influxes. To identify whether calcium or G signals dominate, we added the reactions pertaining to the Gs-mediated cAMP production (43) in our model. Interestingly, we found that the G pathway had very little effect on the conclusions of our model. The main reason for this is that calcium mediates both AC1 and PDE1 activation, thereby creating a strong incoherent feedforward loop for cAMP (123), whereas G linearly activates AC1, and PDE1 depends only on PKA feedback. Of course, the addition of other feedback loops and regulatory pathways in future efforts will dissect this cross talk better.

The calcium input to our model was inspired by measurements of spontaneous calcium dynamics in dendritic spines (Fig. 6). Ideally, we would like to use imaging to simultaneously measure spontaneous calcium and cAMP/PKA activity in dendritic spines to test our model predictions. However, despite the recent developments in fluorescence biosensing methods, imaging and coimaging of multiple-second messenger, protein kinase, and phosphatase activities in micron-sized compartments remains challenging (124). This is also true for simultaneously imaging the dynamics of both cAMP and calcium in single dendritic spines (36,125). Recently, Tang et al. measured PKA activity in dendritic spines in response to glutamate uncaging using new fluorescence resonance energy transfer probes that they developed (36). The timescale of PKA dynamics from our model (Fig. 7 A) is in good agreement with the timescale of ∼10 min (Fig. 6 of Tang et al. (36)). This indicates our developed model closely mimics reported PKA activity in subcompartments, in this case, dendritic spines.

Despite the predictive capability of our model, there are certain limitations that must be acknowledged. First, we only consider a deterministic model and have currently ignored stochastic dynamics in the system. It is possible that the frequency modulation effects presented here will be somewhat altered by including stochasticity, but we anticipate that the overall conclusions will not be affected on average if the same network is considered. An additional limitation of our work emerges in the form of challenges associated with coimaging 2 s messengers simultaneously in the same spine. We are working on developing such probes to overcome these challenges. Finally, there are indeed many protein kinases and phosphatases that affect this pathway (126, 127, 128, 129), and ongoing efforts in our group are focused on building up increasingly comprehensive network maps of these signaling reactions. Furthermore, spatial regulations of cAMP signaling are other important factors that need to be considered (130, 131, 132). Enzyme regulation seems to play a significant role in cAMP/PKA pathway dynamics and signal filtration (133, 134, 135). These enzymes can also be localized by A-kinase anchoring proteins and change cAMP/PKA dynamics further (136, 137, 138). The effect of spine size and enzyme localizations are considered in a companion study (62).

Author Contributions

D.O. and P.R. conceived the study. D.O. developed the model and conducted the simulations. D.L.S. conducted neuronal culture experiments. B.C. aided in experimental design. All authors analyzed the data and wrote the manuscript.