Quantifying the dose-dependent impact of intracellular amyloid beta in a mathematical model of calcium regulation in xenopus oocyte

Alzheimer’s disease (AD) is a devastating illness affecting over 40 million people worldwide. Intraneuronal rise of amyloid beta in its oligomeric forms (iAβOs), has been linked to the pathogenesis of AD by disrupting cytosolic Ca2+ homeostasis. However, the specific mechanisms of action are still under debate and intense effort is ongoing to improve our understanding of the crucial steps involved in the mechanisms of AβOs toxicity. We report the development of a mathematical model describing a proposed mechanism by which stimulation of Phospholipase C (PLC) by iAβO, triggers production of IP3 with consequent abnormal release of Ca2+ from the endoplasmic reticulum (ER) through activation of IP3 receptor (IP3R) Ca2+ channels. After validating the model using experimental data, we quantify the effects of intracellular rise in iAβOs on model solutions. Our model validates a dose-dependent influence of iAβOs on IP3-mediated Ca2+ signaling. We investigate Ca2+ signaling patterns for small and large iAβOs doses and study the role of various parameters on Ca2+ signals. Uncertainty quantification and partial rank correlation coefficients are used to better understand how the model behaves under various parameter regimes. Our model predicts that iAβO alter IP3R sensitivity to IP3 for large doses. Our analysis also shows that the upstream production of IP3 can influence Aβ-driven solution patterns in a dose-dependent manner. Model results illustrate and confirm the detrimental impact of iAβOs on IP3 signaling.


Introduction
Alzheimer's disease (AD) is a devastating neurological illness affecting around 40 million people worldwide. AD is the leading cause of dementia, and while the prevalence is estimated to triple by 2050 [1], no cure currently exists. The progressive accumulation of intracellular Aβ in its soluble oligomeric forms iAβOs has been indicated as the leading event in the pathogenesis of AD [2][3][4]. Aβ is a 36-43 amino-acid-long peptide cleaved from the amyloid precursor protein (APP) by β-and γ-secretase. In neurons, cleavage of APP takes place when γ-secretase forms a complex with presenilin (PS) within the ER membrane, where production of Aβ 42 is more likely to occur [5]. Aβ monomers tend to aggregate into soluble oligomers, fibrils, and plaques [6]. This aggregation occurs as the production of Aβ increases faster than can be degraded naturally [7,8]. Aβ accumulation has been shown to occur as a result of multiple factors including overproduction of Aβ and aging-related changes in its clearance mechanisms; both by neuroglia and the lymphatic system [9,10]. Importantly, the accumulation of intracellular Aβ has been shown to precede the appearance of extracellular amyloid plaques and intracellular neurofibrillar tangles associated with tau proteins, suggesting an early role of soluble Aβ during the progression of AD [7,[11][12][13]. The ability of extracellular applied Aβ oligomers to induce cytosolic Ca 2+ fluxes generated from both extracellular and intracellular sources has been shown using cultured mammalian cells [14][15][16]. We have subsequently characterized these two mechanisms as occurring by: i) formation of plasma membrane Ca 2+ permeable pores [17], and ii) permeation of Aβ oligomers into the cytosol and inducing a PLC-dependent Ca 2+ release from the ER [18]. As a critical secondary messenger, Ca 2+ mediates the signaling pathways that control several neuronal processes including neurotransmitter release, gene expression, metabolism, plasticity, development, proliferation, and cell death [19,20]. Furthermore, accumulation of Aβ in neurons has been shown to disrupt intracellular Ca 2+ homeostasis inducing mitochondrial stress [21,22]. Because Aβ accumulation has been shown to alter intracellular Ca 2+ levels, studying its impact on Ca 2+ regulatory mechanisms is critical for better understanding the pathogenesis of AD.
Intracellular Ca 2+ regulation involves many distinct mechanisms working together. In the presence of Aβ, these Ca 2+ regulatory mechanisms begin to fail [22,23]. For example, the presence of Aβ has been shown to increase Ca 2+ liberation from the ER through 1,4,5-Inositoltriphosphate receptors (IP 3 Rs) and ryanodine receptors (RyRs) [15,24]. Aβ can also spontaneously form Ca 2+ -permeable pores in the plasma membrane [20,25] creating uncontrolled influx of Ca 2+ through the membrane. These alterations can cause stress on the ER that can further lead to dysregulation of Ca 2+ in a feed-forward cyclical pattern [15,22,26,27]. Such breakdowns in regulation can create aberrant or sustained elevated Ca 2+ signals that can lead to cell death [14,18].
As Aβ has been shown to affect numerous intracellular pathways, it is difficult, if not impossible, for experimentalists to investigate independently and simultaneously each of these pathways in a complex neuronal environment. Mathematical and computational approaches can offer a supplementary approach to studying the pathology of AD and the impact of Aβ on cellular mechasims. Theoretical models that can consider the impact of Aβ on multiple pathways simultaneously and independently can provide valuable information for designing future experiments and possibly suggesting therapeutic targets. However, before such models can be constructed, developing dedicated models to investigate each proposed pathway involved in Aβ toxicity is crucial. To this point, our goal is to construct a data validated model that can quantify how Aβ interacts with the IP 3 signaling cascade and its consequential disruption of intracellular Ca 2+ homeostasis. Our single cell model provides important advantages toward the development of a whole-cell model, specifically allowing the study of Aβ in a cause and effect manner.
In our previous study, we have shown that intracellular injection of synthetic Aβ 42 oligomers (Aβ 42 Os) into Xenopus oocytes triggered a PLC-dependent activation of IP 3 Rs in the ER membrane causing cytosolic Ca 2+ rise [18]. However, experimental limitations make it difficult to precisely describe the molecular mechanisms involved. As such, we develop a mathematical model to identify and quantify the molecular mechanisms by which Aβ affects IP 3 production and subsequent Ca 2+ release through IP 3 Rs. We first build a computational model capable of tracking intracellular changes in Ca 2+ concentration as a function of time. We assume that intracellular Aβ 42 Os (iAβ 42 Os) have a direct impact on G protein activation and PLC-mediated IP 3 production. The experimental results in [18] provide data to calibrate our mathematical model and to test our modeling assumptions. We show that increasing iAβ 42 Os from small to large doses causes significant changes in the impact of Aβ on certain cellular mechanisms. Our model analysis substantiates that iAβ 42 Os have a widespread effect on IP 3 -mediated Ca 2+ signaling.
Because experimental recordings of Ca 2+ signals are typically expressed as a ratio of fluorescence relative to the resting fluorescence before stimulation (Δf/f 0 ), we use the conversion methodology outlined in Maravall et al. (2000) [28] to directly compare our simulation results with experimental data. We further explore the implications of such conversion on model solutions and provide a detailed analysis of the impact of various model parameters along with predictions showing how the upstream mechanisms in IP 3 production impacts Ca 2+ signaling. Because model kinetics and parameters are linked to certain biophysical mechanisms, we use the model to study how changes in G protein and PLC activation rates impact Ca 2+ signals. We also explore how large doses of iAβ 42 Os alter the sensitivity of IP 3 Rs. Our results provide insight into which cellular mechanisms could become potential therapeutic targets for treating AD. Although Aβ can take many forms, in this work, we solely focus on iAβ 42 Os, positively recognize by OC antibody and simply refer to them as Aβ for simplicity [6,14].

