Age-related deterioration of pacemaker function has been documented in mammals, including humans. In aged isolated sinoatrial node tissues and cells, reduction in the spontaneous action potential (AP) firing rate was associated with deterioration of intracellular and membrane mechanisms; however, their relative contribution to age-associated deficient pacemaker function is not known. Interestingly, pharmacological interventions that increase posttranslation modification signaling activities can restore the basal and maximal AP firing rate, but the identities of the protein targets responsible for AP firing rate restoration are not known. Here, we developed a numerical model that simulates the function of a single mouse pacemaker cell. In addition to describing membrane and intracellular mechanisms, the model includes descriptions of autonomic receptor activation pathways and posttranslation modification signaling cascades. The numerical model shows that age-related deterioration of pacemaker function is related to impaired intracellular and membrane mechanisms: HCN4, T-type channels, and phospholamban functions, as well as the node connecting these mechanisms, i.e., intracellular Ca2+ and posttranslation modification signaling. To explain the restored maximal beating rate in response to maximal phosphodiesterase (PDE) inhibition, autonomic receptor stimulation, or infused cyclic adenosine monophosphate (cAMP), the model predicts that phospholamban phosphorylation by protein kinase A (PKA) and HCN4 sensitivity to cAMP are altered in advanced age. Moreover, alteration in PKA and cAMP sensitivity can also explain age-reduced sensitivity to PDE inhibition and autonomic receptor stimulation. Finally, the numerical model suggests two pharmacological approaches and one gene manipulation method to restore the basal beating rate of aged pacemaker cells to that of normal adult cells. In conclusion, our numerical model shows that impaired membrane and intracellular mechanisms and the nodes that couple them can lead to deteriorated pacemaker function. By increasing posttranslation modification signaling, the deteriorated basal and maximal age-associated beating rate can be restored to adult levels.

## Introduction

The U.S. Administration on Aging estimates that by the year 2050, 20% of the population will be 65 y and older. By the year 2025, the number of individuals age 85 and older is projected to double and reach 8 million in the United States alone. At the same time, chronic cardiovascular diseases, e.g., coronary artery atherosclerosis, hypertension, and chronic heart failure, have reached epidemic proportions, and their incidence is growing exponentially, along with risks of disability with increasing age (Ogawa et al., 1992; Ferrari, 2002). Importantly, sinoatrial node (SAN) dysfunction increases exponentially with age, reaching an incidence of 1 of every 600 people ages 65 and older (Kusumoto and Goldschlager, 1996).

The SAN rate is controlled by both the autonomic nervous system (ANS) and intrinsic pacemaker mechanisms forming the coupled-clock system (Yaniv et al., 2015b). The ANS controls the beating rate by balancing the ratio of sympathetic to parasympathetic stimulation of G protein–coupled receptors: the sympathetic system stimulates β-adrenergic receptors (β-ARs), activating adenylyl cyclase (AC); the parasympathetic system stimulates cholinergic receptors, inactivating AC. In the cell, ATP is converted to cAMP by AC. The cAMP level controls PKA activity, and both cAMP and PKA determine the activation level of intrinsic pacemaker mechanisms: increased cAMP/PKA activity leads to an increase in the beating rate and vice versa (Yaniv et al., 2015a; Behar et al., 2016). In pacemaker cells, other types of AC exist (AC1 and AC8); these ACs are activated by intracellular Ca2+, which is balanced by the interaction between membrane channels, exchangers, and pumps (membrane clock or M clock) and by internal Ca2+ storage (Ca2+ clock; Yaniv et al., 2013b). Both cAMP/PKA and intracellular Ca2+ serve as the node connecting the two clocks. Thus, deterioration in either clock’s molecules or in ANS receptors can affect pacemaker function.

Measurements of intrinsic (measured during autonomic blockade, precluding ANS effects) and maximal beating/heart rate (HR) provide evidence for deterioration in pacemaker function associated with aging. Although the basal HR (i.e., at rest) does not change in advanced age, the intrinsic HR declines with age in humans (Jose and Collison, 1970) and other mammals (Yaniv et al., 2016; Fig. 1 A). Because the intrinsic heart rate precludes ANS effects, its reduction reflects impaired activity of intrinsic pacemaker mechanisms: cAMP-PKA/Ca2+ signaling (Liu et al., 2014) and ionic channel currents (e.g., hyperpolarization-activated cyclic nucleotide–gated [HCN] channel and T-type and L-type Ca2+ channels; Larson et al., 2013; Fig. 1, B and C). However, the relative contribution of each mechanism to age-associated deterioration of pacemaker function is debatable. In contrast to the in vivo basal HR, the maximal HR declines in advanced age (Heath et al., 1981; Hagberg et al., 1985). Indeed recent data on pacemaker tissue (Liu et al., 2014; Yaniv et al., 2016) and data at the heart level (Brodde and Leineweber, 2004) have shown that the sensitivity to β-AR stimulation or phosphodiesterase (PDE) inhibition (using 3-isobutyl-1-methylxanthine [IBMX]) was reduced in advanced age. However, data have also shown that, in response to maximal autonomic receptor stimulation or maximal increase in cAMP/PKA by PDE inhibition, the age-associated reduction in maximal beating rate can be compensated for. Nonetheless, the underlying biophysical mechanisms that allow this are not known. Numerical modeling is very useful, and sometimes it is the only feasible tool to explore such biophysical mechanisms. However, to date no aged pacemaker cell model exists.

