http://orcid.org/0000-0002-6210-1523
Figures
Abstract
The accumulation of Alzheimer’s disease (AD) associated Amyloid beta (Aβ) oligomers can trigger aberrant intracellular calcium (Ca2+) levels by disrupting the intrinsic Ca2+ regulatory mechanism within cells. These disruptions can cause changes in homeostasis levels that can have detrimental effects on cell function and survival. Although studies have shown that Aβ can interfere with various Ca2+ fluxes, the complexity of these interactions remains elusive. We have constructed a mathematical model that simulates Ca2+ patterns under the influence of Aβ. Our simulations shows that Aβ can increase regions of mixed-mode oscillations leading to aberrant signals under various conditions. We investigate how Aβ affects individual flux contributions through inositol triphosphate (IP3) receptors, ryanodine receptors, and membrane pores. We demonstrate that controlling for the ryanodine receptor’s maximal kinetic reaction rate may provide a biophysical way of managing aberrant Ca2+ signals. The influence of a dynamic model for IP3 production is also investigated under various conditions as well as the impact of changes in membrane potential. Our model is one of the first to investigate the effects of Aβ on a variety of cellular mechanisms providing a base modeling scheme from which further studies can draw on to better understand Ca2+ regulation in an AD environment.
Citation: Latulippe J, Lotito D, Murby D (2018) A mathematical model for the effects of amyloid beta on intracellular calcium. PLoS ONE 13(8): e0202503. https://doi.org/10.1371/journal.pone.0202503
Editor: Laszlo Csernoch, University of Debrecen, HUNGARY
Received: November 6, 2017; Accepted: August 3, 2018; Published: August 22, 2018
This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Data Availability: All relevant data are within the paper.
Funding: Research reported in this project was supported by an Institutional Development Award (IDeA) from the National Institute of General Medical Sciences of the National Institutes of Health under grant number P20GM103449. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of NIGMS or NIH. Part of this study was also funded by a Faculty Development grant from Norwich University.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Alzheimer’s Disease (AD) is a devastating neurodegenerative illness affecting over 40 million people worldwide. AD is the leading cause of dementia and is characterized by a progressive and irreversible decline in memory and cognitive skills [1]. The prevalence of AD and associated dementia is estimated to double in the next 20 years, and as such, there is a critical need to better understand this disease. While the appearance of extracellular hydrophobic amyloid plaques and intracellular neurofibrillar tangles associated with tau proteins have become the hallmarks of the disease, the cause of AD remains unknown. Abnormal intracellular Ca2+ levels have been observed in AD brains even before the presentation of clinical symptoms and amyloid plaques [2, 3]. Because it has been shown that Amyloid beta peptide (Aβ) accumulation can lead to increased intracellular Ca2+ levels [4, 5], studying its effect on intracellular mechanisms is important for understanding its impact on neuronal functions. Intracellular accumulation of Aβ can cause an increase release of Ca2+ from internal stores such as the Endoplasmic Reticulum (ER) [4, 6–8]. As such, Aβ may lead to both local and global Ca2+ proliferation that are persistent and cytotoxic. Sustained Ca2+ disregulation can trigger apoptosis leading to premature neuronal death, a characteristic feature seen in AD. Although the accumulation of Aβ has been linked to the progression of AD by altering Ca2+ signaling processes within neurons and neuroglia, the mechanisms for how and why this occurs are not fully understood.
Our goal is to use mathematical modeling to describe various conditions under which Aβ can lead to aberrant Ca2+ signaling. The amyloid hypothesis suggests that the accumulation of Aβ in the brain is the primary driving force of AD pathogenesis [9–12]. In this hypothesis, the formation of Aβ plaques and fibrils are a consequence of the imbalance between the formation and sequestration of Aβ. The slow accumulation of Aβ peptides can alter Ca2+ signaling processes leading to synaptic failure and neuronal death. Although it is unclear how Aβ disrupts intracellular Ca2+ homeostasis, there is growing evidence that Aβ directly affects the production of inositol triphosphate (IP3) [7], calcium-induced calcium release (CICR) through the ryanodine receptor (RyR) [13, 14], and the plasma membrane [15, 16]. We use the results of these works to make simplifying assumptions for how Aβ affects various Ca2+ signaling mechanisms in a simplified whole-cell model.
In this study we present a theoretical approach to better understand the driving mechanisms for various Ca2+ oscillatory patterns within an AD environment. By developing a mathematical model for intracellular Ca2+ regulation, we can begin to study how Aβ affects Ca2+ flux through various individual channels and pumps. Investigating model solutions can also provide important information on the impact of Aβ on Ca2+ basal levels over various timescales. Due to current experimental limitations, mathematical and computational models can provide insights for targeting specific mechanisms in order to restore neuronal function, and to suggest symptomatic improvement strategies. In fact, according to Liang et al. (2015) there is growing momentum to study Ca2+ dynamics in AD, specifically through computational and mathematical modeling, and this work outlines such an approach.
Methods
Calcium model formulation
In order to study the effects of Aβ on intracellular Ca2+, we first build a simplified whole-cell Ca2+ model by making use of the vast array of work on modeling Ca2+ dynamics and the calcium signaling “toolkit” (see [17–20] for example). Once this model is developed, we add the influence of Aβ by altering various components of the model. We then analyze model solutions by investigating the dynamical structure for various parameter regimes in order to draw out conditions that lead to changes in basal Ca2+ levels and aberrant signals. For our purposes, we characterize aberrant signals as non-periodic oscillations over a timescale of about 100-200 seconds.
We model Ca2+ dynamics using traditional methods by tracking the flux in and out of the cytoplasm. Let c denote the concentration of free Ca2+ ions in the cell cytoplasm, then the rate of change in intracellular Ca2+ is governed by
(1)
where
denotes flux. We assume a spatially homogeneous cell whose volume is fixed and track Ca2+ concentration changes in time. As such, the ER and cytoplasm coexist at each point of the cell. Although these simplifying assumptions make the model limited, such an approach has been extremely useful in quantifying and identifying key mechanisms behind certain Ca2+ signaling patterns. We take this position to study the influence of Aβ on Ca2+ dynamics using a simplified whole-cell model. Many different mathematical approaches have been used to better understand Ca2+-mediated neuroglia function (see [21] for an overview of these) and we have utilized this body of work to develop our model below.
Our model structure was selected to account for the major components that have been shown to be influenced by Aβ and that can lead to Ca2+ dynamics on the timescale of seconds. We assume that intracellular Ca2+ in-fluxes (into the cytoplasm) are those corresponding to IP3 receptors (IPR), RyRs, a general membrane leak Jin, and we include a fast Voltage Gated Ca2+ Channel (VGCC) in Jvca to account for membrane permeability. The out-fluxes (out of the cell or sequestered into Ca2+ pools such as the ER) are modeled using a plasma membrane pump Jpm, and a sarco-endoplasmic reticulum Ca2+ ATPase (SERCA) pump JSERCA. A diagram of the major fluxes outlined in the model is given in Fig 1. Also included in the diagram are the model assumptions of the interaction of Aβ with respect to individual flux term.
This figure shows the critical flux terms utilized in the formulation of the model. In addition, the relevant fluxes affected by Aβ are highlighted.
We assume that the ER is homogeneously distributed throughout the interior of the cell and that the fluxes of Ca2+ through the IPR and RyR are proportional to the difference in the concentration of Ca2+ between the ER and the cytoplasm. Under these standard assumptions, using conservation of fluxes, our general Ca2+ model takes the form
(2)
(3)
where ce denote the concentration of Ca2+ in the ER, and γ is the ratio of cytoplasmic volume to the ER volume. The individual contributions of each flux can vary from a simple Hill function to more complicated forms involving numerous parameters and additional terms. We provide a description of each of the flux terms we use in our model below. We also provide background information on their development and why it may be well suited for our purposes. Each term was chosen in order to balance meaningful biophysical quantities while maintaining a tractable mathematical structure.
The term Jvca in (2) links membrane potentials with other Ca2+ signaling mechanisms. Ca2+ regulation and changes in membrane potential do depend on each other and a model that connects their influence may be critical for advancing our understanding of the long term effects of Aβ on Ca2+ signaling. To address this, we include a section in our results that incorporates the effects of membrane potentials into our model and provide some examples of the impact of Aβ on the dynamics. Although a full exploration of the role of membrane potential on Ca2+ dynamics is beyond the scope of the current study, we do provide the base structure from which to build a more complete whole-cell model.
Ca2+ regulation in non-excitable cells, such as astrocytes, is extremely complex and influenced not only by membrane potential but also Ca2+ buffering. Cytoplasmic Ca2+ buffering plays a significant role in calcium’s ability to move throughout the cell. Almost all of the available Ca2+ is bound to buffers and free cytoplasmic Ca2+ cannot move very far before being bound [18]). Because Ca2+ diffusion plays an important role in spatiotemporal signaling, nonlinear buffering may provide some insights behind certain types of oscillatory patterns and Ca2+ waves. However, compared to the temporal release and uptake of Ca2+ into internal stores, buffering occurs on a much faster timescale. In our model we have assumed that Ca2+ buffering is fast, immobile, and has low affinity. Using these assumptions, we have scaled the model to account for fast and linear buffers even though an explicit description of buffering is not provided (further details on buffering can be found in [18]).
IP3 receptor model.
IP3 is a ligand produced by phospholipase C (PLC) when activated by a G-protein coupled receptor reacting to an external stimulus such as neurotransmitters, hormones, and growth factors [22]. As a key secondary messenger, IP3 regulates many important cellular functions, including the release of Ca2+ from the ER through the IP3 receptor [23]. We assume that the flux from the IP3 receptor follows a saturating binding rate model of the form found in [24, 25]. Thus, we write
(4)
where kf controls the density of the IP3 receptor, Jer is a leak from the ER to the cytoplasm, and P0 is the open probability of the IPR. The leak term is necessary to balance the ATPase flux at steady state. To model P0, we use a version of the Sneyd and Dufour (2002) model that is based on previous models found in [26–28]. The model assumes that the receptor can be in one of six states, with R the resting state, O the open state, A the active state, S the shut state, and I1 and I2 both represent inactive states. The transitions between states can depend on Ca2+ and IP3 concentration p. As such, the IP3 receptor will be in various states depending on the value of p.
In this model, the receptor open probability is given by
(5)
where α1 and α2 are parameters that control the individual contributions of the open and active states of the receptor. As in [24], we assume that α1 = 0.1 and α2 = 0.9. In the model, it is assumed that the receptor has four identical and independent subunits and that Ca2+ flows when all four subunits are in either the O or A state. The equations for the transitions between the states are given in [24], but have also been included in the Supplementary Appendix for completeness. This model was chosen since it does respond reasonably well to changes in Ca2+ concentration and IP3 concentrations [25].
RyR model.
To model the contribution of the RyR, we utilize the algebraic model of Friel (1995). The receptor is modeled as a simple leak channel, with a flux through the channel proportional to the concentration difference between the ER and cytosol. This model was used to investigate Ca2+ oscillations in a sympathetic neuron. Thus, the flux through the RyR is given by
(6)
To more accurately model CICR, the rate constant k3 is defined as an increasing sigmoidal function of intracellular Ca2+ concentration and takes the form
(7)
where k1, k2, and kd are parameters. We use n = 3 to match the experimental results obtained in [29]. The parameter k1 in (7) is the zero calcium concentration level leak. This term is often used to ensure a physiologically significant resting Ca2+ level [18]. The term kd corresponds to the RyR channel sensitivity for CICR, and k2 is the maximal rate of the channel. The simple nature of the Friel model makes it a viable choice especially since data for the contributions of Ca2+ flux through the RyR in the presence of Aβ are minimal.
Leak, membrane pump, and SERCA.
The membrane leak Jin is modeled using a linearly increasing function of IP3 concentration [30]. Although this increase may be due to various mechanisms, here we only include a linearly increasing contribution to make sure that steady-state Ca2+ concentration depend on p. Thus,
(8)
where a1 and a2 are parameters.
When modeling Ca2+ dynamics both the SERCA pump and the plasma membrane pump play an important role in maintaining concentration gradients. Different models of these pumps exist and vary in complexity. Here we model the plasma membrane pump using a simple Hill equation of the form
(9)
where Vpm is the maximal velocity and Kpm is the channel sensitivity. We have chosen a standard Hill coefficient of 2 as is found in [31].
To describe the behavior of the SERCA pump we utilize a simplified four state bidirectional Markov model (as found in [20]) of the form
(10)
to model the SERCA pump. In (10), Ki for i = 1, …, 5 are constants. Notice that this choice of model is more complicated than a simple Hill equation. Experimental evidence show that the rate of the pump may be modulated by the level of Ca2+ in the ER [32, 33], and this model provides a way to account for both c and ce.
Voltage dependent calcium channel.
The flux due to changes in membrane potential is given in our model by Jvca. When the volume of the cell is constant, Ca2+ flux and current are related by the equation
(11)
where Ica is the current through the VGCC, F is Faraday’s constant, and w is the cell volume. Since Ica depends on the membrane potential V and Ca2+, we need a way to track membrane potential. Here, we assume that the membrane potential follows a Hodgkin and Huxley formulation [34] that satisfies
(12)
where Cm is the cell capacitance, V is the membrane potential, and the sum is over all ionic currents across the cell membrane. In the case of non-excitable astrocytes, current flow across its membrane is characterized by potassium (K+)-selective membrane conductance [35] and voltage dependent channels [36, 37]. As such, our membrane equation takes the form
(13)
where Ikir, Ina, Il, and Ica correspond to an inward rectifying potassium, sodium, leak, and Ca2+ current, respectively, and Iapp is an applied current. In (13), we utilize a formulation for Ikir(V) similar to [38] and write
(14)
where gkir is the maximal channel conductivity, Vka is the Nernst K+ potential, K0 is the extracellular K+ concentration, and Va1, Va2 and Va3 are constants. In the model of [38], K0 is time dependent and depends on proximal neuronal activity (further details of this formulation can be found in [38]). Here, we simplify this and assume that the external compartment is saturated with K+ and as such consider K0 to be constant. Each of the remaining currents have the form Ix = gx(V − Vx), where gx is the conductance and Vx is the Nernst potential of each channel.
Although many alternative formulations for Ica exist (such as the Goldman-Hodgkin-Katz current approximation), we use a T-type like channel (as that found in [39]) of the form
(15)
where
is the maximal conductance, Vca is the Ca2+ Nernst potential, and mcaT and hcaT are gating variables similar to those used for the Na+ channel in the Hodgkin and Huxley model. A T-type Ca2+ current has a low-threshold activation with an inactivating gating variable that exhibit sub-threshold oscillations and inhibitory rebound bursts. Although other types of Ca2+ currents (such as L-type, Ca2+-dependent potassium channels, etc.) can also be included in the model, here we simply illustrate how one can link membrane potential with other Ca2+ regulatory mechanisms together into a single model and use that model to investigate the role of Aβ on Ca2+ dynamics. As such, the T-type Ca2+ channel described in (15) is sufficient for studying how changes in membrane potentials could impact intracellular Ca2+. Full details of the membrane potential model are provided in the Supplementary Appendix.
Amyloid beta assumptions
The influence of Aβ on individual pumps, channels, and exchangers remains largely unknown. However, some studies provide insight on the effects of Aβ on Ca2+ regulation and we utilize these findings in making our assumptions about the effects of Aβ. For example, De Caluwé and Dupont (2013) formalized a theoretical model to describe a possible feedback loop between Ca2+ and Aβ. Their model showed that the existence of a bistable steady-state region can exists in the presence Aβ. The model suggests that over time, cytosolic Ca2+ proliferation as a result of Aβ could trigger the onset of AD. In this study, we are interested in the generation of aberrant Ca2+ signals on a relatively short timescale. Since the accumulation of Aβ can occur over months, years, and even decades, Aβ concentrations changes occur on a very long timescale compared to changes in Ca2+. As such, we assume that Aβ may be present in the environment, but make no attempt to track changes in Aβ concentration over time.
The formation of Aβ plasma membrane pores can alter Ca2+ signaling by creating additional influx into the cytoplasm [15, 16, 40, 41]. In order to incorporate the possible influence of Aβ generated plasma membrane pores, we use a similar mechanism as found in [42]. Let a represent a fixed level of Aβ concentration present in the environment. Then, we include the term kβ am in Jin and write an altered membrane leak flux as
(16)
where m represents a cooperativity coefficient, and kβ is a rate constant (see [42] for details).
Although it is well known that Aβ can disrupt RyR-regulated Ca2+ signals, the mechanisms for how this happens remains controversial [43–45]. Several studies have addressed the role of RyR-regulated Ca2+ disruptions in AD [46–49]. Paula Lima et al. (2011) report that Aβ can generate prolonged Ca2+ signals in vitro through RyR in primary hippocampal neurons of rat embryos as a result of NMDAR-dependent Ca2+ entry through the plasma membrane. These results agree with studies showing that Aβ can cause substantial Ca2+ influx through NMDAR [50, 51]. In our model these types of fluxes are combined and accounted for in Jin, although in a somewhat crude fashion. Futhermore, Aβ can increase RyR channel open probability on a short timescale [4, 14]. In order to account for this, we alter the RyR channel sensitivity term and assume that it is affected by Aβ. Thus, our altered RyR model takes the form
(17)
where kα is positive and controls the strength of the influence of Aβ. Notice that an increase in a corresponds to an increase in the RyR sensitivity.
Calcium model formulation with amyloid beta.
The simplified whole-cell Ca2+ model described above tracks changes in cytosolic Ca2+ concentration as a function of time in the presence of Aβ. By putting the different channels, exchangers, and pumps together, our Ca2+ model takes the form
(18)
(19)
Notice that this model has twelve dynamic equations, with five of those controlling the IP3 receptor (R, O, A, I1, and I2), three controlling the membrane potential (V, m, and h), and an additional two variables controlling the gating variables mcaT and hcaT (a description of these additional equations are provided in the Supplementary Appendix). We have also added a scaling and control parameter ps to simplify the contribution from the VGCC flux. When ps = 0, the impact of membrane potential on Ca2+ signaling are ignored. When ps > 0, changes in membrane potential do influence the dynamics of intracellular Ca2+. The large set of parameter values was selected to closely match those of previous studies whenever possible and can be found in Table 1.
In the formulation of our model we have assumed that the presence of Aβ can increase membrane permeability by creating plasma membrane Aβ pores, and have used the work of [42] to alter our membrane in-flux. Additionally, we have assumed that the presence of Aβ increases RyR sensitivity and have incorporated additional terms in the RyR component of the model. Currently in our model, IP3 acts as the primary agonist that can then trigger Ca2+ release from IP3 receptors. In addition to affecting the IP3 receptor open probability, Ca2+ release activates RyR and subsequent CICR. In the next chapter we separate our analysis for the model under three assumptions. We first characterize the solutions of the model for constant levels of IP3 and when the effects of membrane potential are ignored. Second, also in the absence of membrane potentials, we expand our model to include a dynamic variable for IP3 production. Lastly, we include the effects of membrane potential on Ca2+ dynamics by incorporating flux through VGCC under various levels of fixed IP3. We further look to characterize the impact of VGCC on model solutions in the presence of Aβ.
Numerical solutions of the differential equations systems were obtained using explicit Runge-Kutta (4,5) algorithms implemented in Matlab [52]. Bifurcation analyses were done using the AUTO [53] extension of XPPAUT [54]. AUTO is a software program used to analyze the bifurcation structures of systems of ordinary differential equations such as (18) and (19). In the bifurcation diagram illustrated in this manuscript, the bifurcation parameter is plotted on the x-axis and the projection of stable and unstable steady-state solutions, upper and lower branches of periodic orbits, and critical transitions points (such as Hopf and period doubling) are labeled accordingly. Hopf bifurcations occur when a steady-state solution loses it’s stability as a result of small changes in a parameter. Mathematically, this occurs when complex eigenvalues (of the linearization) cross the imaginary axis away from the origin. Period doubling bifurcations occur when a small change in a bifurcation parameter causes the period of an oscillatory system to double. These and other types of bifurcations can be used to describe systemic changes in dynamics and in our case, transitions from single-mode oscillations to mixed-mode oscillations (MMOs). The program AUTO allows us to readily compute the location of these bifurcations as model parameters are varied. Initial conditions for all simulations were set at c0 = 0.05 and ce = 10 with R0 = 1 and all other IP3 receptor initial conditions are set to zero. Initial conditions for the membrane potential are V0 = −65, m0 = 0.05, n0 = 0.32, h0 = 0.6, mcaT0 = 0.29, and hcaT0 = 0.01.
Results
In all three sections below, the model dynamics show that aberrant Ca2+ can emerge in the presence of Aβ. These aberrant signals can occur under various conditions suggesting a complex link between the influence of Aβ and the model components. As such, we break down the dynamics of the model by tracking solutions as a result of altering one or two parameters within a specific signaling component. We first look at model solutions when IP3 concentration is fixed and persists and is not influenced by membrane potential. We then allow for IP3 to change dynamically and simulate the response in Ca2+ with and without the influence of Aβ. Lastly, we incorporate the influence of membrane potentials and investigate model solutions for various levels of Aβ when IP3 concentration is fixed.
Calcium model with constant IP3 and no membrane potential
Simulations of (18) and (19) with ps = 0 show that incorporating Aβ in a simplified Ca2+ model can lead to aberrant Ca2+ signaling through dynamical transitions into MMOs. Particularly, we look to track the influence of Aβ on the location of Hopf points (labeled HB), period doubling points (labeled PD), and regions of MMOs (shaded) in order to better understand Ca2+ steady-states and oscillatory patterns. From a phenomenological perspective, a transition through MMO may explain why Aβ can trigger persistent aberrant Ca2+ signals. Although the long-term accumulation and influence of Aβ remains difficult to model, our approach illustrates that aberrant Ca2+ signals can occur under the influence of Aβ for a variety of parameter regimes.
IP3 influence on model dynamics.
Experimentally, IP3 can be photoreleased simultaneously through out a cell. In such conditions, IP3 diffusion is minimized and can be treated as constant. By varying the amount of IP3 available in the cytosol, the model can exhibit Ca2+ oscillations as typically found in various cell types. These oscillatory patterns are critical in order for cells to maintain appropriate concentration gradients and to re-establish homeostasis levels after a triggering event. In the absence of Aβ, model Ca2+ oscillations appear and disappear as a result of transitions through Hopf bifurcations as the parameter p increases. Although these stable oscillatory patterns are predictable for a large domain of the parameter p, the model does exhibit MMOs for small parameter regions near the Hopf points. It is these MMO patterns that we are most interested in. Our results show that in the presence of Aβ, the regions of MMOs can grow and become larger. This leads to an increased possibility for aberrant Ca2+ signals. We hypothesize that aberrant Ca2+ signals observed experimentally can be described as dynamic transitions through MMOs. The complicated patterns in MMOs are often characterized by sub-threshold oscillations interspersed within large relaxation types of oscillations [55–57].
Illustrated in Fig 2A and 2B are numerical simulations of the model, with a = 0 that exhibit sustained Ca2+ oscillations as the amount of p is increased from p = 5 to p = 10. An increase in frequency occurs as p increases within the predictable region and can be seen by comparing Fig 2A and 2B. In Fig 2C, a partial bifurcation diagram showing the stable maximum and minimum periodic amplitudes is given along with the regions of MMOs and the critical transitions points. Fig 2D shows the solution of the model for p = 18.5. For this p-value, the model undergoes MMOs with varying oscillatory patterns of mixed amplitude, but we do not characterize it as aberrant Ca2+ signaling since it is periodic. Notice that no Ca2+ oscillations occur for small and large values of p. This is due to the assumption that the IP3 receptor is unlikely to open at high and low concentrations of IP3. Although our model does not exhibit a typical Hopf bubble, it has the overall qualitative patterns of a robust Ca2+ model as described in [18].
Properties of the Ca2+ model when IP3 concentrations are fixed and no Aβ is present (i.e., a = 0). A and B show the numerical solution with p = 5 and p = 10, respectively. For these parameter values, the frequency of oscillations increases as p increases. C illustrates the partial bifurcation diagram for the model with p as the bifurcation parameter. The middle curve crossing the diagram correspond to the steady-state values (solid for stable, dashed for unstable). Also shown are the maximum and minimum amplitudes of the periodic orbits of the model (solid curves above and below the steady-state curve). Key bifurcation points are labeled HBp1 and HBp2 (Hopf), and PD (period-doubling). The shaded area corresponds to the regions where MMOs are present. D shows MMOs for p = 18.5.
The presence of amyloid beta.
We consider an AD environment where some accumulation of Aβ has occurred but assume that this amount is fixed during the timescale of our simulations. As such, we assume that a is constant and that the value correlates to the concentration of Aβ in μM. That is, the value of a in the model may be thought of as reflecting the stage of the progression of the disease. For example, a small value of a corresponds to a small amount of Aβ that may have accumulated in the early stages of AD, while a larger value of a may correspond to a late stage AD. As the parameter a is varied, we can then study how the dynamics of the model change as the level of Aβ changes. Our ultimate goal is to better understand the implications of various levels of Aβ on Ca2+ so that we may gain knowledge about the possible progression of AD.
The effects of Aβ on calcium steady-state levels.
The accumulation of Aβ can lead to increased steady-state levels as well as large amplitude oscillations even in the absence of IP3 signaling. To better understand the effects of Aβ on model solutions, we simulate the model and look at the effects of changing the parameter a on steady-state levels. Various experimental studies have used Aβ levels of 0.5 and 1.0 [13, 50, 51]. As such, we consider a range of values for a that matches these types of levels.
In the case when no IP3 is present, the accumulation of Aβ increases the steady-state of Ca2+ and can even lead to large-amplitude Ca2+ oscillations. Illustrated in Fig 3A is a bifurcation diagram generated by changing the parameter a when p = 0. Notice that the steady-state value slowly increases as a increases until solutions transition into stable periodic oscillations between the two Hopf points HBa1 and HBa2. As a increases towards the value of 1.293, the steady-state asymptotes and solutions quickly become non-physical. Fig 3B shows a solution where large amplitude oscillations exist for a = 1.15. Fig 3C shows the solution for a = 1.276. Notice that in this case, the solution first undergoes a large Ca2+ increase before settling in on a steady-state value. When the model is presented with levels of Aβ greater than a = 1.29, solutions become intractable and are not analyzed.
This figure shows the effects of Aβ on the steady-state Ca2+ levels in the absence of IP3. A shows a bifurcation diagram with a as the bifurcation parameter. Two Hopf bifurcations, labeled HBa1 and HBa2, give rise to a Hopf bubble with relatively large oscillation amplitude. The steady-state level quickly becomes unphysical as the amount of Aβ increases towards a = 1.3. B shows stable Ca2+ oscillations for a = 1.15. C shows the solution when a = 1.276.
Two Parameter Analysis.
Since both a and p can affect the dynamics of the model, we look at their effects together. Given in Fig 4A is a partial two-parameter bifurcation set for the whole-cell model (18) and (19) when ps = 0 using p and a as the bifurcation parameters. The solid lines correspond to the Hopf bifurcation manifolds, while the dashed lines correspond to the period-doubling manifolds. The shaded area in Fig 4A corresponds to the parameter values that elicit MMOs. Notice that in addition to the large MMO region between the dashed lines, there is a small thin region near the Hopf bifurcation manifold labeled HB2 m. These MMOs vanish for values of a > 0.275. This bifurcation diagram has a number of interesting features that can bring to light important behavior in Ca2+ regulation. For example, when 0.415 < a < 0.472 the bifurcation diagram has four Hopf points. This leads to small amplitude oscillations for certain ranges in the parameter p. Fig 4B shows the bifurcation diagram when a = 0.45. In this diagram, three of the Hopf bifurcations have been labeled as HB1 since all three points lie on the Hopf manifold labeled HB1m in Fig 4A.
A shows a two parameter bifurcation diagram when p and a are varied. The solid curves in the diagram correspond to the Hopf bifurcation manifolds and are labeled HB1m and HB2m. The dashed lines correspond to manifolds of period doubling points. The shaded regions (between the dashed lines and near the bottom of HB2 m) correspond to regions where MMOs occur. B shows a bifurcation diagram for a = 0.45. Notice that for this value of a, four Hopf bifurcations exist and three of them have been labeled with HB1 while the last is labeled HB2. This figure also illustrates the intermediate region of MMOs between the labels PD1 and PD2.
To help us understand this bifurcation structure, various solutions are plotted in Fig 5 and compared to the bifurcation structure given in Fig 4. In Fig 5, a = 0.45 is kept fixed while various values of p are used. The corresponding model solutions show different types of behaviors ranging from steady-state Ca2+ levels, stable periodic solutions, aberrant Ca2+ signals, and MMOs. Illustrated in Fig 5A is solution that exhibits small amplitude oscillations with p = 5. This solution corresponds to the region where the small Hopf bubble exists as shown in Fig 4B. Fig 5B shows a solution that reaches a stable steady-state when p = 20. Fig 5C and 5D show two different solutions that exhibit MMOs for p = 26 and p = 45.5, respectively. Fig 5D shows aberrant Ca2+ signals when p = 45.8. Notice that this solution enters aberrant oscillations as the bifurcation structure transitions through period doubling points within the MMO region. Fig 5D shows sustained Ca2+ oscillations when p = 50.
This figure shows various solutions when a = 0.45. A shows small amplitude oscillations when p = 5 while B shows a stable-steady state solution when p = 20. C and D both show MMOs with multiple sub-threshold oscillations when p = 26 and p = 45.5, respectively. E shows aberrant Ca2+ signals when p = 45.8. F shows sustained Ca2+ oscillations when p = 50.
The results of the two parameter bifurcation analysis seem to suggest that Aβ not only increases the steady-state level, but also influences the regions of MMOs. Particularly, the internal range where MMOs emerge between the Hopf manifolds increases for a large parameter range of a. As such, the model suggests that Aβ may drive aberrant Ca2+ signals as transitions through MMOs. However, as the simulated amount of Aβ increases towards a = 1, the region of MMOs decreases and stable periodic orbits exist for most of the p parameter range. This condition is illustrated in Fig 6 where two solutions are shown for relatively large values of a. Fig 6A shows the solution for a = 1 with p = 30, while B shows the solution for a = 1.2 and p = 20. In these cases, Ca2+ oscillations have much larger amplitudes than in the absence of Aβ.
This figures shows two examples of sustained Ca2+ oscillations for two sets of parameter selections of a and p. A shows Ca2+ oscillations with peak amplitudes around 2 corresponding to a = 1 and p = 30. B shows similar oscillations with peak amplitude closer to 3 corresponding to a = 1.2 and p = 20.
Contribution of Aβ on calcium signaling through the ryanodine receptor.
To better understand the influence of Aβ on Ca2+ signaling, we now turn our attention to the contributions of the RyR. Recall that in our model we assume that Aβ affects the RyR by altering the receptor’s sensitivity for CICR. To determine the specific contribution of this altered sensitivity, we first simulate the model for fixed a and p and describe the dynamics of varying the parameter ka. Our simulations show that in the presence of Aβ, changing the parameter kα can produce aberrant Ca2+ signals. These aberrant signals result as transitions through MMOs for fixed Aβ levels. However, we also show that we can control these signals by increasing k2, the maximal kinetic rate of the RyR.
Fig 7 shows four Ca2+ traces with kα = 0.5, 0.9, 1, and 1.25 when a = 0.25 and p = 10 are fixed. Notice that for kα = 1 aberrant Ca2+ signals emerge. Also included in Fig 7A is the corresponding partial bifurcation diagram for a = 0.25 using kα as the bifurcation parameter. Notice that MMOs exist for a large parameter set of kα. In Fig 7A there are four period doubling bifurcations that each have been labeled as PD. These occur around the parameter values of kα = {0.5893, 0.7599, 1.01, 1.026}, respectively. The single Hopf bifurcation occurs around the parameter value kα = 1.313. The bifurcation diagram in A provides us with a way to predict the Ca2+ signals for a large set of the parameter kα. Notice that stable periodic solutions occur for the parameter intervals kαs = (0, 0.5893) ∪ (1.026, 1.313). We have MMOs for the parameter interval kαmmo = (0.5893, 1.026) with various sub-threshold oscillations patterns. Fig 7D shows aberrant Ca2+ signals as the parameter kα remains close to the period doubling value 1.01. Fig 7E shows small stable periodic oscillations for the value kα = 1.25 near the Hopf point.
This figure shows a bifurcation diagram and four Ca2+ traces. A shows a partial bifurcation diagram for the parameter kα. In this figure, a number of period doubling points have been labeled as PD. The shaded region corresponds to MMOs. The single Hopf bifurcation point has been labeled with HB. B-E show four solutions for the parameter values kα = 0.5, 0.9, 1, and 1.25, respectively. D shows aberrant Ca2+ signals.
Each trace illustrated in Fig 7 shows a different type of Ca2+ response based on the value of kα. One possible way of dealing with the aberrant signals and MMOs is to drive the maximal reaction rate of the RyR up. This corresponds to increasing the value of k2 in the model. Fig 8A shows a two-parameter bifurcation diagram with kα on the x-axis and k2 on the y-axis. What this diagram illustrates is the region of MMOs produced by the model for the parameter ranges of kα and k2 given in the figure. Suppose that for a fixed kα value, Ca2+ signals undergo aberrant signals or MMOs. By increasing the parameter k2 sufficiently, passage into stable periodic solutions will occur. This suggests that the maximal rate of the RyR kinetics may help to control both aberrant Ca2+ signals and MMOs. To see this, Fig 8B shows the solution corresponding to those in Fig 7C with kα = 0.9 when k2 = 0.18 is increased to k2 = 0.5. Similarly, Fig 8C shows the solution corresponding to those in Fig 7D with kα = 1 when k2 = 0.18 is increased to k2 = 0.65. Notice that although solutions do settle into stable periodic orbits that the amplitude of the signals increase. Thus, increasing the parameter k2 may help to stabilizes aberrant oscillations at the cost of increasing Ca2+ oscillation amplitude. A similar stable region also exists below the region of MMOs. This region would also help control signals but may be more difficult to precisely isolate the appropriate range of k2.
This figures shows the effects of changing the parameter k2 on model solutions. A shows a two parameter bifurcation diagram with kα and k2 as the bifurcation parameters. The shaded region corresponds to the region of MMOs. The solid curve corresponds to the manifold of Hopf bifurcation points. The two dots on the left side represents the locations of the parameter values of k2 used to generate the Ca2+ traces in Fig 7C and Fig 8B. The two dots on the right side represents the locations of the parameter values of k2 used to generate the Ca2+ traces in Fig 7D and Fig 8C. B and C show the stable periodic oscillations that occur when k2 is increased to k2 = 0.5 and k2 = 0.65, respectively.
Dynamic levels of IP3 with no membrane potential.
The model dynamics exhibited in the figures above are triggered by a constant value of IP3 in the cytosol. Recall that in vivo, IP3 production is typically a result of an agonist activated G-couple protein. The rates of IP3 production and degradation are both modulated by intracellular Ca2+, and as such, the release of Ca2+ through the IP3 receptor can directly alter the IP3 signaling pathway. Furthermore, there is growing evidence that Aβ also affects the IP3 signaling pathway [7]. Thus, we extend our model to include a dynamic variable for IP3 production and degradation, and look to include the influence of Aβ on this signaling mechanism when ps = 0. More specifically, we include an additional equation for IP3 and model the influence of Aβ on the production of IP3.
We make use of the hybrid model formulated by Politi et al. (2006) to track IP3 production and degradation. Their model takes the form
(20)
where
The first term on the right hand side of (20) corresponds to the production of IP3 and the second term is the degradation. In (20), VPLC is the maximal production rate, KPLC characterizes the sensitivity of PLC, K3K is the half saturation constant for the degradation term. The constants k3K and k5P correspond to the IP3 phosphorylation and dephosphorylation rates, respectively. KPLC and η are parameters used to adjust the positive and negative feedback of Ca2+, respectively. One of the advantages for using this hybrid model is that it can easily be altered to reproduce both class I and class II mechanisms (see [18] for details). Another advantage, is that this model breaks the IP3 process into two components: a production and a degradation. This will make it easier for us to incorporate the effects of Aβ on the IP3 production process.
In their experiment, Demuro and Parker (2013) showed that introducing Aβ directly into Xenopus oocytes causes an increase in Ca2+ dependent fluorescence (a measure for the amount of intracellular Ca2+). Even though their experiments are in oocytes, the ubiquitous properties of IP3 signaling may make their results relevant to other cells including neurons. Their findings suggest that Aβ does not interact directly with the IP3 receptor, but instead they propose that intracellular Ca2+ liberation evoked by Aβ involves opening of IP3 receptors as a result of stimulated production of IP3 via G-protein-mediated activation of PLC. As such, in the presence of Aβ, IP3 are actively stimulated and persist for many minutes or hours even though IP3 is metabolized within tens of seconds [58]. Based on these findings, we assume that VPLC takes the form
(21)
and KPLC can be written similarly as
(22)
where the parameters μPLC and κPLC control the strength of the linear influence of Aβ on each term, respectively. For our purposes, we set both of these values to be μPLC = κPLC = 1. Thus, we add the following equation to (18) and (19) and look to determine the impact of Aβ on model solutions
(23)
In our simulation we use the work of [59] to set a number of parameter values. However, we assume that both positive and negative feedback are present simultaneously and as such make parameter adjustments as needed. The parameter values that we use for our simulations are given in Table 2. With these additional contributions, we now have a model that includes the impact of Aβ on multiple Ca2+ signaling mechanisms. Although the model has a large number of variables and parameters, we seek to characterize model solutions by investigating the dynamical properties of the model in the presence of Aβ.
In the absence of Aβ, model solution with initial condition for IP3 = 0.01 shows a small Ca2+ influx followed by a transition to it’s steady-state level close to 0.5 μM. This is consistent with what we expect as the amount of IP3 is dynamic and depends on Ca2+. It will take sometime for enough IP3 to be present in order to trigger a signaling event throught the receptor. Fig 9A and 9B show the model solution for Ca2+ and IP3, respectively in the absence of Aβ. Notice that the amount of IP3 also reaches a steady-state level.
A shows Ca2+ and B shows IP3 as a function of time in the absence of Aβ. A spike of Ca2+ occurs once enough IP3 has accumulated but does not lead to sustained oscillations. The amount of IP3 settles to a steady-state level in B.
Since we are interested in understanding the potential effect of Aβ on Ca2+ signals, we use a as a bifurcation parameter and investigate model solutions. Fig 10A shows a partial bifurcation diagram with a as the bifurcation parameter. This diagram has two oscillatory regions separated by a region of a single stable-steady state. The four Hopf bifurcations are labeled along with a number of period doubling points. Notice that there are four regions of MMOs (shaded regions) that appear close to each Hopf point. In addition to the Hopf bifurcation points, four period doubling points have also been labeled with two limit points (LP). We provide the value of each of these points in Table 3 and use them to identify solution patterns. As the parameter a increases to LP3 = 1.274, solutions become unphysical. Thus, we limit our investigation for values of a between 0 and 1.274. Fig 10B–10D show Ca2+ traces predicted by Fig 10A for various values of a. Notice that aberrant signals are also present for the model with a dynamic equation for IP3.
This figure shows a bifurcation diagram and four Ca2+ traces for the model with dynamic IP3. A shows a partial bifurcation diagram where two oscillatory regions are separated by a single steady-state region. The shaded regions near each of the four Hopf bifurcations correspond to MMOs. Numerous period doubling bifurcations occur around the shaded regions. Two important limit points have been labeled LP1 and LP1. These points are important in the description of solutions. B-D show model solutions for a = 0.75, 0.875, 1, and 1.256, respectively. The various patterns in these figures are predicted by the bifurcation diagram in A.
Calcium model with influence of membrane potential
In our simulations below, we examine the effects of changes to membrane potentials by first incorporating the membrane into the model and by stimulating the membrane with a constant applied current pulse. The impact of such an applied current on the voltage V and the resulting Ca2+ current are shown in Fig 11. A stimulating current of 300 nA applied for t = 0.1 seconds initiated at t = 2 generates the potential response shown in Fig 11A. Fig 11B and 11C show the response due to an applied stimulus lasting for one and ten seconds, respectively. These changes in membrane potential cause the inward currents through the VGCC illustrated in Fig 11D–11F, for t = 0.1, t = 1, and t = 10 seconds, respectively. Notice that in Fig 11C the membrane voltage saturates as the duration of the applied current is increased. Although the mechanisms for generating these signals is fairly simplistic, our results align well with experimental data showing similar saturating levels in astrocytes [38]. The effects of the inward Ca2+ through the VGCC are included in the flux term Jvca and as such, allows us to study the impact of membrane potentials on Ca2+ signals. In all subsequent figures, we have applied a constant current pulse at t = 100 for a duration of 50 seconds. Such an applied current does not capture an in vivo-like representation of membrane potentials, but it does offer a way to link the model with typical voltage clamp experiments where the amplitude and duration of the applied current can be controlled.
This figures shows the Ca2+ current Ica in response to changes in membrane potential. A-C show changes in the membrane potential of the astrocytic model when a current of 300 nA is applied for a duration of 0.1 second (A), 1 second (B), and 10 seconds (C). The current was applied at t = 2 seconds and is represented by the black bar at the bottom of each figure. D-F show the corresponding VGCC current Ica as defined by (15).
To investigate the impact of a constant current pulse on Ca2+ signaling, we first simulate (18) and (19) with constant values of IP3. In these simulations we set ps = 1 and plot Ca2+ concentrations as a function of time. Tracking the membrane potential and including it into the model will have an effect on Ca2+ signals. Specifically, the membrane will alter the dynamics of the effect of p on model solutions. In order to illustrate this point we have plotted the model responses when no Aβ is present for various values of p in Fig 12. When these solutions are compared to Fig 2, we can see that the inclusion of the membrane potential increases the response frequency and MMOs occur for smaller values of p (for example p = 13 here instead of p = 17 in Fig 2).
This figures shows the impact of membrane potential on Ca2+ signaling. In each figure, membrane potentials are included as a response to a sustained applied current of 300 nA lasting from t = 100 to t = 150. Figs A, B, and C show intracellular Ca2+ signals when p = 5, p = 10, and p = 12 with a = 0, respectively. Figs D, E, and F show the response when p = 13, p = 13.5, and p = 14, respectively. Note that the inclusion of the membrane potential filters MMOs and can help establish transient single-mode oscillations upon the termination of the applied current (D,E).
Notice that in Fig 12D–12F Ca2+ signals seem to stabilize in single-mode oscillations upon the release of the stimulus before transitioning back into MMOs. This may provide one possible path for stabilizing Ca2+ signaling. In an in-vivo like environment, membrane potentials vary based on intrinsic response mechanisms to a stimulus. Our results show that artificially triggering membrane stimulation could potentially help to stabilize Ca2+ signals. Furthermore, Fig 12E shows that the model solution can exhibit a large number of sub-threshold oscillations before triggering a larger spike. This type of MMO is different from those illustrated previously and shows that membrane potentials can have an impact on global Ca2+ signals and play an important role in Ca2+ regulation. A full analysis for understanding these transitions is beyond the scope of this work but could prove useful for predicting how membrane stimulation may be used to control and/or stabilize aberrant Ca2+ signals.
Model solution patterns are not only linked to the amount of IP3 and the inclusion of membrane potential, but also by the amount of Aβ present in the model. To better understand the roles of IP3, membrane potentials, and Aβ on Ca2+ signaling, we simulate the model and provide the solutions when p = 10 and p = 15 with various levels of Aβ for ps = 1. Illustrated in Fig 13 are six model solutions that illustrate the impact of Aβ and membrane potentials on Ca2+ signaling. When p = 10 is fixed and a is altered, model behavior is directly impacted by the application and removal of a constant current pulse. This is illustrated in Fig 13A–13C where a = 0.2, a = 0.25, and a = 0.28. Specifically, Fig 13C shows that upon the termination of the applied stimulus, Ca2+ signals do not go back to MMO patterns but instead transition to stable single-mode oscillations. This shows that stimulation of the membrane can alter intrinsic dynamical patterns and be used to stabilize various types of Ca2+ signals. Fig 13D–13F show model solutions for p = 15 with a = 0.2, a = 0.25, and a = 0.28, respectively. It is interesting to note that when we apply a constant current pulse, Ca2+ solutions can transition into steady MMOs as the level of Aβ increased towards a = 0.28. Although we do not analyze the bifurcation structure, the inclusion of membrane potentials in the model appears to have altered the parameter dependance where regions of MMOs can occur as well as transitions from single- to mixed-mode oscillations. Further analysis may be beneficial for drawing out the underlying mechanisms in the stable oscillatory patterns when a constant current is applied in the model.
This figures shows the impact of membrane potential on Ca2+ signaling. In each figure, changes in membrane potential are included as a response to a sustained applied current of 300 lasting from t = 100 to t = 150. Figs A, B, and C show intracellular Ca2+ signals when p = 10 and when a = 0.2, a = 0.25, and a = 0.28, respectively. Figs D, E, and F show the response when p = 15 and when a = 0.2, a = 0.25, and a = 0.28, respectively.
Because the impact of membrane potentials are controlled by ps, we have also provided model simulations when this parameter is altered. Small values of ps correspond to small influence of membrane potential while larger values can be used drive the amplitude of Ca2+ signals. Fig 14 shows four simulations for four different values of ps. When ps = 0.5, Fig 14A shows that the effects of membrane potentials are small and no significant changes in Ca2+ are observed other than a slight increase in the subtreshold oscillations. When ps = 0.75, Fig 14B suggests that the effects of membrane potentials are large enough to alter the amplitude of Ca2+ oscillations during the applied current. Fig 14C shows that when p = 1, contributions from membrane potentials increase overall Ca2+ signaling amplitude. With these values of p and a, the solution transitions from single mode oscillations with amplitude around 0.45 to MMO with an increased amplitude to around 0.7. Fig 14D shows that a large influence of membrane potentials (when p = 1.5) does not alter the signal amplitude or oscillatory mode significantly.
This figures shows various model solutions when the scaling parameter ps is altered. In each figure, p = 10, and the amount of Aβ is fixed at a = 0.258. A-D show the response of (18) and (19) when ps = 0.5, ps = 0.75, ps = 1, and ps = 1.5, respectively. Notice that B shows that Ca2+ can enter aberrant oscillatory patterns when the membrane is stimulated by a constant applied current. C shows that Ca2+ signals can enter MMOs with altered amplitudes when the applied current stimulus is turned off.
Although much analysis remains to fully understand the dynamics of the model, the results of our simulations suggest that Aβ can alter Ca2+ regulatory mechanisms in a way that leads to both MMOs and aberrant signaling. We have shown that bifurcation regions for dynamical transitions from single mode to MMOs can increase in the presence of Aβ. By altering RyR receptor dynamics, we show that transitions from MMO back to single-mode oscillations can occur. Furthermore, we show that stimulation of the membrane can also be used to control various types of Ca2+ signals. Although we have made a number of simplifying assumptions in the model development, our approach can be easily altered to include other more complex interactions and mechanisms that influence Ca2+ regulation.
Discussion
Intracellular Ca2+ is a critically important second messenger within the nervous system. In neurons, Ca2+ is known to mediate the signaling pathways that control neurotransmitter release, gene expression, metabolism, plasticity, development, proliferation, and cell death [60]. As such, Ca2+ may play a major role in the pathogenesis of AD. Unfortunately, the complexity of Ca2+ signaling makes it difficult to precisely understand how Aβ impacts different intracellular regulation mechanisms and components. Various studies have decoupled particular components and merging these theories together to form a whole-cell computational model can help us better understand intracellular Ca2+ regulation and what leads to aberrant signaling. One of our goals is to model Ca2+ in a simplified whole-cell environment that has predictable qualitative structure so that we can study the effects of Aβ on the dynamics. As such, the qualitative features described in this study give us a way to track Ca2+ patterns using dynamical systems theory.
The simulations presented in this study occur on the order of seconds to minutes while the progression of AD occurs on the timescale of months to years. However, in our model we are using the accumulation of Aβ to describe the potential stage in the evolution of AD. Even before toxic Aβ plaques can aggregate, the slow accumulation of Aβ peptides can trigger alterations in Ca2+ signaling patterns. Our model shows that aberrant signals and changes in homeostasis levels can emerge as the amount of Aβ is increased. In an in vivo environment these changes may be subtle and actually evolve over days or months. Any alterations in intracellular Ca2+ homeostasis can affect the apoptotic signaling cascade. Both the mitochondria and the ER play a significant role in apoptosis and are sensitive to changes in Ca2+ levels. Although we have not considered mitochondrial effects, we do track ER Ca2+. Further analysis that looks at the time evolution of ce could be useful in predicting chronic changes of Ca2+ homeostasis in AD.
Developing a whole-cell Ca2+ model that has predictable qualitative structure in the presence of Aβ is challenging. Although Aβ influences many Ca2+ regulatory mechanisms, the particular way which Aβ affects these mechanisms is generally not known. Additionally, the temporal influence of Aβ on certain mechanisms could occur on the order of milliseconds, seconds, days, months, or years. As such, any computational model will necessarily make a number of simplifying assumptions. Even by exploiting these simplifications, our model includes a large number of parameters that make mathematical analysis limited. Unfortunately, we do not have robust estimates for many of the parameters involved in the model. However, we have attempted to provide justification for many assumptions and parameter choices based on the literature and the experimental data currently available. We do recognize that many of these assumptions may need to be altered as we continue to improve our understanding of the effects of Aβ in an AD environment.
The ubiquitous nature of the Ca2+ regulatory mechanisms used in our model makes it easily adaptable for studying various cell types with spatial components. Specifically, Aβ has been shown to cause complex Ca2+ signals in astrocytes [61, 62]. In these astrocytes, Ca2+ waves and oscillations signals can occur on timescales even slower than those typical of other non-excitable neuroglia. Although our modeling approach does not include spatial components, additional mechanisms can be constructed to account for wave generating behaviors. Furthermore, astrocytes can facilitate synaptic transmission and plasticity through the uptake of neurotransmitter [63] and complex models between neurons and astrocytes have been developed to study these interactions [38, 64, 65]. Both microglia and astrocytes have been described as modulators for Aβ clearance and degradation [66] and our approach may be useful for better understanding these mechanisms.
It is clear that Aβ plays an essential role in the cognitive decline in AD by directly affecting synaptic transmission [11, 67–70]. However, synaptic transmission is typically precipitated by a presynaptic potential which allows Ca2+ ions to flow into the cell through VGCC. The contributions of fast local Ca2+ signals with slow global Ca2+ patterns, especially under the influence of Aβ, may help explain why breakdowns in synaptic efficacy can occur in an AD environment [11, 45, 67, 69, 71, 72]. The accumulation or presence of Aβ may directly, or indirectly, impact various Ca2+ driven mechanisms during synaptic transmission. The simplified whole-cell Ca2+ model presented here could be linked with a synaptic Ca2+ model to investigate how global aberrant Ca2+ signals may impact synaptic transmission on multiple timescales. Simulations over long timescales may help explain how slow global whole-cell Ca2+ signals interfere with fast local Ca2+ signals at the synapse.
Different individual models exist for the various signaling components used in our simplified whole-cell model development. For example, we used the Sneyd and Dufour (2002) formulation for the IP3 receptor model. Although this model is sound and well-suited for our purposes, it does increase the number of necessary variables considerably. One could use a two equation model for IP3 (such as those described in [26, 73]), but the number of parameters will remain large. Similarly, alternative models for the RyR may tease out alternative conclusions when influenced by Aβ (such as using a model as in [74]). As such, we encourage further development of the model as experimental data becomes available. Matching model dynamics with experimentally recorded data can help select the component model best suited for the particular study.
Our model solutions are highly sensitive to certain parameters and the oscillatory responses presented here only occur under certain scenarios. Because of the complexity of Ca2+ regulation along with understanding the impact of Aβ, any computational model would benefit from both local and global sensitivity analysis. Although we have not performed any sensitivity analysis, we do understand that much of the analysis and many of our conclusions may be valid for a small set of parameters. Further, the sensitivity of a particular parameter may influence how the model transitions into aberrant Ca2+ signals. As such, we recommend that sensitivity analysis be performed as a next step in order to better understand the role of parameters on model dynamics.
In conclusion, we have shown that aberrant Ca2+ signals can occur in a simplified whole-cell model under the influence of Aβ. Furthermore, we showed that regions of MMOs can expand as a consequence of increasing the amount of Aβ in the model. This may partially explain how Ca2+ signals are impacted by Aβ from a dynamics perspective within an in vivo like environment. Continued refinement of the model in conjunction with experimental data matching will help make the model more useful. In turn, this can help us determine how to control for both aberrant signals and increased homeostasis Ca2+ levels. The model can then be used to better understand the impact of Aβ on Ca2+ fluxes through individual regulatory components (such as IP3, RyR, and plasma membrane). This computational model can help us study complex cellular behavior in an AD environment by tracking the influence of many interconnected biological mechanisms.
Supporting information
S1 Appendix. A summary of both the IP3 receptor model equations and those governing the membrane potential.
Also provided in the appendix are the model parameters used in the Hodgkin and Huxley formulation of the membrane potential.
https://doi.org/10.1371/journal.pone.0202503.s001
(PDF)
Acknowledgments
We are indebted to the anonymous referees for providing suggestions that helped us improve this manuscript. We are also grateful to Norwich University for their support of this work.
References
- 1. Alzheimer’s Association. 2016 Alzheimer’s disease facts and figures. Alzheimers Dement. 2016;12(4):459–509. pmid:27570871
- 2. Gouras GK, Tsai J, Naslund J, Vincent B, Edgar M, Checler F, et al. Intraneuronal Aβ42 Accumulation in Human Brain. The American Journal of Pathology. 2000;156(1):15–20. pmid:10623648
- 3. Oddo S, Caccamo A, Smith IF, Green KN, LaFerla FM. A dynamic relationship between intracellular and extracellular pools of Aβ. Am J Pathol. 2006;168(1):184–94. pmid:16400022
- 4. Ferreiro E, Oliveira CR, Pereira C. Involvement of endoplasmic reticulum Ca2+ release through ryanodine and inositol 1,4,5-triphosphate receptors in the neurotoxic effects induced by the amyloid-β peptide. J Neurosci Res. 2004;76(6):872–80. pmid:15160398
- 5. Kuchibhotla KV, Goldman ST, Lattarulo CR, Wu HY, Hyman BT, Bacskai BJ. Aβ plaques lead to aberrant regulation of calcium homeostasis in vivo resulting in structural and functional disruption of neuronal networks. Neuron. 2008;59(2):214–25. pmid:18667150
- 6. Demuro A, Mina E, Kayed R, Milton SC, Parker I, Glabe CG. Calcium dysregulation and membrane disruption as a ubiquitous neurotoxic mechanism of soluble amyloid oligomers. J Biol Chem. 2005;280(17):17294–300. pmid:15722360
- 7. Demuro A, Parker I. Cytotoxicity of intracellular Aβ 42 amyloid oligomers involves Ca2+ release from the endoplasmic reticulum by stimulated production of inositol trisphosphate. J Neurosci. 2013;33(9):3824–33. pmid:23447594
- 8. Lopez JR, Lyckman A, Oddo S, Laferla FM, Querfurth HW, Shtifman A. Increased intraneuronal resting [Ca2+] in adult Alzheimer’s disease mice. J Neurochem. 2008;105(1):262–71. pmid:18021291
- 9. Hardy J, Allsop D. Amyloid deposition as the central event in the aetiology of Alzheimer’s disease. Trends Pharmacol Sci. 1991;12(10):383–8. pmid:1763432
- 10. Hardy J, Selkoe DJ. The amyloid hypothesis of Alzheimer’s disease: progress and problems on the road to therapeutics. Science. 2002;297(5580):353–6. pmid:12130773
- 11. Selkoe DJ, Hardy J. The amyloid hypothesis of Alzheimer’s disease at 25 years. EMBO Mol Med. 2016;8(6):595–608. pmid:27025652
- 12. Selkoe DJ. The molecular pathology of Alzheimer’s disease. Neuron. 1991;6(4):487–98. pmid:1673054
- 13. Paula-Lima AC, Adasme T, SanMartin C, Sebollela A, Hetz C, Carrasco MA, et al. Amyloid β-peptide oligomers stimulate RyR-mediated Ca2+ release inducing mitochondrial fragmentation in hippocampal neurons and prevent RyR-mediated dendritic spine remodeling produced by BDNF. Antioxid Redox Signal. 2011;14(7):1209–23. pmid:20712397
- 14. Shtifman A, Ward CW, Laver DR, Bannister ML, Lopez JR, Kitazawa M, et al. Amyloid-β protein impairs Ca2+ release and contractility in skeletal muscle. Neurobiol Aging. 2010;31(12):2080–90. pmid:19108934
- 15. Demuro A, Smith M, Parker I. Single-channel Ca2+ imaging implicates Aβ1-42 amyloid pores in Alzheimer’s disease pathology. J Cell Biol. 2011;195(3):515–24. pmid:22024165
- 16. Ullah G, Demuro A, Parker I, Pearson JE. Analyzing and Modeling the Kinetics of Amyloid Beta Pores Associated with Alzheimer’s Disease Pathology. PLoS One. 2015;10(9):e0137357. pmid:26348728
- 17. Berridge MJ, Bootman MD, Roderick HL. Calcium signalling: dynamics, homeostasis and remodelling. Nat Rev Mol Cell Biol. 2003;4(7):517–29. pmid:12838335
- 18.
Dupont G, Falcke M, Kirk V, Sneyd J. Models of Calcium Signaling. Springer International Publishing Switzerland; 2016.
- 19.
Fall CP, Marland ES, Wagner JM, Tyson JJ. Computational Cell Biology. New York: Springer Science+Business Media, Inc.; 2002.
- 20.
Keener J, Sneyd J. Mathematical Physiology I: Cellular Physiology. 2nd ed. Springer Science+Business Media, LLC; 2009.
- 21. Manninen T, Havela R, Linne ML. Computational Models for Calcium-Mediated Astrocyte Functions. Front Comput Neurosci. 2018;12:14. pmid:29670517
- 22. Berridge MJ. The Inositol Trisphosphate/Calcium Signaling Pathway in Health and Disease. Physiol Rev. 2016;96(4):1261–96. pmid:27512009
- 23. Berridge MJ. Inositol trisphosphate and calcium signalling. Nature. 1993;361(6410):315–25. pmid:8381210
- 24. Sneyd J, Dufour JF. A dynamic model of the type-2 inositol trisphosphate receptor. Proc Natl Acad Sci U S A. 2002;99(4):2398–403. pmid:11842185
- 25. Sneyd J, Tsaneva-Atanasova K, Bruce JI, Straub SV, Giovannucci DR, Yule DI. A model of calcium waves in pancreatic and parotid acinar cells. Biophys J. 2003;85(3):1392–405. pmid:12944257
- 26. De Young GW, Keizer J. A single-pool inositol 1,4,5-trisphosphate-receptor-based model for agonist-stimulated oscillations in Ca2+ concentration. Proceedings of the National Academy of Sciences. 1992;89(20):9895–9899.
- 27.
Othmer H, Tang Y. Oscillations and waves in a model of InsP3-controlled calcium dynamics. In: Othmer H, Murry J, Maini P, editors. Experimental and Theoretical Advances in Biological Pattern Formation. vol. 259. London: Plenum Press; 1993. p. 277–300.
- 28. LeBeau AP, Yule DI, Groblewski GE, Sneyd J. Agonist-dependent phosphorylation of the inositol 1,4,5-trisphosphate receptor: A possible mechanism for agonist-specific calcium oscillations in pancreatic acinar cells. J Gen Physiol. 1999;113(6):851–72 pmid:10352035
- 29. Friel DD. [Ca2+]i oscillations in sympathetic neurons: an experimental test of a theoretical model. Biophys J. 1995;68(5):1752–66. pmid:7612818
- 30. Dupont G, Goldbeter A. One-pool model for Ca2+ oscillations involving Ca2+ and inositol 1,4,5-trisphosphate as co-agonists for Ca2+ release. Cell Calcium. 1993;14(4):311–22. pmid:8370067
- 31. Lytton J, Westlin M, Burk SE, Shull GE, MacLennan DH. Functional comparisons between isoforms of the sarcoplasmic or endoplasmic reticulum family of calcium pumps. J Biol Chem. 1992;267(20):14483–9. pmid:1385815
- 32. Favre CJ, Schrenzel J, Jacquet J, Lew DP, Krause KH. Highly supralinear feedback inhibition of Ca2+ uptake by the Ca2+ load of intracellular stores. J Biol Chem. 1996;271(25):14925–30. pmid:8662967
- 33. Mogami H, Tepikin AV, Petersen OH. Termination of cytosolic Ca2+ signals: Ca2+ reuptake into intracellular stores is regulated by the free Ca2+ concentration in the store lumen. EMBO J. 1998;17(2):435–42. pmid:9430635
- 34. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117(4):500–44. pmid:12991237
- 35. Seifert G, Henneberger C, Steinhäuser C. Diversity of astrocyte potassium channels: An update. Brain Res Bull. 2018;136:26–36. pmid:27965079
- 36. Breslin K, Wade JJ, Wong-Lin K, Harkin J, Flanagan B, Van Zalinge H, et al. Potassium and sodium microdomains in thin astroglial processes: A computational model study. PLoS Comput Biol. 2018;14(5):e1006151. pmid:29775457
- 37. Dallérac G, Chever O, Rouach N. How do astrocytes shape synaptic transmission? Insights from electrophysiology. Front Cell Neurosci. 2013;7:159. pmid:24101894
- 38. Sibille J, Dao Duc K, Holcman D, Rouach N. The neuroglial potassium cycle during neurotransmission: role of Kir4.1 channels. PLoS Comput Biol. 2015;11(3):e1004137. pmid:25826753
- 39. LeBeau AP, Van Goor F, Stojilkovic SS, Sherman A. Modeling of membrane excitability in gonadotropin-releasing hormone-secreting hypothalamic neurons regulated by Ca2+-mobilizing and adenylyl cyclase-coupled receptors. J Neurosci. 2000;20(24):9290–7. pmid:11125008
- 40. Demuro A, Parker I, Stutzmann GE. Calcium signaling and amyloid toxicity in Alzheimer disease. J Biol Chem. 2010;285(17):12463–8. pmid:20212036
- 41. Itkin A, Dupres V, Dufrêne YF, Bechinger B, Ruysschaert JM, Raussens V. Calcium ions promote formation of amyloid β-peptide (1-40) oligomers causally implicated in neuronal toxicity of Alzheimer’s disease. PLoS One. 2011;6(3):e18250. pmid:21464905
- 42. De Caluwé J, Dupont G. The progression towards Alzheimer’s disease described as a bistable switch arising from the positive loop between amyloids and Ca(2+). J Theor Biol. 2013;331:12–8. pmid:23614875
- 43. Briggs CA, Chakroborty S, Stutzmann GE. Emerging pathways driving early synaptic pathology in Alzheimer’s disease. Biochem Biophys Res Commun. 2017;483(4):988–997. pmid:27659710
- 44. Del Prete D, Checler F, Chami M. Ryanodine receptors: physiological function and deregulation in Alzheimer disease. Mol Neurodegener. 2014;9:21. pmid:24902695
- 45. Liang J, Kulasiri D, Samarasinghe S. Ca2+ dysregulation in the endoplasmic reticulum related to Alzheimer’s disease: A review on experimental progress and computational modeling. Biosystems. 2015;134:1–15. pmid:25998697
- 46. Briggs CA, Schneider C, Richardson JC, Stutzmann GE. beta amyloid peptide plaques fail to alter evoked neuronal calcium signals in APP/PS1 Alzheimer’s disease mice. Neurobiol Aging. 2013;34(6):1632–43. pmid:23337342
- 47. Chakroborty S, Stutzmann GE. Calcium channelopathies and Alzheimer’s disease: Insight into therapeutic success and failures. Eur J Pharmacol. 2014;739:83–95. pmid:24316360
- 48. Goussakov I, Miller MB, Stutzmann GE. NMDA-mediated Ca(2+) influx drives aberrant ryanodine receptor activation in dendrites of young Alzheimer’s disease mice. J Neurosci. 2010;30(36):12128–37. pmid:20826675
- 49. Stutzmann GE, Smith I, Caccamo A, Oddo S, Laferla FM, Parker I. Enhanced ryanodine receptor recruitment contributes to Ca2+ disruptions in young, adult, and aged Alzheimer’s disease mice. J Neurosci. 2006;26(19):5180–9. pmid:16687509
- 50. Costa RO, Lacor PN, Ferreira IL, Resende R, Auberson YP, Klein WL, et al. Endoplasmic reticulum stress occurs downstream of GluN2B subunit of N-methyl-d-aspartate receptor in mature hippocampal cultures treated with amyloid-beta oligomers. Aging Cell. 2012;11(5):823–33. pmid:22708890
- 51. Ferreira IL, Bajouco LM, Mota SI, Auberson YP, Oliveira CR, Rego AC. Amyloid beta peptide 1-42 disturbs intracellular calcium homeostasis through activation of GluN2B-containing N-methyl-d-aspartate receptors in cortical cultures. Cell Calcium. 2012;51(2):95–106. pmid:22177709
- 52.
MATLAB 2017b, The MathWorks, Inc., Natick, MA, USA; 2017.
- 53.
Doedel EJ. AUTO: A program for the automatic bifurcation and analysis of autonomous systems. In: Proc. 10th Manitoba Conf. Num. Anal. and Comp.; 1981. p. 265–284.
- 54.
Ermentrout B. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. Philadelphia: SIAM; 2002.
- 55. Brøns M, Desroches M, Krupa M. Mixed-Mode Oscillations Due to a Singular Hopf Bifurcation in a Forest Pest Model. Mathematical Population Studies. 2015;22(2):71–79.
- 56. Desroches M, Guckenheimer J, Krauskopf B, Kuehn C, Osinga HM, Wechselberger M. Mixed-Mode Oscillations with Multiple Time Scales. SIAM Review. 2012;54(2):211–288.
- 57. Harvey E, Kirk V, Wechselberger M, Sneyd J. Multiple Timescales, Mixed Mode Oscillations and Canards in Models of Intracellular Calcium Dynamics. Journal of Nonlinear Science. 2011;21(5):639–683.
- 58. Sims CE, Allbritton NL. Metabolism of inositol 1,4,5-trisphosphate and inositol 1,3,4,5-tetrakisphosphate by the oocytes of Xenopus laevis. J Biol Chem. 1998;273(7):4052–8. pmid:9461597
- 59. Politi A, Gaspers LD, Thomas AP, Hofer T. Models of IP3 and Ca2+ oscillations: frequency encoding and identification of underlying feedbacks. Biophys J. 2006;90(9):3120–33. pmid:16500959
- 60. Supnet C, Bezprozvanny I. The dysregulation of intracellular calcium in Alzheimer disease. Cell Calcium. 2010;47(2):183–9. pmid:20080301
- 61. Abramov AY, Canevari L, Duchen MR. Changes in intracellular calcium and glutathione in astrocytes as the primary mechanism of amyloid neurotoxicity. J Neurosci. 2003;23(12):5088–95. pmid:12832532
- 62. Abramov AY, Canevari L, Duchen MR. Calcium signals induced by amyloid beta peptide and their consequences in neurons and astrocytes in culture. Biochim Biophys Acta. 2004;1742(1-3):81–7. pmid:15590058
- 63. Pannasch U, Vargová L, Reingruber J, Ezan P, Holcman D, Giaume C, et al. Astroglial networks scale synaptic activity and plasticity. Proc Natl Acad Sci U S A. 2011;108(20):8467–72. pmid:21536893
- 64. Kenny A, Plank MJ, David T. The role of astrocytic calcium and TRPV4 channels in neurovascular coupling. J Comput Neurosci. 2018;44(1):97–114. pmid:29152668
- 65. Tang J, Liu T, Luo J, Yang X. Effect of calcium channel noise in astrocytes on neuronal transmission. Communications in Nonlinear Science and Numerical Simulation. 2016;32:262–272.
- 66. Ries M, Sastre M. Mechanisms of Aβ Clearance and Degradation by Glial Cells. Frontiers in Aging Neuroscience. 2016;8(160). pmid:27458370
- 67. Palop JJ, Mucke L. Amyloid-beta-induced neuronal dysfunction in Alzheimer’s disease: from synapses toward neural networks. Nat Neurosci. 2010;13(7):812–8. pmid:20581818
- 68. Ripoli C, Cocco S, Li Puma DD, Piacentini R, Mastrodonato A, Scala F, et al. Intracellular accumulation of amyloid-β (Aβ) protein plays a major role in Aβ-induced alterations of glutamatergic synaptic transmission and plasticity. J Neurosci. 2014;34(38):12893–903. pmid:25232124
- 69. Ripoli C, Piacentini R, Riccardi E, Leone L, Li Puma DD, Bitan G, et al. Effects of different amyloid β-protein analogues on synaptic function. Neurobiol Aging. 2013;34(4):1032–44. pmid:23046860
- 70. Selkoe DJ. Alzheimer’s disease is a synaptic failure. Science. 2002;298(5594):789–91. pmid:12399581
- 71.
Dehkordy SR, Bahrami F, Janahmadi M. Computational study of the role of calcium in late long-term potentiation induction on the basis of tripartite synapse structure. In: Proceeding from (ICEE) Electrical Engineering 22nd Iranian Conference; 2014.
- 72. Popugaeva E, Bezprozvanny I. Role of endoplasmic reticulum Ca2+ signaling in the pathogenesis of Alzheimer disease. Front Mol Neurosci. 2013;6:29. pmid:24065882
- 73. Li YX, Rinzel J. Equations for InsP3 receptor-mediated [Ca2+]i oscillations derived from a detailed kinetic model: a Hodgkin-Huxley like formalism. J Theor Biol. 1994;166(4):461–73. pmid:8176949
- 74. Keizer J, Levine L. Ryanodine receptor adaptation and Ca2+(-)induced Ca2+ release-dependent Ca2+ oscillations. Biophys J. 1996;71(6):3477–87. pmid:8968617