The closed-cell model development
To investigate the impact of Aβ on Ca 2+ regulation, we make use of the vast literature on calcium dynamics including the Ca 2+ signaling "toolkit" [29][30][31][32]. We use experimental conditions and data from Xenopus oocytes to build a Ca 2+ model using traditional methods of 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+ can be modeled by where J denotes flux across internal and external membranes. While various pumps and channels exist between the ER and cytosol in neuronal and glial cells, intracellular Ca 2+ signaling in Xenopus oocytes is mostly due to IP 3 Rs as oocytes are deficient in RyRs. In an in vivo environment, both the Na + /Ca 2+ exchanger and the plasma membrane Ca 2+ ATPase pumps affect Ca 2+ removal from the cytosol while receptor-operated Ca 2+ channels lead to Ca 2+ entry into the cytosol from external sources. The experimental data on which we build the model are performed by monitoring the temporal evolution of the fluorescence signal generated by the bounding of cytosolic Ca 2+ to the Ca 2+ -dependent fluorescent dye. As such, the data extracted from our experiments intrinsically take into account the endogenous activity of the Na + /Ca 2+ exchanger, the plasma membrane and SERCA Ca 2+ ATPase pumps in the absence of specific blockers.
Based on these conditions, we write channel (such as a Receptor Operated Channel), and Plasma Membrane pump, respectively. The constant α is typically used to control the rate of transport of Ca 2+ across the membrane to that across the ER.
Let c e denote the concentration of ER calcium. With this, we assume a Ca 2+ model of the form where γ is the ratio of cytoplasmic volume to ER volume. Note that we do not explicitly consider the effects of Ca 2+ buffers. We assume that Ca 2+ buffers are fast, immobile, and of low affinity (see [30,32,33] for further details on buffering). As such, Ca 2+ buffering is implicitly included in the model by assuming that all fluxes are effective fluxes.
In our modeling analysis we assume that the contributions of J IN and J PM are small compared to the contributions of the ER. As such, we set J IN − J PM � 0 and reduce the model to a closed-cell model where Ca 2+ transport only occurs between the ER and cytosol. Understanding that stable Ca 2+ oscillations in Xenopus oocytes occur in the absence of external Ca 2+ suggests that Ca 2+ exchange with the extracellular environment plays a minor role in the dynamics. However, this simplification does affect the biological implications and the model's ability to describe Ca 2+ regulation in general, and specifically in glial cells and neurons. For example, the direct exclusion of specific contributions from J IN and J PM may over-simplify Ca 2+ solutions as the cell moves away from steady-state conditions. Furthermore, as cells are injected with Aβ, the contributions of the membrane transport mechanisms will certainly affect cytosolic Ca 2+ concentration even in the absence of extracellular Ca 2+ . Regardless, the simplified deterministic model does allow us to illustrate important dynamical properties of Ca 2+ signaling patterns with minimal components. Accordingly, our closed-cell model assumes that Ca 2+ flux into the cytosol is only due to the IP 3 R on the ER and flux out of the cytosol is due to an ATPase SERCA pump back into the ER. This simplified system allows us to model Ca 2+ flux as a mean-field approximation process that considers an average over a large number of IP 3 Rs. While such a model can provide a macroscopic perspective across the whole cell, it cannot capture the stochastic nature of individual channel dynamics. However, such a model is appropriate for our goal of analyzing the influence of Aβ on the IP 3 signaling cascade.
The flux terms in Eqs (1) and (2) can be modeled using various formulations, such as a saturating binding rate model for IP 3 R [34,35] and Markov models [36][37][38]. For our purposes, we assume that the flux from IP 3 Rs follows a formulation based on previous models found in [39][40][41]. Thus, we write where k f controls the density of IP 3 Rs, J ER is the leak from the ER into the cytoplasm, and P o is the open probability of the IP 3 R. In Eq (3), the leak term is necessary to balance the ATPase flux at steady state.
Recall that in our experiments, [18], individual cells were bathed in a Ca 2+ free solution. As such, we assume a closed-cell environment with Ca 2+ fluxes occurring only between the cytosol and the ER and set where c t is the total number of moles in the cell divided by the cytoplasmic volume [32]. We then replace the term (c e − c) in Eq (3) To model P o , we use the Li and Rinzel [42] simplification of the De Young and Keizer [39] formulation for the open probability of the IP 3 R where y is the proportion of inactivated IP 3 Rs and p is the concentration of IP 3 present in the cytosol. To model the SERCA pump, we use a Hill function of degree two. Replacing the fluxes in Eqs (1) and (2), we have where K i , for i = 1, . . .5 and k −4 and k −2 are parameters associated with the transition rates between various quasi-steady-states of the IP 3 R (see [32,42], and [30] for details), and V s and K s are the parameters associated with the SERCA pump.
The parameter values used for these equations are given in Table 1 and are similar to those used by De Young and Keizer [39] with modifications to the cellular and SERCA parameters. The choices for the the cellular and SERCA parameters were obtained by fitting the model to various experimental data illustrated in [18]. For these parameters, the dynamics of Eqs (6) and (7) are illustrated in Fig 1A where the steady-state values are shown as a function of p. As p increases, the dynamics illustrate the classic Hopf bubble and transitions from stable steadystates into periodic oscillations dynamics then back to stable steady-states through the Hopf bifurcation points labeled HB1 and HB2. The top and bottom branches of the bubble give the max and min values of the oscillations as a function of p. Shown in Fig 1B are the nullclines corresponding to dc/dt = 0 (in red) and dy/dt = 0 (in green) along with the trace of the solution when p = 0.325. The dashed lines correspond to the nullclines when p = 0 while the labeled solid red and green curves are the nullclines when p = 0.325. The temporal Ca 2+ solution showing periodic oscillations when p = 0.325 is shown in Fig 1C. Also illustrated there is the variable y in red.
In [18], changes in Ca 2+ concentration occur as a consequence of the intracellular injection of Aβ. As such, new IP 3 is synthesized within the cell during the experimental procedure. To account for the IP 3 dynamics, we use the hybrid formulation of Politi et al. [43]. Let p denote the concentration of IP 3 present in the cytosol, then we write where � V PLC is the maximal rate of IP 3 production and depends on agonist concentration, K PLC characterizes the sensitivity of PLC to Ca 2+ , τ p = 1/(k 3k + k 5p ) represents the characteristic time of IP 3 turnover where k 3k is the maximum rate of 3-kinase and k 5p is the dephosphorylation rate, K ip3 k is the half-activation constant for 3-kinase, and η = k 3k /(k 3k + k 5p ). Both K ip3 k and η are used to tune the positive and negative feedback Ca 2+ in the IP 3 metabolism [30]. The term � V PLC will depend on the amount of activated PLC available and we alter the model by writing where PLC is the fraction of activated PLC complexes, and V PLC is the IP 3 maximal rate of production, to account for time evolving active PLC.
To model PLC and G-protein activation, we use a kinematic model due to Bennett et al. [44] and Lemon et al. (2003) [45]. We assume that PLC is the fraction of activated PLC complexes that drive IP 3 production and that G is the fraction of activated G-protein complexes and write where PLC tot and G tot are the total amount of available PLC and G-proteins (assumed fixed), ρ governs the production of active G-proteins, δ is used as a control for background activity, and k a , k b , k c , and k d are rate constants. Notice that the kinetic formulations above are a simplification of the model constructed by Mahama and Linderman, [46], where a more complex set of equations that account for the hydrolysis of GTP to GDP. A summary of that model can be found in [30].