Here, we used a numerical model of mouse pacemaker cells (Kharche et al., 2011) adapted to describe the effects of autonomic receptor stimulation. The model was used to answer three specific questions: (a) What is the relative contribution of each clock to the age-associated deterioration in pacemaker function? (b) What age-associated biophysical modulation makes it possible to compensate for the reduction in maximal beating rate obtained in response to maximal PDE inhibition or β-AR stimulation? (c) Can the age-associated deterioration in pacemaker function be reversed by drug intervention or gene manipulation?

## Materials and methods

### Existing experimental results

Several groups have measured the reduction in the basal beating rate of either SAN tissue or single pacemaker cells from aged mice (Larson et al., 2013; Liu et al., 2014; Yaniv et al., 2016; Sharpe et al., 2017). Different mechanisms can cause the reduced pacemaker function: (a) reduced sarcoplasmic reticulum Ca2+ ATPase (SERCA) activity (Liu et al., 2014); (b) reduced RyR function (Liu et al., 2014); (c) shift of the half-activation potential of the HCN current (Sharpe et al., 2017); or (d) reduced L-type and T-type currents (Larson et al., 2013). Although the action potential (AP) sensitivity to β-AR (studied using isoproterenol [ISO]) or PDE inhibition (studied using IBMX) was reduced in advanced age (Yaniv et al., 2016), the maximal AP firing rate in response to these drugs (Yaniv et al., 2016) or to cAMP (Sharpe et al., 2017) was not significantly different between aged and young mice.

### Basal numerical model

We adapted the Kharche et al. (2011) model for mouse pacemaker cells. This model was primarily based on experimental data from isolated mouse SAN cells. We updated the model based on Ca2+ measurements from Liu et al. (2014). For some ionic current parameters and biochemical signaling for which no mouse experimental data were available, we used the formulation from rabbit SAN cell models. The Kharche et al. (2011) model was updated to include the equations of cAMP-PKA signaling cascades and autonomic receptor stimulation (Yaniv et al., 2015a; Behar and Yaniv, 2016; Fig. 2 A). In addition, the following changes were made to the original model: the extracellular concentration of Ca2+ was changed to 2 mM according to Maltsev and Lakatta (2009); the intracellular concentration of Na+ was changed from 8.1 to 10 mM as in the Maltsev–Lakatta model; the release rate parameter of the Ca2+ from the SR, ks, was updated to 250,000 ms−1, similar to the Maltsev–Lakatta model; the formulation of jup to describe the Ca2+ pump rate was changed according to Kurata et al. (2002); and the conductance of ICaL was increased by 18% and 38% for gCaL1.2 and for gCaL1.3, respectively. Finally, K+ and Na+ intracellular concentrations were kept constant, and Ist was removed from the model because it has no molecular basis (Maltsev et al., 2014). All model equations can be found in the supplement.

### Model strategy

The effects of aging were simulated by varying the following model parameters: SERCA activity was decreased by 30%, as observed experimentally in Liu et al. (2014) (experimental data are presented in Fig. 1 C and are taken into account in Eqs. 20 and 22), and the conductance of ICaT was decreased by 20% to simulate current density reduction during diastolic depolarization (DD), as found in Larson et al. (2013) (experimental data are presented in Fig. 1 B and are taken into account in Eqs. 23 and 24). In addition, Liu et al. (2014) showed changes in RyR function with a decrease in RyR activity in aged cells and a decrease in Ca2+ release of the local Ca2+ release (LCR) ensemble. To simulate these results, the release rate parameter of the Ca2+ from the SR, ks = 250,000 ms−1, was decreased to ks = 87,500 ms−1, i.e., decreased by a factor of 65%, corresponding to the experimental decrease in Ca2+ release of the LCR ensemble found in Liu et al. (2014) (Eqs. 29 and 30 and Fig. 1 C). The half-activation potential V0.5 was changed in the HCN current equations from 104.2 mV in adult cells to 110 mV in aged cells (Sharpe et al., 2017; taken into account in Eqs. 14 and 17). The resulting decrease in the AP firing rate caused by aging was compared with published experimental data under basal conditions. In addition, the expression for the I–V curve of the HCN current was modified in the aged model (Fig. 2 B) to match the maximal shift found experimentally in the I–V curve under cAMP saturation (Sharpe et al., 2017), which was shown to be different from that of adults; ∼17 mV in aged cells versus ∼7 mV in adult cells (Sharpe et al., 2017). These values defined the maximal activation of V0.5 by cAMP in adult and aged cells in our model (Fig. 2 B).

Cage release of cAMP was modeled as in our previous work (Behar and Yaniv, 2016). The effect of ISO was simulated using the equations developed in Yaniv et al. (2015a) and updated in Behar and Yaniv (2016) (Fig. 2 A). The effect of IBMX was simulated by decreasing the basal degradation rate of cAMP by PDE.

