A mathematical model for the effects of amyloid beta on intracellular calcium

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.


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 Ca 2+ 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 Ca 2+ 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 Ca 2+ from internal stores such as the Endoplasmic Reticulum (ER) [4,[6][7][8]. As such, Aβ may lead to both local and global Ca 2+ proliferation that are persistent and cytotoxic. Sustained Ca 2+ 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 Ca 2+ 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 Ca 2+ signaling. The amyloid hypothesis suggests that the accumulation of Aβ in the brain is the primary driving force of AD pathogenesis [9][10][11][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 Ca 2+ signaling processes leading to synaptic failure and neuronal death. Although it is unclear how Aβ disrupts intracellular Ca 2+ homeostasis, there is growing evidence that Aβ directly affects the production of inositol triphosphate (IP 3 ) [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 Ca 2+ signaling mechanisms in a simplified whole-cell model.
In this study we present a theoretical approach to better understand the driving mechanisms for various Ca 2+ oscillatory patterns within an AD environment. By developing a mathematical model for intracellular Ca 2+ regulation, we can begin to study how Aβ affects Ca 2+ flux through various individual channels and pumps. Investigating model solutions can also provide important information on the impact of Aβ on Ca 2+ 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 Ca 2+ dynamics in AD, specifically through computational and mathematical modeling, and this work outlines such an approach.

Calcium model formulation
In order to study the effects of Aβ on intracellular Ca 2+ , we first build a simplified whole-cell Ca 2+ model by making use of the vast array of work on modeling Ca 2+ dynamics and the calcium signaling "toolkit" (see [17][18][19][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 Ca 2+ 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 Ca 2+ dynamics using traditional methods by tracking the flux in and out of the cytoplasm. Let c denote the concentration of free Ca 2+ ions in the cell cytoplasm, then the rate of change in intracellular Ca 2+ is governed by whereJ denotes flux. We assume a spatially homogeneous cell whose volume is fixed and track Ca 2+ 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 Ca 2+ signaling patterns. We take this position to study the influence of Aβ on Ca 2+ dynamics using a simplified whole-cell model. Many different mathematical approaches have been used to better understand Ca 2+ -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 Ca 2+ dynamics on the timescale of seconds. We assume that intracellular Ca 2+ in-fluxes (into the cytoplasm) are those corresponding to IP 3 receptors (IPR), RyRs, a general membrane leak J in , and we include a fast Voltage Gated Ca 2+ Channel (VGCC) in J vca to account for membrane permeability. The out-fluxes (out of the cell or sequestered into Ca 2+ pools such as the ER) are modeled using a plasma membrane pump J pm , and a sarco-endoplasmic reticulum Ca 2+ ATPase (SERCA) pump J SERCA . 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.
We assume that the ER is homogeneously distributed throughout the interior of the cell and that the fluxes of Ca 2+ through the IPR and RyR are proportional to the difference in the concentration of Ca 2+ between the ER and the cytoplasm. Under these standard assumptions, using conservation of fluxes, our general Ca 2+ model takes the form where c e denote the concentration of Ca 2+ 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 J vca in (2) links membrane potentials with other Ca 2+ signaling mechanisms. Ca 2+ 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 Ca 2+ 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 Ca 2+ 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.
Ca 2+ regulation in non-excitable cells, such as astrocytes, is extremely complex and influenced not only by membrane potential but also Ca 2+ buffering. Cytoplasmic Ca 2+ buffering plays a significant role in calcium's ability to move throughout the cell. Almost all of the available Ca 2+ is bound to buffers and free cytoplasmic Ca 2+ cannot move very far before being bound [18]). Because Ca 2+ diffusion plays an important role in spatiotemporal signaling, nonlinear buffering may provide some insights behind certain types of oscillatory patterns and Ca 2+ waves. However, compared to the temporal release and uptake of Ca 2+ into internal stores, buffering occurs on a much faster timescale. In our model we have assumed that Ca 2+ 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]). IP 3 receptor model. IP 3 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, IP 3 regulates many important cellular functions, including the release of Ca 2+ from the ER through the IP 3 receptor [23]. We assume that the flux from the IP 3 receptor follows a saturating binding rate model of the form found in [24,25]. Thus, we write where k f controls the density of the IP 3 receptor, J er is a leak from the ER to the cytoplasm, and P 0 is the open probability of the IPR. The leak term is necessary to balance the ATPase flux at steady state. To model P 0 , we use a version of the Sneyd and Dufour (2002) model that is based on previous models found in [26][27][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 I 1 and I 2 both represent inactive states. The transitions between states can depend on Ca 2+ and IP 3 concentration p. As such, the IP 3 receptor will be in various states depending on the value of p.
In this model, the receptor open probability is given by 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 Ca 2+ 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 Ca 2+ concentration and IP 3 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 Ca 2+ oscillations in a sympathetic neuron. Thus, the flux through the RyR is given by To more accurately model CICR, the rate constant k 3 is defined as an increasing sigmoidal function of intracellular Ca 2+ concentration and takes the form where k 1 , k 2 , and k d are parameters. We use n = 3 to match the experimental results obtained in [29]. The parameter k 1 in (7) is the zero calcium concentration level leak. This term is often used to ensure a physiologically significant resting Ca 2+ level [18]. The term k d corresponds to the RyR channel sensitivity for CICR, and k 2 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 Ca 2+ flux through the RyR in the presence of Aβ are minimal. Leak, membrane pump, and SERCA. The membrane leak J in is modeled using a linearly increasing function of IP 3 concentration [30]. Although this increase may be due to various mechanisms, here we only include a linearly increasing contribution to make sure that steadystate Ca 2+ concentration depend on p. Thus, where a 1 and a 2 are parameters. When modeling Ca 2+ 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 where V pm is the maximal velocity and K pm 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 to model the SERCA pump. In (10), K i 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 Ca 2+ in the ER [32,33], and this model provides a way to account for both c and c e . Voltage dependent calcium channel. The flux due to changes in membrane potential is given in our model by J vca . When the volume of the cell is constant, Ca 2+ flux and current are related by the equation where I ca is the current through the VGCC, F is Faraday's constant, and w is the cell volume.
Since I ca depends on the membrane potential V and Ca 2+ , we need a way to track membrane potential. Here, we assume that the membrane potential follows a Hodgkin and Huxley formulation [34] that satisfies where C m 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 C m dV dt ¼ À I kir ðVÞ À I na ðVÞ À I l ðVÞ À I ca ðVÞ þ I app ; ð13Þ where I kir , I na , I l , and I ca correspond to an inward rectifying potassium, sodium, leak, and Ca 2+ current, respectively, and I app is an applied current. In (13), we utilize a formulation for I kir (V) similar to [38] and write where g kir is the maximal channel conductivity, V ka is the Nernst K + potential, K 0 is the extracellular K + concentration, and V a1 , V a2 and V a3 are constants. In the model of [38], K 0 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 K 0 to be constant. Each of the remaining currents have the form where g x is the conductance and V x is the Nernst potential of each channel. Although many alternative formulations for I ca exist (such as the Goldman-Hodgkin-Katz current approximation), we use a T-type like channel (as that found in [39]) of the form where " g caT is the maximal conductance, V ca is the Ca 2+ Nernst potential, and m caT and h caT are gating variables similar to those used for the Na + channel in the Hodgkin and Huxley model.
A T-type Ca 2+ current has a low-threshold activation with an inactivating gating variable that exhibit sub-threshold oscillations and inhibitory rebound bursts. Although other types of Ca 2+ currents (such as L-type, Ca 2+ -dependent potassium channels, etc.) can also be included in the model, here we simply illustrate how one can link membrane potential with other Ca 2+ regulatory mechanisms together into a single model and use that model to investigate the role of Aβ on Ca 2+ dynamics. As such, the T-type Ca 2+ channel described in (15) is sufficient for studying how changes in membrane potentials could impact intracellular Ca 2+ . 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ proliferation as a result of Aβ could trigger the onset of AD. In this study, we are interested in the generation of aberrant Ca 2+ 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 Ca 2+ . 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 Ca 2+ 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 β a m in J in and write an altered membrane leak flux as 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 Ca 2+ signals, the mechanisms for how this happens remains controversial [43][44][45]. Several studies have addressed the role of RyR-regulated Ca 2+ disruptions in AD [46][47][48][49]. Paula Lima et al. (2011) report that Aβ can generate prolonged Ca 2+ signals in vitro through RyR in primary hippocampal neurons of rat embryos as a result of NMDAR-dependent Ca 2+ entry through the plasma membrane. These results agree with studies showing that Aβ can cause substantial Ca 2+ influx through NMDAR [50,51]. In our model these types of fluxes are combined and accounted for in J in , 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 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 Ca 2+ model described above tracks changes in cytosolic Ca 2+ concentration as a function of time in the presence of Aβ. By putting the different channels, exchangers, and pumps together, our Ca 2+ model takes the form Notice that this model has twelve dynamic equations, with five of those controlling the IP 3 receptor (R, O, A, I 1 , and I 2 ), three controlling the membrane potential (V, m, and h), and an additional two variables controlling the gating variables m caT and h caT (a description of these additional equations are provided in the Supplementary Appendix). We have also added a scaling and control parameter p s to simplify the contribution from the VGCC flux. When p s = 0, the impact of membrane potential on Ca 2+ signaling are ignored. When p s > 0, changes in membrane potential do influence the dynamics of intracellular Ca 2+ . 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 Table 1. Parameter values of the Ca 2+ model (18) and (19).
RyR 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 IP 3 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 IP 3 production. Lastly, we include the effects of membrane potential on Ca 2+ dynamics by incorporating flux through VGCC under various levels of fixed IP 3 . 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 c 0 = 0.05 and c e = 10 with R 0 = 1 and all other IP 3 receptor initial conditions are set to zero. Initial conditions for the membrane potential are V 0 = −65, m 0 = 0.05, n 0 = 0.32, h 0 = 0.6, m caT0 = 0.29, and h caT0 = 0.01.

Results
In all three sections below, the model dynamics show that aberrant Ca 2+ 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 IP 3 concentration is fixed and persists and is not influenced by membrane potential. We then allow for IP 3 to change dynamically and simulate the response in Ca 2+ 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 IP 3 concentration is fixed.

Calcium model with constant IP 3 and no membrane potential
Simulations of (18) and (19) with p s = 0 show that incorporating Aβ in a simplified Ca 2+ model can lead to aberrant Ca 2+ 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 Ca 2+ steady-states and oscillatory patterns. From a phenomenological perspective, a transition through MMO may explain why Aβ can trigger persistent aberrant Ca 2+ signals. Although the long-term accumulation and influence of Aβ remains difficult to model, our approach illustrates that aberrant Ca 2+ signals can occur under the influence of Aβ for a variety of parameter regimes.
IP 3 influence on model dynamics. Experimentally, IP 3 can be photoreleased simultaneously through out a cell. In such conditions, IP 3 diffusion is minimized and can be treated as constant. By varying the amount of IP 3 available in the cytosol, the model can exhibit Ca 2+ 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 Ca 2+ 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 Ca 2+ signals. We hypothesize that aberrant Ca 2+ 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][56][57].
Illustrated in Fig 2A and 2B are numerical simulations of the model, with a = 0 that exhibit sustained Ca 2+ 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 Ca 2+ signaling since it is periodic. Notice that no Ca 2+ oscillations occur for small and large values of p. This is due to the assumption that the IP 3 receptor is unlikely to open at high and low concentrations of IP 3 . Although our model does not exhibit a typical Hopf bubble, it has the overall qualitative patterns of a robust Ca 2+ model as described in [18].
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 Ca 2+ 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 IP 3 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 IP 3 is present, the accumulation of Aβ increases the steady-state of Ca 2+ and can even lead to large-amplitude Ca 2+ 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 HB a1 and HB a2 . 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 Ca 2+ 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.
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 p s = 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     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 Ca 2+ 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, Ca 2+ oscillations have much larger amplitudes than in the absence of Aβ.  Contribution of Aβ on calcium signaling through the ryanodine receptor. To better understand the influence of Aβ on Ca 2+ 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 k a . Our simulations show that in the presence of Aβ, changing the parameter k α can produce aberrant Ca 2+ 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 k 2 , the maximal kinetic rate of the RyR.   Fig 7D with k α = 1 when k 2 = 0.18 is increased to k 2 = 0.65. Notice that although solutions do settle into stable periodic orbits that the amplitude of the signals increase. Thus, increasing the parameter k 2 may help to stabilizes aberrant oscillations at the cost of increasing Ca 2+ 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 k 2 .
Dynamic levels of IP 3 with no membrane potential. The model dynamics exhibited in the figures above are triggered by a constant value of IP 3 in the cytosol. Recall that in vivo, IP 3 production is typically a result of an agonist activated G-couple protein. The rates of IP 3 production and degradation are both modulated by intracellular Ca 2+ , and as such, the release of Ca 2+ through the IP 3 receptor can directly alter the IP 3 signaling pathway. Furthermore, there is growing evidence that Aβ also affects the IP 3 signaling pathway [7]. Thus, we extend our model to include a dynamic variable for IP 3 production and degradation, and look to include the influence of Aβ on this signaling mechanism when p s = 0. More specifically, we include an additional equation for IP 3 and model the influence of Aβ on the production of IP 3 .
We make use of the hybrid model formulated by Politi et al. (2006) to track IP 3 production and degradation. Their model takes the form where The first term on the right hand side of (20) corresponds to the production of IP 3 and the second term is the degradation. In (20), V PLC is the maximal production rate, K PLC characterizes The two dots on the right side represents the locations of the parameter values of k 2 used to generate the Ca 2+ traces in Fig 7D and Fig 8C. B and C show the stable periodic oscillations that occur when k 2 is increased to k 2 = 0.5 and k 2 = 0.65, respectively. https://doi.org/10.1371/journal.pone.0202503.g008 Model of amyloid beta on calcium the sensitivity of PLC, K 3K is the half saturation constant for the degradation term. The constants k 3K and k 5P correspond to the IP 3 phosphorylation and dephosphorylation rates, respectively. K PLC and η are parameters used to adjust the positive and negative feedback of Ca 2+ , 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 IP 3 process into two components: a production and a degradation. This will make it easier for us to incorporate the effects of Aβ on the IP 3 production process.
In their experiment, Demuro and Parker (2013) showed that introducing Aβ directly into Xenopus oocytes causes an increase in Ca 2+ dependent fluorescence (a measure for the amount of intracellular Ca 2+ ). Even though their experiments are in oocytes, the ubiquitous properties of IP 3 signaling may make their results relevant to other cells including neurons. Their findings suggest that Aβ does not interact directly with the IP 3 receptor, but instead they propose that intracellular Ca 2+ liberation evoked by Aβ involves opening of IP 3 receptors as a result of stimulated production of IP 3 via G-protein-mediated activation of PLC. As such, in the presence of Aβ, IP 3 are actively stimulated and persist for many minutes or hours even though IP 3 is metabolized within tens of seconds [58]. Based on these findings, we assume that V PLC takes the form and K PLC can be written similarly as 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 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 Ca 2+ 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 IP 3 = 0.01 shows a small Ca 2+ 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 IP 3 is dynamic and depends on Ca 2+ . It will take sometime for enough IP 3 to be present in order to trigger a signaling event throught the receptor. Fig 9A and 9B show the model solution for Ca 2+ and IP 3 , respectively in the absence of Aβ. Notice that the amount of IP 3 also reaches a steady-state level. Since we are interested in understanding the potential effect of Aβ on Ca 2+ 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

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 Ca 2+ 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 Ca 2+ through the VGCC are included in the flux term J vca and as such, allows us to study the impact of membrane potentials on Ca 2+ 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.
To investigate the impact of a constant current pulse on Ca 2+ signaling, we first simulate (18) and (19) with constant values of IP 3 . In these simulations we set p s = 1 and plot Ca 2+ concentrations as a function of time. Tracking the membrane potential and including it into the model will have an effect on Ca 2+ 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).
Notice that in Fig 12D-12F Ca 2+ signals seem to stabilize in single-mode oscillations upon the release of the stimulus before transitioning back into MMOs. This may provide one  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, Ca 2+ 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 Ca 2+ 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, Ca 2+ 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 singleto 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.
Because the impact of membrane potentials are controlled by p s , we have also provided model simulations when this parameter is altered. Small values of p s correspond to small influence of membrane potential while larger values can be used drive the amplitude of Ca 2+ signals. suggests that the effects of membrane potentials are large enough to alter the amplitude of Ca 2+ oscillations during the applied current. Fig 14C shows that when p = 1, contributions from membrane potentials increase overall Ca 2+ 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.
Although much analysis remains to fully understand the dynamics of the model, the results of our simulations suggest that Aβ can alter Ca 2+ 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 Ca 2+ 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 Ca 2+ regulation.

Discussion
Intracellular Ca 2+ is a critically important second messenger within the nervous system. In neurons, Ca 2+ is known to mediate the signaling pathways that control neurotransmitter  (18) and (19) when p s = 0.5, p s = 0.75, p s = 1, and p s = 1.5, respectively. Notice that B shows that Ca 2+ can enter aberrant oscillatory patterns when the membrane is stimulated by a constant applied current. C shows that Ca 2+ signals can enter MMOs with altered amplitudes when the applied current stimulus is turned off.
https://doi.org/10.1371/journal.pone.0202503.g014 release, gene expression, metabolism, plasticity, development, proliferation, and cell death [60]. As such, Ca 2+ may play a major role in the pathogenesis of AD. Unfortunately, the complexity of Ca 2+ 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 Ca 2+ regulation and what leads to aberrant signaling. One of our goals is to model Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ levels. Although we have not considered mitochondrial effects, we do track ER Ca 2+ . Further analysis that looks at the time evolution of c e could be useful in predicting chronic changes of Ca 2+ homeostasis in AD.
Developing a whole-cell Ca 2+ model that has predictable qualitative structure in the presence of Aβ is challenging. Although Aβ influences many Ca 2+ 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 Ca 2+ 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 Ca 2+ signals in astrocytes [61,62]. In these astrocytes, Ca 2+ 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][68][69][70]. However, synaptic transmission is typically precipitated by a presynaptic potential which allows Ca 2+ ions to flow into the cell through VGCC. The contributions of fast local Ca 2+ signals with slow global Ca 2+ 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 Ca 2+ driven mechanisms during synaptic transmission. The simplified whole-cell Ca 2+ model presented here could be linked with a synaptic Ca 2+ model to investigate how global aberrant Ca 2+ signals may impact synaptic transmission on multiple timescales. Simulations over long timescales may help explain how slow global whole-cell Ca 2+ signals interfere with fast local Ca 2+ 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 IP 3 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 IP 3 (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 Ca 2+ 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 Ca 2+ 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 Ca 2+ signals can occur in a simplified wholecell 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 Ca 2+ 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 Ca 2+ levels. The model can then be used to better understand the impact of Aβ on Ca 2+ fluxes through individual regulatory components (such as IP 3 , 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.