The effects of Aβ
Although Aβ is clearly implicated in the disruption of intracellular Ca 2+ homeostasis, its interaction with individual pumps, channels, and exchangers remains difficult to quantify. In our previous experiments [18], we performed intracellular injections of Aβ oligomers at various concentrations levels. We also show that the injection of Aβ causes rise of cytotoxic levels of Ca 2+ that carry on over time. This cytotoxicity may be due to stress caused by persistent Ca 2+ release through IP 3 Rs. Of particular interest are the spatiotemporal patterns of fluorescence Ca 2+ signals evoked by Aβ at dose of 1 μg/ml. Recordings in different oocytes showed that Aβ led to various Ca 2+ signaling with ranging patterns from slowly increasing to steady oscillations ( Fig 1C-1E in [18]). Furthermore, when concentration levels of 3 μg/ml, 10 μg/ml, and 30 μg/ml were utilized, the time courses of the fluorescence level of Ca 2+ show that the amplitude of the Ca 2+ signals increases, and the latency to onset and peak response time decreases as the amount of Aβ is increased [18]. In addition, the behavior of Ca 2+ signals for the doses above 1 μg/ml exhibited a prolonged time dependence with an increasing rapid decay as the amount of Aβ is increased. To capture the disparate Ca 2+ signals evoked by various doses of Aβ, our model considers both "small" (1 μg/ml or less) and "large" (greater than 1 μg/ml) doses of Aβ. We utilize these results to hypothesize how Aβ impacts various cellular mechanisms in a dose-dependent manner, and how to incorporate Aβ into the model.  model's ability to reproduce dose-dependent experimental results and are discussed in greater detail below.
Our closed-cell model must be able to reproduce slow monotonic increases in Ca 2+ as a result of the introduction of Aβ, as well as give rise to repetitive oscillations and baseline spikes for doses of 1μg/ml. The model must also be able to reproduce and explain how Aβ leads to increasing signaling peak and a decreasing latency to peak of the response for "large" doses ranging from 3-30 μg/ml. To determine the precise mechanisms by which Aβ affects the cellular machineries that regulate cytosolic Ca 2+ , using several antagonists, we suggest that Aβ acts upstream of IP 3 Rs and hypothesize that Aβ stimulates IP 3 production by PLC in a G-proteindependent manner [17]. Our modeling assumptions for incorporating Aβ were developed through a Monte Carlo Filtering process aimed to isolate the impact of Aβ within our model structure. First, we assume that Aβ acts as an agonist for G-protein activation and write where V R is a scalar, K R is the Aβ concentration producing half activation. The term q represents the effects of a current injection at time t = t 1 of Aβ at concentration a and has the form where H is the Heaviside function and e −r(t−t 1 ) represents the decay of Aβ over time. To match the timeframe of the experimental injections, we set t 1 = 2. In [18], Aβ responses were still evident after 10-15 minutes and as such, we assume a slow decay rate for Aβ and fix r = 0.001 in the model. In this representation, we are assuming that Aβ is acting like a G-protein agonist in a similar way as is expressed in [44].
Our second assumption is to alter the maximal rate of PLC mediated IP 3 production to depend on Aβ as follows where V 0 accounts for PLC mediated IP 3 production under normal conditions, V Q accounts for influence of Aβ on PLC-mediated IP 3 production, and K Q is the dissociation constant. The exponent in V PLC corresponds to a Hill coefficient of 2. A key finding based on this model formulation is that in order to match experimental results, PLC activation needed to be tied to Aβ concentrations. This assumption was determined critical for altering the amplitude of Ca 2+ signals in coordination with the time to peak in our filtering process. Various alternative structures for V PLC were explored numerically but those structures were deemed insufficient for generating the experimental Ca 2+ signaling patterns outlined in [18]. As such, we have assumed that the maximal rate of PLC mediated IP 3 production takes the form of Eq (14), but more data is needed to determine whether this assumption actually captures how Aβ alters PLC-mediated IP 3 production.
Altogether, our closed cell model consists of five differential equations with Aβ input driving the system. In summary, where the term q given in Eq (13) simulates the intracellular injection at time t = t 1 of Aβ at concentration a. Base parameters for the IP 3 , PLC, and G-protein equations are given in Table 2. The parameters are separated by the dose of Aβ used in the model. We characterize a dose of 1 μg/ml and smaller as "small" and doses above 1 μg/ml "large". We explain the distinction and need to separate the parameter space based on Aβ dosage below.

Closed-cell model for small doses
In this section we investigate model solutions in relation to the experimental results described in [18] where a small amount of Aβ is used. A current injection of Aβ at dose of 1 μg/ml gives rise to various spatio-temporal patterns in different cells ranging from a steady increase to periodic solutions. Although we are considering 1 μg/ml a small dose, it was sufficient for evoking local puffs and global responses. Our ODE model cannot capture the traveling waves exhibited in the experiments, but we do show temporal Ca 2+ oscillations that form the basis of wave activity. When the model given by Eqs (15)- (19) is simulated using the parameter values given in the Small Doses column of Table 2 with a = 1 μg/ml of Aβ, we are able to reproduce many of the qualitative features illustrated in Fig 1 in [18]. For example, in some oocytes, a dose of 1 μg/ml leads to a slow and steady increase in Ca 2+ signals that persists. Other cells exhibit amplitude increasing oscillations or steady spike-like responses. These types of responses are captured by the model for baseline parameters with slight variation in the cellular parameters. Because Ca 2+ recordings in Fig 1 in [18] come from different oocytes, we justify slight alterations to cellular parameters as a way to account for variations between individual cells. Note that IP 3 R parameters may also vary, but for now we simply focus on the SERCA parameter K s .
where K D = 0.3 is a dissociation constant that depends on indicator properties, and c 0 is the resting Ca 2+ concentration. In addition, to set the initial condition c 0 (and those of the other variables) we first calculate the steady-state value for the parameter set when a = 0. As such, initial conditions for each of the solutions shown in