### Model equations

cAMP activity

(1)

Our numerical model predicts that kiso and kibmx are different for the adult and aged cases (see Results for further details).

$kiso,adult=0.1599⋅[ISO]0.623876.54410.6238+[ISO]0.6238,$
(2)
$kibmx,adult=−0.8730⋅[IBMX]0.83954.05500.8781+[IBMX]0.8395+1.$
(3)

Aged

$kiso,aged=0.1434⋅[ISO]2.0704334.27992.0704+[ISO]2.0704,$
(4)
$kibmx,aged=−0.9118⋅[IBMX]0.878147.67050.8781+[IBMX]0.8781+1,$
(5)
$kcch=0.0146⋅[CCh]1.440251.73311.4402+[CCh]1.4402,$
(6)
$k1=KAC,I+KAC1+exp((KCa−kbCM⋅fCMikfCM⋅(1−fCMi))/KAC,Ca),$
(7)
$k2=kibmx⋅265.3512⋅[cAMP]5.734324.72906.7343+[cAMP]6.7343,$
(8)
$k3=kPKA⋅[cAMP](nPKA−1)kPKA,cAMPnPKA+[cAMP]nPKA.$
(9)

In addition, in the aged-model equations that follow, two additional changes were made to the cAMP activation curve for If and the phospholamban (PLB) sensitivity to PKA phosphorylation (see Results and Fig. 2 C for further details).

Hyperpolarization activated “funny” current, If

$IfNa=0.3833⋅gIf⋅(Vm−ENa)⋅y,$
(10)
(11)

(12)
$Vshift=Kif⋅[cAMP]nifK0.5ifnif+[cAMP]nif−18.1040,$
(13)
$y∞=11+exp((Vm+104.2−Vshift)/16.3).$
(14)

Aged

(15)
$Vshift=Kif⋅[cAMP]nifK0.5ifnif+[cAMP]nif−18.0457,$
(16)
$y∞=11+exp((Vm+110−Vshift)/16.3),$
(17)
$τy=1.5049/(exp(−(Vm+590.3)⋅0.01094)+exp((Vm−85.1)/17.2)).$
(18)

The rate of Ca2+ uptake (pumping) by the SR, jup

$f([PLBp])=2.9102⋅[PLBp]9.55170.27639.5517+[PLBp]9.5517+0.4998,$
(19)
$jup=Pup,basal⋅f([PLBp])1+Kup/[Ca2+]i.$
(20)

Aged

$f([PLBp])=4.0152⋅[PLBp]9.94630.28539.9463+[PLBp]9.9463+0.4999,$
(21)
$jup=0.7∗Pup,basal⋅f([PLBp])1+Kup/[Ca2+]i.$
(22)

T-type Ca2+ current, ICaT

$ICaT=gCaT⋅(Vm−ECaT)⋅dT⋅fT.$
(23)

Aged

$ICaT=0.8⋅gCaT⋅(Vm−ECaT)⋅dT⋅fT,$
(24)
$dT∞=11+exp(−(Vm+26)/6),$
(25)
(26)
$τdT=1/(1.068⋅exp((Vm+26.3)/30)+1.068⋅exp(−(Vm+26.3)/30)),$
(27)
$τfT=1/(0.0153⋅exp(−(Vm+61.7)/83.3)+0.015⋅exp((Vm+61.7)/15.38)).$
(28)

RyR function

$jSRCarel=ks⋅OO⋅([Ca2+]jSR−[Ca2+]sub).$
(29)

Aged

$jSRCarel=0.35⋅ks⋅OO⋅([Ca2+]jSR−[Ca2+]sub),$
(30)
$kCaSR = MaxSR − MaxSR−MinSR1+(EC50SR/[Ca2+]jSR)HSR,$
(31)
$koSRCa=koCa/kCaSR,$
(32)
$kiSRCa=kiCa⋅kCaSR,$
(33)
$dR/dt=(kim⋅RI−kiSRCa⋅[Ca2+]sub⋅R)−(koSRCa⋅[Ca2+]sub2⋅R−kom⋅OO),$
(34)
$dOO/dt=(koSRCa⋅[Ca2+]sub2⋅R−kom⋅OO)−(kiSRCa⋅[Ca2+]sub⋅OO−kim⋅S),$
(35)
$dS/dt=(kiSRCa⋅[Ca2+]sub⋅OO−kim⋅S)−(kom⋅S−koSRCa⋅[Ca2+]sub2⋅RI),$
(36)
$dRI/dt=(kom⋅S−koSRCa⋅[Ca2+]sub2⋅RI)−(kim⋅RI−kiSRCa⋅[Ca2+]sub⋅R).$
(37)

The software was developed in Matlab (MathWorks). Numerical integration was performed using the Matlab ode15s stiff solver, and the model simulations were run for 2,000 s to ensure that steady state was reached. Computation was performed on an Intel Core i7-4790 CPU at 3.60 GHz with 8 GB of RAM.

All model equations and parameters (Tables S1 and S2) appear in the supplement. The source code of the numerical model is available at http://bioelectric-bioenergetic-lab.net.technion.ac.il/programs.