PLOS ONE
Dose-dependent impact of amyloid beta on calcium regulation Fig 3D-3E show spike-like Ca 2+ pattern that have a smaller frequency when K s = 0.11 and a decreasing amplitude oscillatory solution when K s = 0.125, respectively. A partial bifurcation diagram with K s as the bifurcation parameter is provided in Fig 3F. As the parameter K s decreases from K s = 0.15, a transition from stable fixed points into periodic orbits occurs through a Hopf bifurcation point around HB3 � 0.1242. As K s continues to decrease, solutions will exhibit sustained oscillations with increased amplitude and a decrease in frequency. The dynamics of model solutions are much more intricate than the partial bifurcation diagram in Fig 3F suggests, especially around K s = 0.09 where multiple limit points and Hopf bifurcations exist. However, our goal is not to fully examine the model dynamics but to merely show that by altering a single model parameter, we can generate solutions that are similar to experimental recordings. A complete description of the dynamics in this region is beyond the scope of investigation and is not included in our analysis. K s is the dissociation constant for the SERCA pump and is the Ca 2+ concentration occupying half of the binding sites of the pump. A smaller K s value corresponds to needing less Ca 2+ to attain 50% of the maximal response for the pump. Whether changes in K s are due to Aβ or simply through chance variation in cells remains debatable. Here, we argue that it is alterations in cellular structures modeled through differences in parameters that is causing the changes in Ca 2+ signals and not because of Aβ's direct impact on the SERCA pump. However, more A shows an increasing Ca 2+ signal that settles to an increased steady-state when K s = 0.15. B shows that oscillations in Ca 2+ can exhibit increasing amplitudes such as those found in Fig 1D in [18]. C and D show that as the value of K s decreases, the oscillatory patterns of the model reproduce the spike-like Ca 2+ signals observed in Fig 1E in [18]. E illustrates an oscillatory solution with an increased steady-state Ca 2+ homeostasis level when K s is just above the Hopf point. Both D and E show the traces for c s , y, and p in blue, red, and black, respectively. F shows a simplification of the scaled bifurcation diagram with the bifurcation parameter K s . Notice that as K s decreases from the base value of 0.15, a transition from stable steady-states into periodic oscillations occurs through a Hopf bifurcation around HB3�0.1242. The dynamics around K s = 0.09 include multiple Hopf bifurcations and has more complex structure than what is presented here. analysis is needed to fully understand how different doses of Aβ may influence the generation of various Ca 2+ signals.
When looking at the model, the subsystem given by Eqs (18) and (19) is driving IP 3 through the inclusion of the PLC term in Eq (17). As such, we can investigate the subsystem given by Eqs (15)-(17) by treating PLC as a parameter and fixing a = 1. Illustrated in Fig 4A are the general dynamics of Ca 2+ using PLC as a bifurcation parameter for the subsystem given in Eqs (15)- (17). Notice that the dynamics of Ca 2+ will transition from stable steady-states (solid black curve), at the Hopf bifurcation point HB4�0.0043 (labeled in blue), into periodic solutions until transitioning back to stable steady-states at Hopf bifurcation point HB5�0.0078 (labeled in red). The green Hopf bubble captures the maximum and minimum values of the Ca 2+ oscillations. Fig 4B shows the subsystem solution when PLC = 0.005 for the base parameters given in Table 2 for small doses. Intracellular Ca 2+ signal, the proportion of inactivated IP 3 Rs, and the concentration of IP 3 are given by the blue, red, and black traces, respectively. Since both PLC and a drive the responses of Eqs (15)-(17), a two parameter bifurcation diagram where the location of the Hopf points have been tracked as a function of PLC and a is given in Fig 4C. Notice that as the dose of Aβ gets closer to zero, the location of the Hopf bubble shifts to the right. This implies that in order to observe oscillatory behavior when a is close to zero, the amount of active PLC needs to be greater. In Fig 3, we showed that the model solutions will behave differently as the cellular parameter K s is varied. Here, we also look at the impact of varying K s on the solutions of the subsystem Eqs (15)- (17). Illustrated in Fig 4D is the bifurcation diagram when K s = 0.11. In this case, the bifurcation diagram shows an increased region of oscillations accompanied with increased amplitudes for most of the values of PLC in the range of the plot. Depending on the value of PLC, the oscillations will take more of a spiking form than sinusoidal oscillations, which is important as the signals observed experimentally correspond to spike-like signals of local puffs and global Ca 2+ spikes. The red dashed curves in this figure correspond to unstable oscillations. The two Hopf bifurcations are given by HB6 � 0.0002 and HB7 � 0.02115. Fig 4E shows the subsystem solution when PLC = 0.01 and K s = 0.11. To better understand the impact of changes in both PLC and K s on the dynamics of the subsystem Eqs (15)-(17), a two parameter bifurcation diagram is given in Fig 4F. In this figure, the parameter space is separated into a region where periodic orbits exist and a region where the model has no periodic orbits. The blue curve in this figure tracks the location of the Hopf points generated by the subsystem. For values of K s between approximately 0.1127 and 0.1511 the bifurcation diagram will have a Hopf-like bubble between two Hopf bifurcations (as those illustrated in Fig 4A and 4D). Although the complexities of these bifurcation diagrams varies, the two parameter bifurcation diagram helps us understand the oscillatory nature of solutions when variations in PLC and K s occur. The red cross and triangle shown in Fig 4F correspond to the location of the parameter values for the diagrams generated in Fig 4A and 4B, and Fig 4C and 4D, respectively.
To further investigate the behavior of the small doses parameters, we decouple the PLC and G subsystem given by Eqs (18) and (19) and look at the time evolution of the fraction of active PLC and G. As a is varied for small doses, both PLC and G quickly (on the order of seconds) settle to their new steady-states values. Illustrated in Fig 5A and 5B are the temporal solutions of the subsystem (18) and (19) for a = 0.1 (black), a = 0.5 (magenta), and a = 1 (blue), respectively. Fig 5C shows (18) and (19) will generate oscillatory responses when K s is decreased from K s = 0.15. Although additional analysis can be done for various parameters in the model, we now turn our attention to how altering the dose of Aβ impacts the model solutions.