### Online supplemental material

Table S1 shows coupled differential equation variables in the initial conditions. Table S2 shows model constants. The supplemental text shows all model equations.

## Results

### Mechanisms associated with both the M and Ca2+ clock contribute to the age-deteriorated pacemaker function

To explore the relative role of M and Ca2+ clock components in age-associated deterioration of pacemaker function, we used our numerical model. Fig. 3 shows the coupled-clock function of adult mice, aged mice when only changes in the M clock mechanisms were taken into account (Fig. 1 B), aged mice when only changes in the Ca2+ clock mechanisms were taken into account (Fig. 1 C), and aged mice when changes in both M and Ca2+ clock mechanisms were included. The numerical model simulations show that only when mechanisms associated with both the M and Ca2+ clocks are taken into account does the reduction in spontaneous AP firing rate match the experimental results (Larson et al., 2013). Fig. 3 shows that age-reduced pacemaker function was associated with reduced release of Ca2+ from the RyR (reduced jSRCarel), small elevation in SR load (less release compared with sequestration of Ca2+), and reduction in If, INCX, ICa,L, and ICa,T. Quantitatively, the peak amplitude of jSRCarel changed from 0.291 µM/ms in adult cells to 0.20 µM/ms in aged cells, and the amount of Ca2+ ejected over one cycle was 8.14 µM in adult cells and 7.69 µM in aged cells.

Although we changed SERCA activity, HCN current conductivity, and RyR kinetics in accordance with experimental data (Fig. 1), we also performed parameter sensitivity analysis to prove that the model predictions are not dependent on only one set of parameters. Fig. 4 shows the changes in AP firing rate in the aged model in response to changes in V1/2 of the HCN channel or in response to changes in Pup of the SERCA. Moreover, the range we chose for Pup allowed us not only to reproduce the experimental results, but also to provide a stable beating rate (Fig. 4).

### Reduction in INCX, INa, IK, or INaK can deteriorate age-dependent pacemaker function

Although age-dependent changes of RyR, SERCA, PLB, ICa,L, ICa,T, and If have been reported, no data exist on changes of INCX, INa, IK, or INaK. However, the numerical model can illustrate the indirect result of other pacemaker functions caused by the coupled-clock effect. Fig. 5 illustrates all currents for the basal adult and aged cases. INCX amplitude decreases significantly with aging (Fig. 5 J), IKr and INaK density (Fig. 5, D and K) are higher during early DD in aged cells, and the density of INa (Fig. 5, G and H) is lower during DD, but with a peak amplitude more shifted to the early DD. The net effect of all coupled-clock changes is a slower AP upstroke during DD in aged cells than in adult cells, thus leading to a lower AP firing rate.

### Age-associated modification in If sensitivity to cAMP and PLB sensitivity to PKA can explain the restoration of age-deteriorated pacemaker function

We next simulated the changes in spontaneous AP firing rate in response to maximal ISO, IBMX, or cAMP concentration. Fig. 6 shows that in response to any of these drugs (which bring the beating rate to its maximal value), the spontaneous AP firing rate in aged pacemaker cells is lower than in adult cells. However, these results contradict recently published data (Liu et al., 2014; Yaniv et al., 2016). The model was modified to fit the experimental data: we searched for modifications in clock mechanisms that could bring the maximal spontaneous AP firing rate of aged cells to the one simulated for adults. We focused on two targets, age-dependent modification in HCN4 sensitivity to cAMP and PLB sensitivity to PKA. Fig. 6 shows that modification of HCN4 sensitivity to cAMP in aged cells can partially restore the maximal age-deteriorated pacemaker AP firing rate to the adult level, but not to the level documented experimentally. We modified PLB sensitivity to PKA so that the maximal heart rate in response to ISO or IBMX stimulation will fit the experimentally published levels (Liu et al., 2014; Yaniv et al., 2016). Only when the modification in PLB sensitivity to PKA (see Materials and methods, Eqs. 19, 20, 21, and 22) was also taken into account could the maximal AP firing rate of aged pacemaker cells be restored to the adult level. Note that although there is a slight difference between the final predicted maximal heart rate of adult and aged pacemaker cells, this difference is not significant and cannot not be measured experimentally.

Figs. 7 and 8 show how the coupled-clock mechanisms lead to restored maximal pacemaker function. In response to maximal concentration of ISO (Fig. 7), cAMP and PKA are restored, the release of Ca2+ is compensated for, more Ca2+ is accumulated in the SR because the rate of Ca2+ sequestration is higher than the rate of release, and If, INCX, ICa,L, and ICa,T are all completely or partly restored. In response to maximal concentration of IBMX (Fig. 8), similar changes in coupled-clock mechanisms occur.

### Age-associated modification in If sensitivity to cAMP and PLB sensitivity to PKA can explain the reduced sensitivity of pacemaker function to external and internal stimuli