Dose response relationship between amplitude and latency
Dose-response experiments in Xenopus oocytes demonstrate two major effects on Ca 2+ fluxes following increasing doses of Aβ: the amplitude of the Ca 2+ signals increases with the amount of Aβ and the latency of the maximum peak time decreases as the amount of dose increases. We can test model against the experimental data starting with the small dose parameters to determine how the amplitude and latency of solutions vary as the doses of Aβ are increased. Illustrated in Fig 6A are scaled model solutions for Aβ doses of a = 1 μg/ml (black), a = 3 μg/ ml (blue), a = 10 μg/ml (red), and a = 30 μg/ml (green) using the small doses parameters in Table 2. Notice that as a increases, the model captures both the amplitude increases and the decrease in latency to peak but is insufficient for reproducing the observed Ca 2+ signals for large doses. Using the small dose parameters to study to explore model solutions and investigate the long term behavior of the model is helpful even though our analysis suggests needing two different dose-dependent parameter sets in order to match key experimental observations.
In the short term (on the order of minutes), solutions of the model with the small dose parameters tend to an apparent new homeostasis level. However, since the amount of Aβ introduced in the model through Eq (13) will eventually decay towards zero, the solution will tend back to the original steady-state value. This can be seen in Fig 6B where the model solutions are shown on a timescale of hours with r = 0.001 (the initial peak of solutions have been removed to better illustrate the long-term behavior). Whether Aβ decays naturally or persists in cells may depend on many factors. The Calcium hypothesis for AD suggests that the amyloidogenic pathway remodels the neuronal Ca 2+ signaling pathway responsible for cognition [13,47,48]. As such, a slow accumulation of Aβ may increase the cytosolic Ca 2+ level of cells leading to toxic stress and in turn can feed back into the hydrolysis of the amyloid precursor protein in a vicious cycle. In an in vivo environment, Aβ may slowly transition from small to large concentrations over timescales of months to years. Although any long-term analysis is beyond the current model, this model shows that if Aβ persists in the model (i.e., when r = 0), solutions would tend to new higher dose-dependent steady-state values as indicated by the dashed lines in Fig 6B. As expected, increasing r in the model will cause the solutions to decrease back to the original steady-state more rapidly.
To understand the impact of variations in the parameters on the amplitude and latency of solutions, Fig 6C shows the location of the solution peak as the parameters k f , J ER , and γ are uniformly varied 10% from base values for 100 trials. The dashed black curve in this figure corresponds to the amplitude and latency for the base small doses parameters in Table 2 for Aβ doses between a = 0.1 and a = 40. Notice that the amplitude is more variable than the latency for the dose of a = 30 while the opposite occurs for smaller values of a. Fig 6C is intended to illustrate that the model can capture some of the effects of "large" doses of Aβ as observed in Fig 1G of [18] while using the small doses parameters given in Table 2.
Interestingly, although the model given by Eqs (15)- (19) with the small doses parameters can capture many of the qualitative behaviors observed experimentally, it lacks some important features when large doses of Aβ are introduced. For example, the recorded average fluorescence response for doses of 3, 10, and 30 μg/ml, have a much longer time dependence and display an increasingly rapid decay (see Fig 1G of [18]). These Ca 2+ signals differ from responses evoked by a dose of 1 μg/ml (such as those illustrated in Fig 3). The model solutions shown in Fig 6A do not capture these behaviors and as such cannot fully represent the impact of Aβ on cellular mechanisms (at least for large doses). Although we do not fully understand how large doses of Aβ affects the Ca 2+ signaling cascade, our goal is to use the model to better understand how Aβ may be impacting individual cellular mechanisms through appropriate parameter selection. To do this, we alter model parameters to match the experimental data in Fig 1G of [18] and then use those results to describe the possible role that large doses of Aβ plays in Ca 2+ signaling. In essence, in order to reproduce the observe experimental data when various doses of Aβ are used, we distinguish model behavior through the selection of smalland large-doses parameter sets.

Large doses parameter fitting
The model developed in the previous section tracks Ca 2+ concentration as a function of time. The experimental data in [18] tracks changes in Ca 2+ as a ratio of changes in fluorescence intensity with baseline fluorescence levels. This is often written as δf = (f−f 0 )/f 0 = Δf/f 0 with f 0 representing the fluorescence intensity at resting Ca 2+ concentration. To better understand the impact of Aβ on Ca 2+ dynamics through modeling, we first rescale fluorescence measurements to Ca 2+ concentrations. According to Maravall et al. [28], changes in Ca 2+ concentration are associated with changes in fluorescence through the equation where f max is the intensity of the dye at maximum Ca 2+ concentration, R f = f max /f min is the indicator's dynamic range with f min being the intensity at minimum Ca 2+ concentration, δf max is the saturation of the Ca 2+ indicator, and f m = f max /f 0 . We use Eq (21) to rescale the experimental fluorescence data found in Fig 1G of [18]. Further details regarding the rescaling procedure are provided in the Appendix. With the rescaling procedure described in the Appendix, we now have a way to convert the experimental fluorescence data in [18] to Ca 2+ concentrations and link model solutions with experimental data. We first fix the scaling parameters K D = 0.3, R f = 100, f m = 40 and then determine the value of the model parameters that will evoke the appropriate Ca 2+ signals. The parameters used for the large doses of Aβ are given in Table 2 under the Large Doses column and were determined by fitting solutions to the converted experimental data for each level of Aβ. Starting with the small doses value, each parameter was stochastically chosen from an individual parameter distribution and a least-squares fitting procedure was used to identify a model parameter set corresponding to an approximate minimum of our objective function. We used a random sampling procedure to draw a parameter set q s from an admissible parameter space Q 2 R p (where p is the number of model parameters). The distribution of each parameter was chosen to match those of previous studies whenever possible. We then minimized the objective function where sed(i) is the scaled experimental data value at i, and ic s (i, q) is the corresponding (interpolated) scaled Ca 2+ solution at i. Our minimization technique uses a random sampling procedure with a random walk process when local minima are found. That is, we randomly select parameter values and compute Err. If Err is less than some threshold, we then perform a random walk around the parameter values that generated the local minimum error Err to locate a local minimizer. While minimizing the objective function for a large number of parameter selections provides potentially good estimates for model parameters, we did not analyze the parameter space with the intention of finding a global minimum. Regardless, the minimization technique does provide a way to establish parameter values that otherwise would be difficult to estimate.
To further understand the impact of parameters on model solutions, we also implemented an additional minimization technique where we took individual parameter subsets from Table 2, varied those, then compared the results with the experimental data. For example, starting with the small doses parameters, we only varied the PLC parameters to determine whether changes in those parameters could capture the large doses experimental results, and so on. This process was conducted for many parameter subset combinations starting with the small doses parameter set. Illustrated in Fig 7 are two "best fit" scaled model solutions c s (smooth curve) shown on top of the scaled experimental data (dashed curve) for each of the three Aβ concentrations. Fig 7A shows a best fit solution when Cellular, SERCA, and the IP 3 R parameters are keep fixed. Observe that the "best" fit parameters are not those listed in Table 2 for either the small or large doses since we are varying some parameters and keeping others fixed. Note that the best fit solutions illustrated here are not much different from the solutions shown in Fig 6A. This suggests that alterations in some of the Cellular, SERCA, and IP 3 R parameters appear to be necessary to capture the observed behavior for large doses of Aβ. Similarly, Fig 7B shows a best fit solution when all but the IP 3 R parameters are varied. This simulation is included to clearly illustrate the need for altering all model parameters, particularly the IP 3 R parameters. The results of these simulations demonstrate that Aβ has a pervasive effect on the entire cell structure in large doses since matching the experimental data did require variation in every set of cellular mechanisms included in the model. Illustrated in Fig 8 are model solutions when a = 3, a = 10, and a = 30 using the large doses parameters given in Table 2. We also included the solution for a = 1 to illustrate how this large doses model behaves for the dose of 1 μg/ml for comparison. Fig 8A shows the scaled model solution (smooth curve) on top of the scaled experimental data (dashed curve) for each of the three Aβ concentrations. Fig 8D shows the unscaled Ca 2+ concentration c illustrating that the scaling procedure does not effect the model's ability to capture the general behavior of the Ca 2+ signals observed experimentally. Solutions for p are plotted in Fig 8B and 8E using two different timescales. Again, in our model Aβ decays exponentially and over the course of a couple of hours, the model solutions settle back to their original steady-states. Fig 8C and 8F show the evolution of y using two different timescales. These two figures show that the proportion of IP 3 Rs that are inactivated by Ca 2+ remains fairly high over the course of hours acting to suppress Ca 2+ spikes over time.
All model parameters used in the simulations illustrated in Fig 8 are given in Table 2 under the Large Doses column. Notice that the differences in each solution (as given by the different colors) of Fig 8A is only driven by changes in the value of a. In all simulations, initial conditions were found using the steady-state values when a = 0. Noteworthy, our large doses model is efficiently capable of capturing the increase in amplitude of the Ca 2+ concentration signal and the decrease in latency to peak onset as well as increasingly rapid decay as the Aβ concentration a is increased, agreeing well with high suitability the experimental data for large doses of Aβ. Furthermore, the model with this parameter set is able to capture the slowly increasing Ca 2+ response seen in some oocytes with a dose of 1 μg/ml (such as responses similar to those shown in Fig 3A), but it cannot reproduce the various oscillatory and spiking behavior through small variations in parameters (such as those shown in Fig 3B-3E). The model with the small doses parameters cannot capture the increasingly rapid decay based on Aβ nor the extended time dependence, underscoring the need for two different parameter sets. The difference in parameter values between the two sets suggests that Aβ has a pervasive impact that permeates throughout a cell over time and gives credence that Aβ may indeed be affecting multiple cellular mechanisms simultaneously.