Experimental data show that the sensitivity to either PDE inhibition (by IBMX) or β-AR receptor stimulation (by ISO) is reduced in advanced age (Yaniv et al., 2016). Our numerical model could reproduce these results only when kiso,aged and kibmx,aged (Eqs. 4 and 5) were modified with respect to their adult expression. Thus, the model predicts that in aged cells, the response to ISO or IBMX stimulation is different than in adult pacemaker cells. This age-related modification was generated to match the drug dose–dependent results that were reported experimentally (Yaniv et al., 2016). Fig. 9 shows that higher concentrations of ISO or IBMX are needed to achieve the same percentage of increase in spontaneous AP firing rate. Based on the numerical model results, the age-associated EC50 was ∼203 nM for ISO and 26 µM for IBMX compared with ∼11 nM and ∼6.3 µM in adults in response to ISO and IBMX, respectively. Note that these results can be obtained only by age-associated modification in HCN4 sensitivity to cAMP and PLB sensitivity to PKA.

### The age-associated basal AP firing rate can be reversed

The numerical model was used to predict whether the age-associated decrease in basal beating rate can be reversed using drug perturbations or gene manipulation. Fig. 10 A shows that 167 nM ISO can reverse the aged basal AP firing rate to that of adult mouse pacemaker cells. Fig. 10 B shows that 50 µM IBMX can also reverse the aged basal AP firing rate to that of adult mouse pacemaker cells. Thus, using pharmacological interventions, the old heart can be made adult again. A recent patent (Maltsev et al., 2016) suggests that overexpression of AC1 or AC8 can lead to restoration of normal beating rate. Fig. 11 shows that increasing the presence of AC1 and AC8 (by increasing KAC) reverses the age-associated decrease in beating rate and restores it to the adult level.

## Discussion

Our first major conclusion is that both the M and Ca2+ clocks contribute to the deterioration of basal pacemaker function in advanced age (Fig. 3). This conclusion is in agreement with experimental data from mouse SAN tissue and pacemaker cells, which have shown deterioration in ionic current (decrease in HCN, L-type, and T-type currents; Larson et al., 2013), Ca2+-related mechanisms (decrease in the amount of RyR and SERCA; Liu et al., 2014), and the node signaling that coupled the clocks (Ca2+ global and local release parameters; Liu et al., 2014). Note that only when changes in SERCA, If, and RyR (i.e., both M-related and Ca2+-related clock mechanisms) are taken into account can we simulate both the age-associated decrease in basal mean AP firing rate and the changes in other experimentally measured coupled-clock mechanisms (Fig. 1). Because the two clocks are coupled (Yaniv et al., 2013b), changes in one mechanism related to one clock can indirectly affect mechanisms related to the other clock (Behar et al., 2016). Thus, even when changes in only one clock mechanism are taken into account (either M or Ca2+ in Fig. 3), the clock entrainment leads to indirect changes in the other. Consequently, an explicit change in the numerical model equations for a specific mechanism is not always necessary to describe known experimental modifications observed in that mechanism’s function. For example, the Na+/Ca2+ exchanger (NCX) current decreases in aged simulations as a consequence of changes in other mechanisms related to Ca2+ signaling (Fig. 3). Importantly, entrainment between the clocks means that restoring one mechanism may partially restore the others; thus, even if a particular drug is not known to restore a specific mechanism, it can restore it indirectly by acting on other mechanisms.

Our second major conclusion is that the sensitivity to β-AR stimulation or PDE inhibition is reduced in advanced age (Fig. 9). These results are in accordance with recently published data (Liu et al., 2014; Yaniv et al., 2016). This reduction in sensitivity can be the result of a reduction in the number of receptors or in the sensitivity of each receptor. Thus, a higher dose of drugs is needed to induce the desired response.

Our third major conclusion is that the reduced maximal beating rate can be reversed using maximal ANS stimulation, cAMP application, or PDE inhibition (Figs. 6, 7, and 8). These results are in accordance with recent publications in which the reduced maximal HR of aged pacemaker tissue can be reversed by maximal β-AR stimulation by using ISO or a high level of PDE inhibition using IBMX (Liu et al., 2014; Yaniv et al., 2016). In addition, the maximal AP firing rate could be restored in single isolated pacemaker cells by cAMP (Sharpe et al., 2017). Our model predicts that the age-associated reduction in maximal beating rate can be reversed only if both the relationship between cAMP and HCN4 and the relationship between PKA and PLB phosphorylation are different in aged mice than in adult mice (Figs. 6, 7, and 8). Our model also predicts that, at maximal β-AR stimulation or PDE inhibition, the cAMP/PKA activity of aged pacemaker tissue is similar to that of adults. It is important to note that restoring the maximal beating rate is not equivalent to reversing the aged cell back to adult cell function because, although the AP firing rate was restored (i.e., the end effect “maximal beating rate”), the intrinsic currents and ion channel conductances are not all the same between adult cells and aged cells that were treated with one of these drugs (Figs. 7 and 8). In other words, a restored maximal beating rate is not equivalent to a restoration of the coupled-clock mechanisms to the adult beating rate but only to its end effect.