Uncertainty quantification and partial rank coefficient correlation for large doses
As with any experimental procedure, uncertainty in measurement naturally arises within the environment. These variations mean that finding exact values for model parameters is unrealistic. Performing uncertainty quantification allows us to determine how changes in parameter inputs affect model solutions. For example, in [18] Ca 2+ responses are categorized by the change in fluorescent signaling and results are given as an average of 4-5 cells. Responses from individual cells can also change from cell to cell and as such, there could be natural variations in output.
To account for these uncertainty principles we vary the large doses parameters stochastically within 5% and 10% of baseline using a uniform distribution and generate n = 100, 000 solutions to the model. This type of simulation allows us to better understand the robustness of the model and provides some way to assess the influence of parameter selection on model results (see [49] for details on method). With the collection of n sample solution paths, we then compute the mean and standard deviation at each time t.   in Table 2 is not optimal in minimizing our objective function, it does provide a reasonable set even under small perturbations. As such, our simulations convey evidence that the modeling assumptions may help capture how Aβ influences the cellular mechanisms involved in PLCmediated IP 3 production.
To better understand how each parameter impacts model solutions, we use sensitivity analysis based on partial rank correlation coefficients (PRCC). This allows us to determine the statistical relationship between model parameters and the resulting Ca 2+ dynamics [50]. To do this, we characterize the resulting Ca 2+ dynamics with two quantities: the peak Ca 2+ concentration achieved during the simulation and the time at which this peak occurs. The PRCC measures the strength of the linear relationship between each model parameter and the model outcome after correcting for the linear effects of all other model parameters. The resulting PRCC scores take values between −1 and 1 with a negative value indicating that the model outcome decreases as the parameter increases and a positive value indicating that the model outcome increases as the parameter increases. The strength of the relationship between the model parameter and model output is indicated by the magnitude of the score.
The results of the PRCC are given in Tables 3 and 4. Table 3 shows the correlation to peak Ca 2+ concentration while Table 4 shows the correlation of the time of peak. The tables list the Table 3. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the maximum Ca 2+ concentrations for each of the three levels of Aβ. Results are with 10% variation in parameters values. � indicates the correlation coefficient is not significant at the p = 0.05 level. correlation coefficients for each parameter when a = 3, 10, and 30. The ranking of the parameters was done by taking the average of the PRCC for the three doses of Aβ. As such, the parameters that most decrease the peak amplitude of Ca 2+ solutions are the parameters k d , k b , and K 1 while the parameters that most increase the amplitude are V Q , k c and k a , as a is increased. Similarly, the parameters that most decrease the time of peak of Ca 2+ solutions are the parameters V Q , k a , and k c while the parameters that most increase the time peak are K 1 , k d and k b . Although these parameters exhibit the strongest effect, we note that most other parameters exhibit a smaller but significant effect. Our intention is not to give a complete analysis for each model parameter, however we do analyze some interesting behaviors pertaining to specific parameters below. When looking at the PRCC analysis, it appears that the PLC and G-protein rate constants k a , k b , k c , and k d all have a large impact on the solution patterns in terms of solution peak and time to peak. Recall that, k a and k c are the activation rates for PLC and G-proteins, respectively. As the activation rates increase, this will lead to an increase in IP 3 production and you will see the peak of the Ca 2+ signal occur sooner. On the other hand, k b and k d correspond to the inactivation of PLC and G-proteins, respectively. A higher inactivation rate for both PLC and Gproteins will decrease IP 3 production and thus lower the peak amplitude of Ca 2+ responses. Table 4. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the time to peak when maximum Ca 2+ concentration was reached for each of the three levels of large doses of Aβ concentration. � indicates the correlation coefficient is not significant at the p = 0.05 level. From a biological perspective this makes sense, once PLC is activated, the production of IP 3 occurs through hydrolysis of phosphatidylinositol-4,5-biphosphate (PIP2). Thus, as the amount of active PLC increases, we should see an increase in the amplitude peak and a decrease in the time to peak in Ca 2+ responses. Conversely, as the amount of active PLC decreases, we should see a decrease in the amplitude peak but an increase in the time to peak in Ca 2+ responses as fewer IP 3 are available for binding to the IP 3 R. Even though these results align with what one might suspect occurs from a biological perspective, these behaviors are directly linked to how the model was constructed. Specifically, recall that the subsystem given by Eqs (15)- (17) is solely driven by PLC and Aβ. Changes in PLC will play a major role in the amount of IP 3 available for IP 3 R binding. Further analysis on the impact of these parameters is provided below. As noted above, the parameter K 1 also plays a major role in solution patterns. As adapted from the De Young and Keizer (1992) model, this parameter corresponds to the effective binding rate of IP 3 to one of the IP 3 R model subunits when no inactivating Ca 2+ is present. As such, this parameter helps drive the IP 3 R dynamics. In the model, an increase in K 1 has an inactivating effect on the IP 3 R since either the unbinding rate of IP 3 to receptor binding site is increased or the binding rate is decreased. In either case, this would decrease the opportunity for the receptor to remain in an active and open state. The PRCC analysis highlights that K 1 is critical for understanding the Ca 2+ patterns of the model. Because of the influence of this parameter on model solutions, this suggests that the IP 3 R dynamics does contribute to the observed Ca 2+ patterns. We analyze the model below to further expand on the influence of K 1 on model solutions. As the model suggests that changes in K 1 may be dependent on Aβ levels, further investigations on the connection between the IP 3 R and Aβ may be necessary.

Correlation to peak of signal
The PRCC also highlights additional interesting information regarding the influence of specific parameters on model solutions. Interestingly, the dependence of solution amplitude peak with respect to the parameter K s appears to be tied with the size of a. More specifically, the PRCC for K s when a = 3, 10, and 30 are 0.712, 0.514, and −0.147, respectively. This implies that as a increases, altering K s has a different effect on model amplitude. Namely, the amplitude increases for a = 3 and a = 10, but decreases when a = 30. Notice that similar results occur for the parameters K PLC and k 5p but in the opposite direction. The dependence of solution time to peak with respect to the parameter δ also appears to be linked to the value of a. In this case, the PRCC for δ when a = 3, 10, and 30 are −0.494, −0.212, and −0.056, respectively. Although the sign of the PRCC is negative in each case, the disparity of the correlation coefficient may indicate that Aβ is affecting the intrinsic background production of active G-proteins differently as the doses vary. The dependence of a on these parameters suggests that Aβ is impacting the mechanisms differently as the amount of Aβ is altered. Further exploration of these parameters may tease out additional information about the influence of Aβ on cellular mechanisms but is beyond the scope of this study.