Our numerical model equations of the membrane and Ca2+ components are in general based on the model of Kharche et al. (2011) (see details in Materials and methods under Basal numerical model). However, to study whether age-deteriorated AP firing is caused by membrane or Ca2+ molecules, a description of the node signaling that couples them is critical. Our numerical model is the first to include a description of AC-cAMP-PKA signaling and updated Ca2+ dynamics based on recent data (Liu et al., 2014). Thus, only our adaptation to the Kharche et al. (2011) model can be used to answer this study’s major question: “What is the relative contribution of each clock to the age-associated deterioration in pacemaker function?” To further elaborate on this question and prove that our conclusions are robust against the choice of the model parameters, we performed parameter sensitivity analysis (Fig. 4). The AP firing rate in the aged model depends on both membrane and Ca2+ clock components when changes in V1/2 of the HCN channel or in response to changes in Pup of the SERCA are taken into account.

Recent experimental data demonstrate the inability to restore the maximal beating rate of isolated mouse pacemaker cells by means of PDE inhibition (Sharpe et al., 2017). However, these data contradict experimental results at the tissue level (Liu et al., 2014; Yaniv et al., 2016) and our numerical model predictions. Moreover, it is well known that heterogeneity among single SAN cells exists (Boyett et al., 2000) and, importantly, that the heterogeneity of the beating rate increases in the aged model (Lowsky et al., 2014), specifically at the single-cell level (Glauche et al., 2011). Furthermore, during isolation, some receptors and molecules can be impaired, especially in the aged model. Thus, the validity of experiments performed on single cells isolated from a small sample of aged mammals is debatable. To strengthen the findings from single cells, data must also be collected at the SAN tissue level. Note, however, that although pacemaker cell heterogeneity exists in the SAN tissue, the neighborhood effect (Michaels et al., 1986) synchronizes all the cells to follow the strongest one (i.e., the cell that beats the fastest).

### Limitations

Our model does not incorporate the previously documented age-associated reduction in the number of NCXs (Liu et al., 2014). Although a decrease in the number of NCXs should result in a lower entrainment level in the coupled-clock model and lead to a decrease in AP firing rate, the numerical model that we used here causes the AP firing rate to increase: reducing the number of NCXs allows more Ca2+ to accumulate in the cytosol and thus to increase the AP firing rate. Thus, we did not decrease the conductance of NCX in the numerical model. However, because of the reduction in SERCA, RyR, and HCN channel function, there is a significant decrease in NCX current, which reduces pacemaker function. This is caused by the coupling between the clocks (see above).

Because of the lack of experimental data on Ca2+/calmodulin-dependent protein kinase II (CaMKII) activity in both adult and aged mice (Maltsev et al., 2014; Yaniv and Maltsev, 2014), we did not take into account age-related changes in CaMKII activity. Future experimental data on CaMKII activity in adult and aged pacemaker cells will allow us to explore its function. Moreover, no data exist for both adult and aged mice on local Ca2+ concentration and cAMP concentrations. Future experimental data will allow us to develop a local numerical model instead of our common pool model.

We showed in Fig. 6 B that the model predicted an 18% decrease in AP firing rate in aged pacemaker cells. Thus, the model prediction is on the lower bound of the experimental data. This suggests that other mechanisms, although playing a more minor role than the ones modeled, may be involved in further reducing the AP firing rate in aged cells.

The model was designed to describe mouse pacemaker function. Because of its short life span, the mouse is the ideal animal model to explore aging. Thus, the large majority of data on age-deteriorated pacemaker function were reported for mice. However, the end effect of aging in humans and mice is the same: in both models, there is no change in the basal beating rate, but only in the intrinsic beating rate (Jose and Collison, 1970). Future data from SAN tissue and signal pacemaker cells from larger aged mammals, including humans, will make it possible to develop computational models for mammals larger than mice.

The effect of PKA on IKs was documented in ventricular and atrial myocytes (Sampson et al., 2008). However, because no experimental data exist for sinoatrial cells, it is ambiguous to model PKA effect on IKs. Therefore, the current numerical model does not include the PKA effect on IKs.

### Clinical aspects

Maintaining the basal heart rate when there is a decrease in the intrinsic beating rate of SAN cells implies that the ANS system must overcompensate for the deteriorated pacemaker function. Chronic stimulation of the β-AR can lead to cardiac hypertrophy and other damage to the heart (Kudej et al., 1997). We showed here that the PDE inhibitor can reverse the age-associated intrinsic beating rate, and thus may be a therapeutic tool if applied specifically only to the SAN area to prevent chronic stimulation of β-AR.

Overexpression of AC1 and/or AC8 has been suggested as a method to restore normal heart rate in the case of bradycardia (Maltsev et al., 2016). We showed here (Fig. 11) that overexpression of Ca2+-activated AC can reverse the age-associated intrinsic beating rate. Thus, our numerical model shows for the first time the feasibility of such an approach for age-associated deficient pacemaker function.

Although we showed here that the maximal intrinsic beating rate can be restored under in vivo conditions, additional mechanisms, such as the preload and afterload, also affect the heart rate. Because the skeletal muscles contribute to the afterload and preload, they, too, affect the heart rate changes. Thus, even if the maximal action potential firing rate can be restored, aging-associated changes in skeletal muscle may prevent the restoration of maximal HR. Moreover, the ability of the heart to reach its maximal rate in response to pharmacological interventions implies that these drugs also restore the energetic balance in the heart (Yaniv et al., 2013a).