Impact of Aβ on IP 3 R for large doses
The impact of Aβ on the IP 3 signaling cascade appears to be concentration dependent. Not surprising, the PRCC analysis suggests that the rates k a , k b , k c , and k d play a significant role on the amplitude of responses and the peak time. These parameters directly influence the amount of PLC that feeds into the subsystem given by Eqs (15)- (17) and small variations in these parameters will greatly affect the solutions of the model. Instead of looking specifically at these parameters, we can alternatively investigate the impact of changes in V Q . Recall that V Q controls the influence of Aβ on PLC-mediated IP 3 production. As such, it is no surprise that V Q also plays a significant role in the solution patterns. Although the PRCC identifies the parameter K 1 for example, as playing a significant role on model solution's amplitude and time to peak, the PRCC analysis cannot capture how variations in a single parameter will affect model solutions in general. For example, it is not evident in the PRCC analysis that the parameter k −4 plays a significant role on solutions and is a critical parameter when considering the large doses Ca 2+ signaling patterns observed experimentally. Varying k −4 has a direct impact on the Ca 2+ signal tail and partly controls the decay of the signals, but does not alter the amplitude or time to peak significantly. Both K 1 and k −4 are parameters that help control the dynamics of IP 3 Rs.
Shown in Fig 11 are two diagrams that illustrate the impact of Aβ on the IP 3 R itself through the parameters K 1 and k −4 when a = 10 (a similar effect occurs for a = 3 and a = 30). Fig 11A  shows the representation of the effects of varying K 1 model solutions. Starting with the large doses parameters, we simulate the model by altering K 1 from the base small doses value of K 1 = 0.13 (bold black trace) and increasing the parameter to the large doses value K 1 = 0.13 (smooth red trace). As is suggested by the PRCC analysis, we see that K 1 is negatively associated with the peak amplitude and positively correlated with respect to the time to peak. The impact of changes to the parameter k −4 is shown in Fig 11B. Similar to Fig 11A, starting with the small doses parameter value k −4 = 0.029 (bold black trace) and decreasing the parameter to the large doses value k −4 = 0.00006 (red trace) shows that k −4 plays a critical role in controlling the decay of Ca 2+ signals. Interestingly, the PRCC does not capture this effect as it was only conducted to track the impact on the amplitude peak and latency of solutions. Altering the other IP 3 R parameters will have various effects on solutions similar to the impact of varying K 1 . Changes to IP 3 R parameters seem necessary in order to capture the increasingly rapid decay and suggests that Aβ for large doses may act to desensitize the IP 3 R. Whether Aβ directly interferes with IP 3 Rs remains debatable but our model suggests that Aβ does indeed alter the receptor dynamics for large doses. There may be some intrinsic threshold on Aβ concentration within the cellular environment for which the sensitivity of IP 3 Rs is affected by Aβ. Of particular interest is the role of the IP 3 R parameters in capturing the observed rapid decay of Ca 2+ signal for large Aβ doses.

Limitations of the model
As with any mathematical model, many limitations exist with the approach presented here. Because of our interest in dissecting the effects of Aβ on the IP 3 signaling cascade, the model development and construction utilized a number of simplifying assumptions. While many of these assumptions are traditional, the simplistic nature of the model cannot fully represent the biological environment in a holistic way. None-the-less, our approach has sought to balance the complex biophysical mechanisms involved in Ca 2+ signaling with that of a mathematical structure that can be useful in identifying key factors involved in generating certain solution patterns. Unfortunately, a lack of data has made it difficult to determine the precise conditions and the validity of many of the modeling assumptions. For example, we acknowledge that the steady-state assumptions and the particular mechanisms for how Aβ may be interfering in the Ca 2+ signaling process need to be explored further. Although these assumptions contributed to model solutions whose behavior and dynamics match experimental results, more data is needed to fully justify these assumptions. Additionally, the inclusion of other Ca 2+ regulatory mechanisms will be necessary to describe whole-cell calcium dynamics in a biologically robust way.
Our model construction assumes that iAβ 42 Os (1) act as an agonist for G-protein activation, and (2) affect the maximal rate of PLC mediated IP 3 production. The second assumption was developed based on the results of a series of Monte Carlo numerical simulations that considered a wide-array of possible sites for including the impact of iAβ 42 Os on cellular mechanisms. These simulations were conducted using a large number of initial parameter sets and a variety of functional representations (such as Hill functions of various degrees). Although we were able to match some of the observed experimental results for large doses without including the assumption given in Eq (14), we could not reproduce the three Ca 2+ signals (a = 3, 10, and 30) with the same parameter set simultaneously. Furthermore, any parameter set that closely matched the changes in amplitude and time to peak for small doses of a could not reproduce any spiking behavior observed through cellular and SERCA parameter variations unless V Q 6 ¼ 0. That led us to incorporate the Aβ-dependent term for the maximal rate of PLC mediated IP 3 production given in Eq (14). Due to the complex dependence on model parameters, it may be that this model assumption does not accurately capture how Aβ interferes with the IP 3 production pathway. However, it proved valuable in reproducing observed data for both the small and large doses and provides a possible avenue for further investigations.
As with any model involving numerous parameters, solutions will vary based on the parameter set utilized. In this work, we first rescaled the experimental data, then fitted our model using a best fit parameter estimation procedure. When alternative scaling parameters are used, the model parameters will necessarily change. However, our results show that the model captures the changes in the amplitude and peak time of the signals in a robust and predictable way for both small and large doses of iAβ 42 Os. The PRCC analysis also provides a structured way for understanding how each individual parameter impacts model solutions. Further analysis of our PRCC results could bring to light additional parameter and Aβ-related dependencies. For example, the PRCC values for some parameters are highly dependent on Aβ concentration. Such parameters may also play an important role in determining the possible kinetic interaction of Aβ within the IP 3 production cascade.