## Acknowledgments

The work was supported by the National Natural Science Foundation of China/Israel Science Foundation Joint Research Program, no. 398/14 (Y. Yaniv), the Ministry of Science and Technology, Israel (Y. Yaniv), an Aly-Kaufman Postdoctoral Fellowship (J. Behar), and the Center for Absorption in Science, Ministry of Aliyah and Immigrant Absorption, State of Israel (J. Behar). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The authors declare no competing financial interests.

Author contributions: Y. Yaniv and J. Behar conceived and designed the research; J. Behar performed simulations and analyzed the simulated data; J. Behar and Y. Yaniv drafted the manuscript; J. Behar prepared the figures; and J. Behar and Y. Yaniv edited and revised the manuscript and approved the final version.

José D. Faraldo-Gómez served as editor.

## References

References
Behar
,
J.
, and
Y.
Yaniv
.
2016
.
Dynamics of PKA phosphorylation and gain-of-function in cardiac pacemaker cells: A computational model analysis
.
Am. J. Physiol. Circ. Physiol.
310
:
H1259
H1266
.
Behar
,
J.
,
A.
Ganesan
,
J.
Zhang
, and
Y.
Yaniv
.
2016
.
The autonomic nervous system regulates the heart rate through cAMP-PKA dependent and independent coupled-clock pacemaker cell mechanisms
.
Front. Physiol.
7
:
419
.
Boyett
,
M.R.
,
H.
Honjo
, and
I.
Kodama
.
2000
.
The sinoatrial node, a heterogeneous pacemaker structure
.
Cardiovasc. Res.
47
:
658
687
.
Brodde
,
O.E.
, and
K.
Leineweber
.
2004
.
Autonomic receptor systems in the failing and aging human heart: Similarities and differences
.
Eur. J. Pharmacol.
500
:
167
176
.
Ferrari
,
A.U.
2002
.
Modifications of the cardiovascular system with aging
.
Am. J. Geriatr. Cardiol.
11
:
30
33
.
Glauche
,
I.
,
L.
Thielecke
, and
I.
Roeder
.
2011
.
Cellular aging leads to functional heterogeneity of hematopoietic stem cells: A modeling perspective
.
Aging Cell.
10
:
457
465
.
Hagberg
,
J.M.
,
W.K.
Allen
,
D.R.
Seals
,
B.F.
Hurley
,
A.A.
Ehsani
, and
J.O.
Holloszy
.
1985
.
A hemodynamic comparison of young and older endurance athletes during exercise
.
J. Appl. Physiol.
58
:
2041
2046
.
Heath
,
G.W.
,
J.M.
Hagberg
,
A.A.
Ehsani
, and
J.O.
Holloszy
.
1981
.
A physiological comparison of young and older endurance athletes
.
J. Appl. Physiol.
51
:
634
640
.
Jose
,
A.D.
, and
D.
Collison
.
1970
.
The normal range and determinants of the intrinsic heart rate in man
.
Cardiovasc. Res.
4
:
160
167
.
Kharche
,
S.
,
J.
Yu
,
M.
Lei
, and
H.
Zhang
.
2011
.
A mathematical model of action potentials of mouse sinoatrial node cells with molecular bases
.
Am. J. Physiol. Circ. Physiol.
301
:
H945
H963
.
Kudej
,
R.K.
,
M.
Iwase
,
M.
Uechi
,
D.E.
Vatner
,
N.
Oka
,
Y.
Ishikawa
,
R.P.
Shannon
,
S.P.
Bishop
, and
S.F.
Vatner
.
1997
.
Effects of chronic β-adrenergic receptor stimulation in mice
.
J. Mol. Cell. Cardiol.
29
:
2735
2746
.
Kurata
,
Y.
,
I.
Hisatome
,
S.
Imanishi
, and
T.
Shibamoto
.
2002
.
Dynamical description of sinoatrial node pacemaking: Improved mathematical model for primary pacemaker cell
.
Am. J. Physiol. Circ. Physiol.
283
:
H2074
H2101
.
Kusumoto
,
F.M.
, and
N.
Goldschlager
.
1996
.
Cardiac pacing
.
N. Engl. J. Med.
334
:
89
99
.
Larson
,
E.D.
,
J.R.
St Clair
,
W.A.
Sumner
,
R.A.
Bannister
,
C.
Proenza
,
J.
Clair
, and
W.A.
Sumner
.
2013
.
Depressed pacemaker activity of sinoatrial node myocytes contributes to the age-dependent decline in maximum heart rate
.
110
:
18011
18016
.
Liu
,
J.
,
S.
Sirenko
,
M.
Juhaszova
,
S.J.
Sollott
,
S.
Shukla
,
Y.
Yaniv
, and
E.G.
Lakatta
.
2014
.
Age-associated abnormalities of intrinsic automaticity of sinoatrial nodal cells are linked to deficient cAMP-PKA-Ca2+ signaling
.
Am. J. Physiol. Circ. Physiol.
306
:
H1385
H1397
.
Lowsky
,
D.J.
,
S.J.
Olshansky
,
J.
Bhattacharya
, and
D.P.
Goldman
.
2014
.
Heterogeneity in healthy aging
.
J. Gerontol. A Biol. Sci. Med. Sci.
69
:
640
649
.
Maltsev
,
V.A.
, and
E.G.
Lakatta
.
2009
.
Synergism of coupled subsarcolemmal Ca2+ clocks and sarcolemmal voltage clocks confers robust and flexible pacemaker function in a novel pacemaker cell model
.
Am. J. Physiol. Circ. Physiol.
296
:
H594
H615
.
Maltsev
,
V.A.
,
Y.
Yaniv
,
A.V.
Maltsev
,
M.D.
Stern
, and
E.G.
Lakatta
.
2014
.
Modern perspectives on numerical modeling of cardiac pacemaker cell
.
J. Pharmacol. Sci.
125
:
6
38
.
Maltsev
,
V.
,
E.
Lakatta
,
I.
Zahanich
, and
S.
Sirenko
.
2016
. Engineered biological pacemakers. US Patent No. 9506032. 29 November 2016.
Michaels
,
D.C.
,
E.P.
Matyas
, and
J.
Jalife
.
1986
.
Dynamic interactions and mutual synchronization of sinoatrial node pacemaker cells. A mathematical model
.
Circ. Res.
58
:
706
720
.
Ogawa
,
T.
,
R.J.
Spina
,
W.H.
Martin
III
,
W.M.
Kohrt
,
K.B.
Schechtman
,
J.O.
Holloszy
, and
A.A.
Ehsani
.
1992
.
Effects of aging, sex, and physical training on cardiovascular responses to exercise
.
Circulation.
86
:
494
503
.
Sampson
,
K.J.
,
C.
Terrenoire
,
D.O.
Cervantes
,
R.A.
Kaba
,
N.S.
Peters
, and
R.S.
Kass
.
2008
.
Adrenergic regulation of a key cardiac potassium channel can contribute to atrial fibrillation: Evidence from an I Ks transgenic mouse
.
J. Physiol.
586
:
627
637
.
Sharpe
,
E.J.
,
E.D.
Larson
, and
C.
Proenza
.
2017
.
Cyclic AMP reverses the effects of aging on pacemaker activity and If in sinoatrial node myocytes
.
J. Gen. Physiol.
149
:
237
247
.
Yaniv
,
Y.
, and
V.A.
Maltsev
.
2014
.
Numerical modeling calcium and CaMKII effects in the SA node
.
Front. Pharmacol.
5
:
58
.
Yaniv
,
Y.
,
M.
Juhaszova
, and
S.J.
Sollott
.
2013
a
.
Age-related changes of myocardial ATP supply and demand mechanisms
.
Trends Endocrinol. Metab.
24
:
495
505
.
Yaniv
,
Y.
,
S.
Sirenko
,
B.D.
Ziman
,
H.A.
Spurgeon
,
V.A.
Maltsev
, and
E.G.
Lakatta
.
2013
b
.
New evidence for coupled clock regulation of the normal automaticity of sinoatrial nodal pacemaker cells: Bradycardic effects of ivabradine are linked to suppression of intracellular Ca2+ cycling
.
J. Mol. Cell. Cardiol.
62
:
80
89
.
Yaniv
,
Y.
,
A.
Ganesan
,
D.
Yang
,
B.D.
Ziman
,
A.E.
Lyashkov
,
A.
Levchenko
,
J.
Zhang
, and
E.G.
Lakatta
.
2015
a
.
Real-time relationship between PKA biochemical signal network dynamics and increased action potential firing rate in heart pacemaker cells: Kinetics of PKA activation in heart pacemaker cells
.
J. Mol. Cell. Cardiol.
86
:
168
178
.
Yaniv
,
Y.
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2015
b
.
From two competing oscillators to one coupled-clock pacemaker cell system
.
Front. Physiol.
6
:
28
.
Yaniv
,
Y.
,
I.
Ahmet
,
K.
Tsutsui
,
J.
Behar
,
J.M.
Moen
,
Y.
Okamoto
,
T.R.
Guiriba
,
J.
Liu
,
R.
Bychkov
, and
E.G.
Lakatta
.
2016
.
Deterioration of autonomic neuronal receptor signaling and mechanisms intrinsic to heart pacemaker cells contribute to age-associated alterations in heart rate variability in vivo
.
Aging Cell.
15
:
716
724
.

Abbreviations used:

• AC

•
• ANS

autonomic nervous system

•
• AP

action potential

•
• β-AR

•
• CaMKII

Ca2+/calmodulin-dependent protein kinase II

•
• DD

diastolic depolarization

•
• HCN channel

hyperpolarization-activated cyclic nucleotide–gated channel

•
• HR

heart rate

•
• IBMX

3-isobutyl-1-methylxanthine

•
• ISO

isoproterenol

•
• LCR

local Ca2+ release

•
• NCX

Na+/Ca2+ exchanger

•
• PDE

phosphodiesterase

•
• PLB

phospholamban

•
• SAN

sinoatrial node

•
• SERCA

sarcoplasmic reticulum Ca2+ ATPase