Discussion
Ca 2+ is one of the most versatile and universal signals in the human body playing a pivotal role in controlling numerous aspects in the physiology and biochemistry of neurons [51]. Accordingly, intracellular Ca 2+ dysregulation has been implicated in a wide variety of immunological disorders and neurodegenerative diseases including Alzheimer's, Parkinson's, and Huntington's disease. In neurons, as in many other cell types, IP 3 -mediated elementary Ca 2+ signals, also referred to as puffs, are the building blocks of cellular Ca 2+ signaling, and arise through the concerted opening of clustered IP 3 Rs coordinated via a Ca 2+ -induced Ca 2+ -release mechanism [52]. Although the cytosolic Ca 2+ dependency of IP 3 Rs has been well characterized, little is known as to how changes in basal cytosolic [Ca 2+ ] would alter the dynamics of IP 3 -evoked Ca 2+ signals in disease cells, such as neuronal cells of Alzheimer's and Parkinson's disease brains. In AD, iAβOs are now believed to play a major role in the early phase of the disease as their intracellular rise correlates well with the symptoms of AD [3,53]. More generally, AβOs have been found to be predictive of cognitive status at death among patients with AD [54]. Various mechanisms have been proposed to correlate the progressive intracellular Ca 2+ elevation with the concomitant increase of iAβOs observed in neurons during the progression of the AD [25]. Among them, the detrimental activity of iAβOs on the normal functioning of the IP 3 -signaling pathway has been indicated as a potential mechanism responsible for alteration of the Ca 2+ homeostasis in AD neurons.
We and others have suggested that a G-protein mediated activation of PLC by iAβ 42 Os is responsible for the overproduction of IP 3 and consequent rise of cytosolic Ca 2+ observed in cells exposed to iAβ 42 Os [14,18]. Moreover, others have suggested that Aβ may cause cytosolic Ca 2+ rise by a mixed mechanisms of PLC-dependent and independent manner [15,16,55]. The effect of iAβ 42 Os on intracellular Ca 2+ fluxes have previously been investigated by developing a computational model to study important intracellular Ca 2+ pathways in normal and in iAβ 42 Os affected conditions [27]. However, no upstream IP 3 production processes were incorporated in the model. Here, we have illustrated a possible mechanistic way for how iAβ 42 Os triggers IP 3 overproduction with consequent rise in cytosolic Ca 2+ by including some mechanisms of upstream IP 3 production in the model. Specifically, we pinpoint two main possible sites of action for iAβ 42 Os to interact in the cascade of events resulting from stimulation of Gprotein in the plasma membrane to the release of Ca 2+ from the ER.
In our previous study [18], we argued that it was unlikely that iAβ 42 Os act on IP 3 Rs in the generation of Aβ-induced Ca 2+ signaling events. The results of the model are consistent with this for the small doses parameters. However, the model also suggests that iAβOs may in-fact be directly affecting the IP 3 Rs when large doses are introduced. The analysis illustrated in Fig  11 helps us understand what happens to Ca 2+ signaling in the presence of iAβ 42 Os as changes to IP 3 Rs occur. The persistent increase of iAβ 42 Os may alter the sensitivity of IP 3 Rs to Ca 2+ over time. For large doses of iAβ 42 Os, IP 3 Rs may become more sensitive to low-or sub-threshold IP 3 levels and in turn trigger local and global Ca 2+ signaling events. The fact that the parameter k −4 appears to play a major role in the decay of observed Ca 2+ signals singles out the potential that iAβ 42 Os do act on the IP 3 R itself, at least for large doses. Our model suggests the need for further investigation on the relationship between iAβ 42 Os and the sensitivity of IP 3 Rs to IP 3 levels.
Our approach provides a precise way to incorporate the effects of iAβ 42 Os on IP 3 signaling mechanisms that does not necessarily depend on the choice of the IP 3 R model. When a saturating binding rate model for the IP 3 R model is used (as that used in [33] instead of the Li-Rinzel formulation), such a model can also capture the changes in amplitude and peak times for large doses using the same upstream modeling assumptions as outlined above (unpublished results J. Latulippe). This provides further justification that the modeling kinetics of the possible interactions of iAβ 42 Os with G proteins and PLC may be sufficiently captured by the model. Additionally, Toglia et al. [27] have also suggested a relationship between IP 3 concentration and iAβ 42 Os. However, their investigation assume that IP 3 concentration levels are impacted by iAβ 42 Os but use a data fitting procedure to do this rather than attributing those changes to upstream mechanisms. As such, we believe that the model presented here is the first to quantify possible mechanisms for how iAβ 42 Os affects the upstream mechanisms in the IP 3 signaling cascade.
Although our model considers the impact of iAβ 42 Os specifically on the IP 3 signaling cascade in oocytes, our results could be useful in more complex models of various cells. Existing astrocyte models (such as [34,[56][57][58]) that incorporate Ca 2+ dynamics could be altered to include the effects of Aβ on IP 3 signaling components described in this study. This would provide a way to test model assumptions and determine whether solution patterns are consistent in different model environments. Furthermore, the current model could be expanded to include additional pumps and channels known to play a role in various cell types. Incorporating data driven models within the Ca 2+ modeling toolbox may prove to be an efficient way to develop whole cell models that can be used to study how Aβ alters various signaling pathways. For example, the ability to express exogenous proteins, including NMDA Receptors, provides a powerful tool as a possible next step in developing increasingly elaborate mathematical models capable of more closely mimicking neuronal behavior.
Because of the complex cross-talk nature of Ca 2+ signaling, our model also provides a way to control for and test various therapeutic strategies in a modeling environment. For example, to mimic the intrinsically slow accumulation of Aβ seen in the pathology of AD, Aβ can be introduced very slowly into the model and solutions simulated accordingly. We can then introduce artificial agonists or antagonists that affect G-protein activation and PLC function to see how they affect Ca 2+ signals over various timescales. Using the model to better understand what happens to Ca 2+ regulation in these simulations can directly influence and suggest how one could control Ca 2+ signaling in the presence of Aβ, and more generally, various AD environments.
The results of this study suggest the need for two different dose-dependent models to incorporate changes in cellular Ca 2+ signaling in the presence of increasing concentrations of iAβ 42 Os. In in vivo environments, it may be the case that in the early phase of AD, slowly accumulating levels of iAβOs remain relatively small. Under such conditions, the small doses model may be better suited than the large doses model. Regardless, our model development and analysis suggests that increasing the amount of iAβ 42 Os present in the cell can have a pervasive impact on numerous cellular mechanisms.
Building computational models can help provide a better understanding for the complex cross-talk between various signaling mechanisms within neurons, something difficult to establish with current experimental capabilities. Through further analysis and development, researchers can use the model to formulate novel experimental procedures and eventually suggest new therapies for treating AD.
As can be seen from Fig 12, R f has little effect on the value of δf max . This is consistent with the idea that for indicators with a large dynamic range, the exact value of R f is insignificant [28]. For indicators such as Fluo 4, K D is often assumed to be between 0.25 and 0.5 μM but some studies suggest that K D may have much greater range [60,61].
Without loss of generality, here we consider a basal Ca 2+ concentration of c 0 = 0.05 μM and set K D = 0.3 μM and R f = 100. Because we have no previous knowledge for the value of f m , we consider a range f m � 1−100 where the exact value depends on the ratio of the maximal intensity and the resting intensity. Using these values, we plot the corresponding Ca 2+ concentrations from the fluorescence data in [18] for various estimates of f m . Fig 13A-13C show the time traces of the converted fluorescence data for the impact of a 10 nl injection of Aβ at concentrations of a = 3 μg/ml, a = 10 μg/ml, and a = 30 μg/ml, respectively. In Fig 13A-13C each dashed plot corresponds to a different value of f m ranging from f m = 1 to f m = 100 (black) with n = 11 (f m = 1, 10, 20, . . ., 100). The maximum value is also highlighted for each conversion plot (circle) and provides the peak time for the three Aβ levels.
To study the impact of the conversion to Ca 2+ concentrations, Fig 14A-14C shows the corresponding maximum value of the concentration as a function of f m and the range of δf max between 6 and 20. These three dimensional plots allow us to better understand the impact of the conversion parameters on the maximum values of the fluorescence data in [18]. Again, because we do not have estimates for f m or c 0 , a true conversion from fluorescence to concentration is elusive. However, in all the profiles illustrated, each conversion does capture the changes in amplitude and latency to peak time observed experimentally as levels of Aβ are increased.