TY - JOUR
AU1 - Minicucci, Joseph
AU2 - Alfond, Molly
AU3 - Demuro, Angelo
AU4 - Gerberry, David
AU5 - Latulippe, Joe
AB - 1 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–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–13]. The ability of extracellular applied Aβ oligomers to induce cytosolic Ca2+ fluxes generated from both extracellular and intracellular sources has been shown using cultured mammalian cells [14–16]. We have subsequently characterized these two mechanisms as occurring by: i) formation of plasma membrane Ca2+ permeable pores [17], and ii) permeation of Aβ oligomers into the cytosol and inducing a PLC-dependent Ca2+ release from the ER [18]. As a critical secondary messenger, Ca2+ 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 Ca2+ homeostasis inducing mitochondrial stress [21, 22]. Because Aβ accumulation has been shown to alter intracellular Ca2+ levels, studying its impact on Ca2+ regulatory mechanisms is critical for better understanding the pathogenesis of AD. Intracellular Ca2+ regulation involves many distinct mechanisms working together. In the presence of Aβ, these Ca2+ regulatory mechanisms begin to fail [22, 23]. For example, the presence of Aβ has been shown to increase Ca2+ liberation from the ER through 1,4,5-Inositol-triphosphate receptors (IP3Rs) and ryanodine receptors (RyRs) [15, 24]. Aβ can also spontaneously form Ca2+-permeable pores in the plasma membrane [20, 25] creating uncontrolled influx of Ca2+ through the membrane. These alterations can cause stress on the ER that can further lead to dysregulation of Ca2+ in a feed-forward cyclical pattern [15, 22, 26, 27]. Such breakdowns in regulation can create aberrant or sustained elevated Ca2+ 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 IP3 signaling cascade and its consequential disruption of intracellular Ca2+ 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β42Os) into Xenopus oocytes triggered a PLC-dependent activation of IP3Rs in the ER membrane causing cytosolic Ca2+ 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 IP3 production and subsequent Ca2+ release through IP3Rs. We first build a computational model capable of tracking intracellular changes in Ca2+ concentration as a function of time. We assume that intracellular Aβ42Os (iAβ42Os) have a direct impact on G protein activation and PLC-mediated IP3 production. The experimental results in [18] provide data to calibrate our mathematical model and to test our modeling assumptions. We show that increasing iAβ42Os from small to large doses causes significant changes in the impact of Aβ on certain cellular mechanisms. Our model analysis substantiates that iAβ42Os have a widespread effect on IP3-mediated Ca2+ signaling. Because experimental recordings of Ca2+ signals are typically expressed as a ratio of fluorescence relative to the resting fluorescence before stimulation (Δf/f0), 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 IP3 production impacts Ca2+ 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 Ca2+ signals. We also explore how large doses of iAβ42Os alter the sensitivity of IP3Rs. 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β42Os, positively recognize by OC antibody and simply refer to them as Aβ for simplicity [6, 14]. 2 Methods 2.1 The closed-cell model development To investigate the impact of Aβ on Ca2+ regulation, we make use of the vast literature on calcium dynamics including the Ca2+ signaling “toolkit” [29–32]. We use experimental conditions and data from Xenopus oocytes to build a Ca2+ model using traditional methods of tracking the flux in and out of the cytoplasm. Let c denote the concentration of free Ca2+ ions in the cell cytoplasm, then the rate of change in intracellular Ca2+ 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 Ca2+ signaling in Xenopus oocytes is mostly due to IP3Rs as oocytes are deficient in RyRs. In an in vivo environment, both the Na+/Ca2+ exchanger and the plasma membrane Ca2+ ATPase pumps affect Ca2+ removal from the cytosol while receptor-operated Ca2+ channels lead to Ca2+ 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 Ca2+ to the Ca2+-dependent fluorescent dye. As such, the data extracted from our experiments intrinsically take into account the endogenous activity of the Na+/Ca2+ exchanger, the plasma membrane and SERCA Ca2+ ATPase pumps in the absence of specific blockers. Based on these conditions, we write where JIPR, JSERCA, JIN, and JPM are the fluxes due to IP3Rs, SERCA pump, a plasma membrane channel (such as a Receptor Operated Channel), and Plasma Membrane pump, respectively. The constant α is typically used to control the rate of transport of Ca2+ across the membrane to that across the ER. Let ce denote the concentration of ER calcium. With this, we assume a Ca2+ model of the form (1) (2) where γ is the ratio of cytoplasmic volume to ER volume. Note that we do not explicitly consider the effects of Ca2+ buffers. We assume that Ca2+ buffers are fast, immobile, and of low affinity (see [30, 32, 33] for further details on buffering). As such, Ca2+ 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 JIN and JPM are small compared to the contributions of the ER. As such, we set JIN − JPM ≈ 0 and reduce the model to a closed-cell model where Ca2+ transport only occurs between the ER and cytosol. Understanding that stable Ca2+ oscillations in Xenopus oocytes occur in the absence of external Ca2+ suggests that Ca2+ 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 Ca2+ regulation in general, and specifically in glial cells and neurons. For example, the direct exclusion of specific contributions from JIN and JPM may over-simplify Ca2+ 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 Ca2+ concentration even in the absence of extracellular Ca2+. Regardless, the simplified deterministic model does allow us to illustrate important dynamical properties of Ca2+ signaling patterns with minimal components. Accordingly, our closed-cell model assumes that Ca2+ flux into the cytosol is only due to the IP3R 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 Ca2+ flux as a mean-field approximation process that considers an average over a large number of IP3Rs. 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 IP3 signaling cascade. The flux terms in Eqs (1) and (2) can be modeled using various formulations, such as a saturating binding rate model for IP3R [34, 35] and Markov models [36–38]. For our purposes, we assume that the flux from IP3Rs follows a formulation based on previous models found in [39–41]. Thus, we write (3) where kf controls the density of IP3Rs, JER is the leak from the ER into the cytoplasm, and Po is the open probability of the IP3R. 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 Ca2+ free solution. As such, we assume a closed-cell environment with Ca2+ fluxes occurring only between the cytosol and the ER and set (4) where ct is the total number of moles in the cell divided by the cytoplasmic volume [32]. We then replace the term (ce − c) in Eq (3) with (γ(ct − c) − c). To model Po, we use the Li and Rinzel [42] simplification of the De Young and Keizer [39] formulation for the open probability of the IP3R (5) where y is the proportion of inactivated IP3Rs and p is the concentration of IP3 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 (6) (7) where Ki, for i = 1, …5 and k−4 and k−2 are parameters associated with the transition rates between various quasi-steady-states of the IP3R (see [32, 42], and [30] for details), and Vs and Ks 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 steady-states 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 Ca2+ solution showing periodic oscillations when p = 0.325 is shown in Fig 1C. Also illustrated there is the variable y in red. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Parameter values of the closed-cell Ca2+ base model. All IP3R parameters were adopted from De Young and Keizer (1992) while the Cellular and SERCA parameters were altered to match experimental results. https://doi.org/10.1371/journal.pone.0246116.t001 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Dynamics and bifurcation structure for constant IP3. A shows the bifurcation structure for Eqs (6) and (7) for constant values of p. A classic Hopf bubble emerges between two Hopf bifurcation points labeled HB1 and HB2. The black solid line corresponds to stable fixed points while the dashed black curve are unstable fixed points. B shows the nullclines when p = 0 (dashed curves) and when p = 0.325 (solid curves). The nullclines corresponding to dc/dt = 0 and dy/dt = 0 are given by the red and green curves, respectively. A trace of the solution when p = 0.325 is given by the blue trajectory and is attracted to a period orbit (dark blue). The oscillating solution to the model when p = 0.325 is shown in C as a function of time with the corresponding solution to the y equation is shown in red. https://doi.org/10.1371/journal.pone.0246116.g001 In [18], changes in Ca2+ concentration occur as a consequence of the intracellular injection of Aβ. As such, new IP3 is synthesized within the cell during the experimental procedure. To account for the IP3 dynamics, we use the hybrid formulation of Politi et al. [43]. Let p denote the concentration of IP3 present in the cytosol, then we write (8) where is the maximal rate of IP3 production and depends on agonist concentration, KPLC characterizes the sensitivity of PLC to Ca2+, τp = 1/(k3k+ k5p) represents the characteristic time of IP3 turnover where k3k is the maximum rate of 3-kinase and k5p is the dephosphorylation rate, Kip3 k is the half-activation constant for 3-kinase, and η = k3k/(k3k+ k5p). Both Kip3 k and η are used to tune the positive and negative feedback Ca2+ in the IP3 metabolism [30]. The term will depend on the amount of activated PLC available and we alter the model by writing (9) where PLC is the fraction of activated PLC complexes, and VPLC is the IP3 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 IP3 production and that G is the fraction of activated G-protein complexes and write (10) (11) where PLCtot and Gtot 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 ka, kb, kc, and kd 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]. 2.2 The effects of Aβ Although Aβ is clearly implicated in the disruption of intracellular Ca2+ 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 Ca2+ that carry on over time. This cytotoxicity may be due to stress caused by persistent Ca2+ release through IP3Rs. Of particular interest are the spatiotemporal patterns of fluorescence Ca2+ signals evoked by Aβ at dose of 1 μg/ml. Recordings in different oocytes showed that Aβ led to various Ca2+ 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 Ca2+ show that the amplitude of the Ca2+ 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 Ca2+ 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 Ca2+ 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. Illustrated in Fig 2 are two diagrams showing the model assumptions for the interaction of Aβ on the IP3 signaling cascade along with the key model components for “small” and “large” doses. The black arrows (solid and dashed) emanating from and going into Ca2+ illustrate the flow of Ca2+ along with feedback mechanisms. The two red arrows emerging from Aβ in Fig 2A show the location of the impact of “small” doses of Aβ within the model structure. The blue arrows emerging from Aβ in Fig 2B show the mechanisms impacted by “large” doses of Aβ. The assumptions for how Aβ alters the mechanisms illustrated in Fig 2 are based on the model’s ability to reproduce dose-dependent experimental results and are discussed in greater detail below. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Model structure and components. Modeling assumptions for the location of impact of Aβ on the production of IP3 with key Ca2+ signaling mechanisms included in the closed-cell model. The key model assumptions for how Aβ impacts the IP3 signaling cascade are illustrated as red arrows for “small” doses in A. The impacted cellular mechanisms for “large” doses of Aβ are highlighted by the blue arrows in B. https://doi.org/10.1371/journal.pone.0246116.g002 Our closed-cell model must be able to reproduce slow monotonic increases in Ca2+ 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 Ca2+, using several antagonists, we suggest that Aβ acts upstream of IP3Rs and hypothesize that Aβ stimulates IP3 production by PLC in a G-protein-dependent 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 (12) where VR is a scalar, KR is the Aβ concentration producing half activation. The term q represents the effects of a current injection at time t = t1 of Aβ at concentration a and has the form (13) where H is the Heaviside function and e−r(t−t1) represents the decay of Aβ over time. To match the timeframe of the experimental injections, we set t1 = 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 IP3 production to depend on Aβ as follows (14) where V0 accounts for PLC mediated IP3 production under normal conditions, VQ accounts for influence of Aβ on PLC-mediated IP3 production, and KQ is the dissociation constant. The exponent in VPLC 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 Ca2+ signals in coordination with the time to peak in our filtering process. Various alternative structures for VPLC were explored numerically but those structures were deemed insufficient for generating the experimental Ca2+ signaling patterns outlined in [18]. As such, we have assumed that the maximal rate of PLC mediated IP3 production takes the form of Eq (14), but more data is needed to determine whether this assumption actually captures how Aβ alters PLC-mediated IP3 production. Altogether, our closed cell model consists of five differential equations with Aβ input driving the system. In summary, (15) (16) (17) (18) (19) where the term q given in Eq (13) simulates the intracellular injection at time t = t1 of Aβ at concentration a. Base parameters for the IP3, 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Parameter values of the closed-cell Ca2+ model. https://doi.org/10.1371/journal.pone.0246116.t002 2.1 The closed-cell model development To investigate the impact of Aβ on Ca2+ regulation, we make use of the vast literature on calcium dynamics including the Ca2+ signaling “toolkit” [29–32]. We use experimental conditions and data from Xenopus oocytes to build a Ca2+ model using traditional methods of tracking the flux in and out of the cytoplasm. Let c denote the concentration of free Ca2+ ions in the cell cytoplasm, then the rate of change in intracellular Ca2+ 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 Ca2+ signaling in Xenopus oocytes is mostly due to IP3Rs as oocytes are deficient in RyRs. In an in vivo environment, both the Na+/Ca2+ exchanger and the plasma membrane Ca2+ ATPase pumps affect Ca2+ removal from the cytosol while receptor-operated Ca2+ channels lead to Ca2+ 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 Ca2+ to the Ca2+-dependent fluorescent dye. As such, the data extracted from our experiments intrinsically take into account the endogenous activity of the Na+/Ca2+ exchanger, the plasma membrane and SERCA Ca2+ ATPase pumps in the absence of specific blockers. Based on these conditions, we write where JIPR, JSERCA, JIN, and JPM are the fluxes due to IP3Rs, SERCA pump, a plasma membrane channel (such as a Receptor Operated Channel), and Plasma Membrane pump, respectively. The constant α is typically used to control the rate of transport of Ca2+ across the membrane to that across the ER. Let ce denote the concentration of ER calcium. With this, we assume a Ca2+ model of the form (1) (2) where γ is the ratio of cytoplasmic volume to ER volume. Note that we do not explicitly consider the effects of Ca2+ buffers. We assume that Ca2+ buffers are fast, immobile, and of low affinity (see [30, 32, 33] for further details on buffering). As such, Ca2+ 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 JIN and JPM are small compared to the contributions of the ER. As such, we set JIN − JPM ≈ 0 and reduce the model to a closed-cell model where Ca2+ transport only occurs between the ER and cytosol. Understanding that stable Ca2+ oscillations in Xenopus oocytes occur in the absence of external Ca2+ suggests that Ca2+ 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 Ca2+ regulation in general, and specifically in glial cells and neurons. For example, the direct exclusion of specific contributions from JIN and JPM may over-simplify Ca2+ 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 Ca2+ concentration even in the absence of extracellular Ca2+. Regardless, the simplified deterministic model does allow us to illustrate important dynamical properties of Ca2+ signaling patterns with minimal components. Accordingly, our closed-cell model assumes that Ca2+ flux into the cytosol is only due to the IP3R 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 Ca2+ flux as a mean-field approximation process that considers an average over a large number of IP3Rs. 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 IP3 signaling cascade. The flux terms in Eqs (1) and (2) can be modeled using various formulations, such as a saturating binding rate model for IP3R [34, 35] and Markov models [36–38]. For our purposes, we assume that the flux from IP3Rs follows a formulation based on previous models found in [39–41]. Thus, we write (3) where kf controls the density of IP3Rs, JER is the leak from the ER into the cytoplasm, and Po is the open probability of the IP3R. 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 Ca2+ free solution. As such, we assume a closed-cell environment with Ca2+ fluxes occurring only between the cytosol and the ER and set (4) where ct is the total number of moles in the cell divided by the cytoplasmic volume [32]. We then replace the term (ce − c) in Eq (3) with (γ(ct − c) − c). To model Po, we use the Li and Rinzel [42] simplification of the De Young and Keizer [39] formulation for the open probability of the IP3R (5) where y is the proportion of inactivated IP3Rs and p is the concentration of IP3 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 (6) (7) where Ki, for i = 1, …5 and k−4 and k−2 are parameters associated with the transition rates between various quasi-steady-states of the IP3R (see [32, 42], and [30] for details), and Vs and Ks 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 steady-states 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 Ca2+ solution showing periodic oscillations when p = 0.325 is shown in Fig 1C. Also illustrated there is the variable y in red. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Parameter values of the closed-cell Ca2+ base model. All IP3R parameters were adopted from De Young and Keizer (1992) while the Cellular and SERCA parameters were altered to match experimental results. https://doi.org/10.1371/journal.pone.0246116.t001 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Dynamics and bifurcation structure for constant IP3. A shows the bifurcation structure for Eqs (6) and (7) for constant values of p. A classic Hopf bubble emerges between two Hopf bifurcation points labeled HB1 and HB2. The black solid line corresponds to stable fixed points while the dashed black curve are unstable fixed points. B shows the nullclines when p = 0 (dashed curves) and when p = 0.325 (solid curves). The nullclines corresponding to dc/dt = 0 and dy/dt = 0 are given by the red and green curves, respectively. A trace of the solution when p = 0.325 is given by the blue trajectory and is attracted to a period orbit (dark blue). The oscillating solution to the model when p = 0.325 is shown in C as a function of time with the corresponding solution to the y equation is shown in red. https://doi.org/10.1371/journal.pone.0246116.g001 In [18], changes in Ca2+ concentration occur as a consequence of the intracellular injection of Aβ. As such, new IP3 is synthesized within the cell during the experimental procedure. To account for the IP3 dynamics, we use the hybrid formulation of Politi et al. [43]. Let p denote the concentration of IP3 present in the cytosol, then we write (8) where is the maximal rate of IP3 production and depends on agonist concentration, KPLC characterizes the sensitivity of PLC to Ca2+, τp = 1/(k3k+ k5p) represents the characteristic time of IP3 turnover where k3k is the maximum rate of 3-kinase and k5p is the dephosphorylation rate, Kip3 k is the half-activation constant for 3-kinase, and η = k3k/(k3k+ k5p). Both Kip3 k and η are used to tune the positive and negative feedback Ca2+ in the IP3 metabolism [30]. The term will depend on the amount of activated PLC available and we alter the model by writing (9) where PLC is the fraction of activated PLC complexes, and VPLC is the IP3 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 IP3 production and that G is the fraction of activated G-protein complexes and write (10) (11) where PLCtot and Gtot 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 ka, kb, kc, and kd 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]. 2.2 The effects of Aβ Although Aβ is clearly implicated in the disruption of intracellular Ca2+ 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 Ca2+ that carry on over time. This cytotoxicity may be due to stress caused by persistent Ca2+ release through IP3Rs. Of particular interest are the spatiotemporal patterns of fluorescence Ca2+ signals evoked by Aβ at dose of 1 μg/ml. Recordings in different oocytes showed that Aβ led to various Ca2+ 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 Ca2+ show that the amplitude of the Ca2+ 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 Ca2+ 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 Ca2+ 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. Illustrated in Fig 2 are two diagrams showing the model assumptions for the interaction of Aβ on the IP3 signaling cascade along with the key model components for “small” and “large” doses. The black arrows (solid and dashed) emanating from and going into Ca2+ illustrate the flow of Ca2+ along with feedback mechanisms. The two red arrows emerging from Aβ in Fig 2A show the location of the impact of “small” doses of Aβ within the model structure. The blue arrows emerging from Aβ in Fig 2B show the mechanisms impacted by “large” doses of Aβ. The assumptions for how Aβ alters the mechanisms illustrated in Fig 2 are based on the model’s ability to reproduce dose-dependent experimental results and are discussed in greater detail below. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Model structure and components. Modeling assumptions for the location of impact of Aβ on the production of IP3 with key Ca2+ signaling mechanisms included in the closed-cell model. The key model assumptions for how Aβ impacts the IP3 signaling cascade are illustrated as red arrows for “small” doses in A. The impacted cellular mechanisms for “large” doses of Aβ are highlighted by the blue arrows in B. https://doi.org/10.1371/journal.pone.0246116.g002 Our closed-cell model must be able to reproduce slow monotonic increases in Ca2+ 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 Ca2+, using several antagonists, we suggest that Aβ acts upstream of IP3Rs and hypothesize that Aβ stimulates IP3 production by PLC in a G-protein-dependent 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 (12) where VR is a scalar, KR is the Aβ concentration producing half activation. The term q represents the effects of a current injection at time t = t1 of Aβ at concentration a and has the form (13) where H is the Heaviside function and e−r(t−t1) represents the decay of Aβ over time. To match the timeframe of the experimental injections, we set t1 = 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 IP3 production to depend on Aβ as follows (14) where V0 accounts for PLC mediated IP3 production under normal conditions, VQ accounts for influence of Aβ on PLC-mediated IP3 production, and KQ is the dissociation constant. The exponent in VPLC 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 Ca2+ signals in coordination with the time to peak in our filtering process. Various alternative structures for VPLC were explored numerically but those structures were deemed insufficient for generating the experimental Ca2+ signaling patterns outlined in [18]. As such, we have assumed that the maximal rate of PLC mediated IP3 production takes the form of Eq (14), but more data is needed to determine whether this assumption actually captures how Aβ alters PLC-mediated IP3 production. Altogether, our closed cell model consists of five differential equations with Aβ input driving the system. In summary, (15) (16) (17) (18) (19) where the term q given in Eq (13) simulates the intracellular injection at time t = t1 of Aβ at concentration a. Base parameters for the IP3, 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Parameter values of the closed-cell Ca2+ model. https://doi.org/10.1371/journal.pone.0246116.t002 3 Model results 3.1 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 IP3R parameters may also vary, but for now we simply focus on the SERCA parameter Ks. Illustrated in Fig 3 are various scaled model solutions along with a partial bifurcation diagram highlighting the key behaviors of the model when Ks is altered. Model solutions illustrated in Fig 3 have been scaled according to the following (20) where KD = 0.3 is a dissociation constant that depends on indicator properties, and c0 is the resting Ca2+ concentration. In addition, to set the initial condition c0 (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 Fig 3 are slightly different as altering Ks also changes the Ca2+ homeostasis level in the model but typically range between (0.01, 0.15). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Model solutions mimic experimental Ca2+ patterns for doses of a = 1 μg/ml of Aβ. The dependence of model solutions for a dose of 1 μg/ml of Aβ on the cellular parameter Ks is investigated. A shows an increasing Ca2+ signal that settles to an increased steady-state when Ks = 0.15. B shows that oscillations in Ca2+ can exhibit increasing amplitudes such as those found in Fig 1D in [18]. C and D show that as the value of Ks decreases, the oscillatory patterns of the model reproduce the spike-like Ca2+ signals observed in Fig 1E in [18]. E illustrates an oscillatory solution with an increased steady-state Ca2+ homeostasis level when Ks is just above the Hopf point. Both D and E show the traces for cs, y, and p in blue, red, and black, respectively. F shows a simplification of the scaled bifurcation diagram with the bifurcation parameter Ks. Notice that as Ks 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 Ks = 0.09 include multiple Hopf bifurcations and has more complex structure than what is presented here. https://doi.org/10.1371/journal.pone.0246116.g003 Fig 3A–3C show responses similar to those in Fig 1 of [18]. More specifically, Fig 3A shows a solution where Ca2+ increases to a new steady-state when Ks = 0.15. Fig 3B illustrates a solution that has increasing amplitude oscillations when Ks = 0.12 while Fig 3C shows repetitive oscillations when Ks = 0.118. In both Fig 3B and 3C, model Ca2+ signals occur between 2-5 minutes, matching the experimental timescale for these types of responses. The responses in Fig 3D–3E show spike-like Ca2+ pattern that have a smaller frequency when Ks = 0.11 and a decreasing amplitude oscillatory solution when Ks = 0.125, respectively. A partial bifurcation diagram with Ks as the bifurcation parameter is provided in Fig 3F. As the parameter Ks decreases from Ks = 0.15, a transition from stable fixed points into periodic orbits occurs through a Hopf bifurcation point around HB3 ≈ 0.1242. As Ks 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 Ks = 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. Ks is the dissociation constant for the SERCA pump and is the Ca2+ concentration occupying half of the binding sites of the pump. A smaller Ks value corresponds to needing less Ca2+ to attain 50% of the maximal response for the pump. Whether changes in Ks 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 Ca2+ signals and not because of Aβ’s direct impact on the SERCA pump. However, more analysis is needed to fully understand how different doses of Aβ may influence the generation of various Ca2+ signals. When looking at the model, the subsystem given by Eqs (18) and (19) is driving IP3 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 Ca2+ using PLC as a bifurcation parameter for the subsystem given in Eqs (15)–(17). Notice that the dynamics of Ca2+ 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 Ca2+ oscillations. Fig 4B shows the subsystem solution when PLC = 0.005 for the base parameters given in Table 2 for small doses. Intracellular Ca2+ signal, the proportion of inactivated IP3Rs, and the concentration of IP3 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Dynamics of Ca2+ using PLC as a parameter for doses of a ≤ 1 μg/ml of Aβ. Model dynamics for the subsystem Eqs (15)–(17) in terms of PLC and Ks. The bifurcation diagram when Ks = 0.15 with the subsystem parameters given in Table 2, is shown in A. The figure shows a typical Hopf bubble between two Hopf bifurcation values labeled HB4 and HB5 as PLC is varied. The subsystem solution when PLC = 0.005 is presented in B, where c, y, and p are shown in the blue, red, and black traces, respectively. C shows the two parameter bifurcation diagram when both PLC and a are varied. Note that only the small doses of a are considered and the region of oscillations shifts to the right as a decreases and PLC increases. D shows the subsystem bifurcation diagram when Ks = 0.011 including two Hopf bifurcation values labeled HB6 and HB7. E shows the subsystem solution when PLC = 0.01 where c, y, and p are shown in the blue, red, and black traces, respectively. F shows the two parameter diagram tracking the location of the Hopf bifurcation points when PLC and Ks are treated as parameters. The parameter space is separated into regions where periodic orbits exist and don’t. The red cross and triangle correspond to the location of the parameter values used to generate the solutions in B and E, respectively. https://doi.org/10.1371/journal.pone.0246116.g004 In Fig 3, we showed that the model solutions will behave differently as the cellular parameter Ks is varied. Here, we also look at the impact of varying Ks on the solutions of the subsystem Eqs (15)–(17). Illustrated in Fig 4D is the bifurcation diagram when Ks = 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 Ca2+ 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 Ks = 0.11. To better understand the impact of changes in both PLC and Ks 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 Ks 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 Ks 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 the phase space solutions for a = 0.1 (black), a = 0.5 (magenta), and a = 1 (blue). The dashed lines in this figure correspond to the nullclines for Eq (18) (red) and Eq (19) (respective color). The new a-dependent steady-state values occur at the intersection of the respective dashed lines for each a value. All three solutions shown in Fig 5C start at the a = 0 equilibrium value of (PLC0, G0) ≈ (2.415 × 10−4, 1.518 × 10−3). When comparing the analysis shown in Fig 4, the values of PLC produced through the subsystem given by Eqs (18) and (19) will generate oscillatory responses when Ks is decreased from Ks = 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Steady-state values for PLC and G for doses of a ≤ 1 μg/ml of Aβ. The steady-state fraction of activated PLC and G-proteins settles to a new value when a = 0.1 (black), a = 0.5 (magenta), and a = 1 (blue) in A and B, respectively. Model solutions for the various small doses a-values are shown on the phase plane for PLC and G in C. The dashed lines correspond to the PLC nullcline (red) and the G nullclines (black, magenta, and blue). https://doi.org/10.1371/journal.pone.0246116.g005 3.2 Dose response relationship between amplitude and latency Dose-response experiments in Xenopus oocytes demonstrate two major effects on Ca2+ fluxes following increasing doses of Aβ: the amplitude of the Ca2+ 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 Ca2+ 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Amplitude and latency of model solutions vary with doses of Aβ. For the small doses parameters, A shows the model captures the increase in Ca2+ signal amplitude as well as the decrease in time to peak onset. B shows the long term impact of Aβ for doses of Aβ corresponding to a = 1, a = 3, a = 10, and a = 30 in black, blue, red, and green, respectively. The dashed lines in B correspond to the steady-state values in the event where the amount of Aβ does not decay and is fixed. C shows the location of 100 stochastically chosen cells under the given a-value. Each color-coded circle corresponds to the location of the solution peak and time of peak when cellular parameters are varied uniformly with 10% variation for the particular a value. The dashed black curve corresponds to the location of the amplitude peak for the small doses parameters for a ranging between (0.1, 40). https://doi.org/10.1371/journal.pone.0246116.g006 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 Ca2+ signaling pathway responsible for cognition [13, 47, 48]. As such, a slow accumulation of Aβ may increase the cytosolic Ca2+ 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 kf, JER, 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 small- and large-doses parameter sets. 3.1 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 IP3R parameters may also vary, but for now we simply focus on the SERCA parameter Ks. Illustrated in Fig 3 are various scaled model solutions along with a partial bifurcation diagram highlighting the key behaviors of the model when Ks is altered. Model solutions illustrated in Fig 3 have been scaled according to the following (20) where KD = 0.3 is a dissociation constant that depends on indicator properties, and c0 is the resting Ca2+ concentration. In addition, to set the initial condition c0 (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 Fig 3 are slightly different as altering Ks also changes the Ca2+ homeostasis level in the model but typically range between (0.01, 0.15). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Model solutions mimic experimental Ca2+ patterns for doses of a = 1 μg/ml of Aβ. The dependence of model solutions for a dose of 1 μg/ml of Aβ on the cellular parameter Ks is investigated. A shows an increasing Ca2+ signal that settles to an increased steady-state when Ks = 0.15. B shows that oscillations in Ca2+ can exhibit increasing amplitudes such as those found in Fig 1D in [18]. C and D show that as the value of Ks decreases, the oscillatory patterns of the model reproduce the spike-like Ca2+ signals observed in Fig 1E in [18]. E illustrates an oscillatory solution with an increased steady-state Ca2+ homeostasis level when Ks is just above the Hopf point. Both D and E show the traces for cs, y, and p in blue, red, and black, respectively. F shows a simplification of the scaled bifurcation diagram with the bifurcation parameter Ks. Notice that as Ks 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 Ks = 0.09 include multiple Hopf bifurcations and has more complex structure than what is presented here. https://doi.org/10.1371/journal.pone.0246116.g003 Fig 3A–3C show responses similar to those in Fig 1 of [18]. More specifically, Fig 3A shows a solution where Ca2+ increases to a new steady-state when Ks = 0.15. Fig 3B illustrates a solution that has increasing amplitude oscillations when Ks = 0.12 while Fig 3C shows repetitive oscillations when Ks = 0.118. In both Fig 3B and 3C, model Ca2+ signals occur between 2-5 minutes, matching the experimental timescale for these types of responses. The responses in Fig 3D–3E show spike-like Ca2+ pattern that have a smaller frequency when Ks = 0.11 and a decreasing amplitude oscillatory solution when Ks = 0.125, respectively. A partial bifurcation diagram with Ks as the bifurcation parameter is provided in Fig 3F. As the parameter Ks decreases from Ks = 0.15, a transition from stable fixed points into periodic orbits occurs through a Hopf bifurcation point around HB3 ≈ 0.1242. As Ks 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 Ks = 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. Ks is the dissociation constant for the SERCA pump and is the Ca2+ concentration occupying half of the binding sites of the pump. A smaller Ks value corresponds to needing less Ca2+ to attain 50% of the maximal response for the pump. Whether changes in Ks 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 Ca2+ signals and not because of Aβ’s direct impact on the SERCA pump. However, more analysis is needed to fully understand how different doses of Aβ may influence the generation of various Ca2+ signals. When looking at the model, the subsystem given by Eqs (18) and (19) is driving IP3 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 Ca2+ using PLC as a bifurcation parameter for the subsystem given in Eqs (15)–(17). Notice that the dynamics of Ca2+ 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 Ca2+ oscillations. Fig 4B shows the subsystem solution when PLC = 0.005 for the base parameters given in Table 2 for small doses. Intracellular Ca2+ signal, the proportion of inactivated IP3Rs, and the concentration of IP3 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Dynamics of Ca2+ using PLC as a parameter for doses of a ≤ 1 μg/ml of Aβ. Model dynamics for the subsystem Eqs (15)–(17) in terms of PLC and Ks. The bifurcation diagram when Ks = 0.15 with the subsystem parameters given in Table 2, is shown in A. The figure shows a typical Hopf bubble between two Hopf bifurcation values labeled HB4 and HB5 as PLC is varied. The subsystem solution when PLC = 0.005 is presented in B, where c, y, and p are shown in the blue, red, and black traces, respectively. C shows the two parameter bifurcation diagram when both PLC and a are varied. Note that only the small doses of a are considered and the region of oscillations shifts to the right as a decreases and PLC increases. D shows the subsystem bifurcation diagram when Ks = 0.011 including two Hopf bifurcation values labeled HB6 and HB7. E shows the subsystem solution when PLC = 0.01 where c, y, and p are shown in the blue, red, and black traces, respectively. F shows the two parameter diagram tracking the location of the Hopf bifurcation points when PLC and Ks are treated as parameters. The parameter space is separated into regions where periodic orbits exist and don’t. The red cross and triangle correspond to the location of the parameter values used to generate the solutions in B and E, respectively. https://doi.org/10.1371/journal.pone.0246116.g004 In Fig 3, we showed that the model solutions will behave differently as the cellular parameter Ks is varied. Here, we also look at the impact of varying Ks on the solutions of the subsystem Eqs (15)–(17). Illustrated in Fig 4D is the bifurcation diagram when Ks = 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 Ca2+ 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 Ks = 0.11. To better understand the impact of changes in both PLC and Ks 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 Ks 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 Ks 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 the phase space solutions for a = 0.1 (black), a = 0.5 (magenta), and a = 1 (blue). The dashed lines in this figure correspond to the nullclines for Eq (18) (red) and Eq (19) (respective color). The new a-dependent steady-state values occur at the intersection of the respective dashed lines for each a value. All three solutions shown in Fig 5C start at the a = 0 equilibrium value of (PLC0, G0) ≈ (2.415 × 10−4, 1.518 × 10−3). When comparing the analysis shown in Fig 4, the values of PLC produced through the subsystem given by Eqs (18) and (19) will generate oscillatory responses when Ks is decreased from Ks = 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Steady-state values for PLC and G for doses of a ≤ 1 μg/ml of Aβ. The steady-state fraction of activated PLC and G-proteins settles to a new value when a = 0.1 (black), a = 0.5 (magenta), and a = 1 (blue) in A and B, respectively. Model solutions for the various small doses a-values are shown on the phase plane for PLC and G in C. The dashed lines correspond to the PLC nullcline (red) and the G nullclines (black, magenta, and blue). https://doi.org/10.1371/journal.pone.0246116.g005 3.2 Dose response relationship between amplitude and latency Dose-response experiments in Xenopus oocytes demonstrate two major effects on Ca2+ fluxes following increasing doses of Aβ: the amplitude of the Ca2+ 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 Ca2+ 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Amplitude and latency of model solutions vary with doses of Aβ. For the small doses parameters, A shows the model captures the increase in Ca2+ signal amplitude as well as the decrease in time to peak onset. B shows the long term impact of Aβ for doses of Aβ corresponding to a = 1, a = 3, a = 10, and a = 30 in black, blue, red, and green, respectively. The dashed lines in B correspond to the steady-state values in the event where the amount of Aβ does not decay and is fixed. C shows the location of 100 stochastically chosen cells under the given a-value. Each color-coded circle corresponds to the location of the solution peak and time of peak when cellular parameters are varied uniformly with 10% variation for the particular a value. The dashed black curve corresponds to the location of the amplitude peak for the small doses parameters for a ranging between (0.1, 40). https://doi.org/10.1371/journal.pone.0246116.g006 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 Ca2+ signaling pathway responsible for cognition [13, 47, 48]. As such, a slow accumulation of Aβ may increase the cytosolic Ca2+ 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 kf, JER, 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 small- and large-doses parameter sets. 4 Large doses parameter fitting The model developed in the previous section tracks Ca2+ concentration as a function of time. The experimental data in [18] tracks changes in Ca2+ as a ratio of changes in fluorescence intensity with baseline fluorescence levels. This is often written as δf = (f−f0)/f0 = Δf/f0 with f0 representing the fluorescence intensity at resting Ca2+ concentration. To better understand the impact of Aβ on Ca2+ dynamics through modeling, we first rescale fluorescence measurements to Ca2+ concentrations. According to Maravall et al. [28], changes in Ca2+ concentration are associated with changes in fluorescence through the equation (21) where fmax is the intensity of the dye at maximum Ca2+ concentration, Rf = fmax/fmin is the indicator’s dynamic range with fmin being the intensity at minimum Ca2+ concentration, δfmax is the saturation of the Ca2+ indicator, and fm = fmax/f0. 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 Ca2+ concentrations and link model solutions with experimental data. We first fix the scaling parameters KD = 0.3, Rf = 100, fm = 40 and then determine the value of the model parameters that will evoke the appropriate Ca2+ 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 qs from an admissible parameter space (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 (22) where sed(i) is the scaled experimental data value at i, and ics(i, q) is the corresponding (interpolated) scaled Ca2+ 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 cs (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 IP3R 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 IP3R 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 IP3R parameters are varied. This simulation is included to clearly illustrate the need for altering all model parameters, particularly the IP3R 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Variation in all parameter subsets required to reproduce the impact of Aβ for large doses. A shows a “best” fit solution when the Cellular, SERCA, and the IP3 Receptor parameters are kept fixed. Notice that altering the remaining parameters (IP3 Model, PLC, and G-Protein) cannot capture the observed Ca2+ signal. B shows a “best” fit solution when only the IP3R parameters are kept fixed. Notice that without altering the IP3R parameters, model solution peaks and decay also do not reproduce the observed experimental behaviors for large doses of Aβ. https://doi.org/10.1371/journal.pone.0246116.g007 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 Ca2+ concentration c illustrating that the scaling procedure does not effect the model’s ability to capture the general behavior of the Ca2+ 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 IP3Rs that are inactivated by Ca2+ remains fairly high over the course of hours acting to suppress Ca2+ spikes over time. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Model matches experimental data for large doses of Aβ. Model solutions for the “large” dose parameter set are illustrated in this figure. A shows model simulation (smooth curve) when a = 3 (blue), a = 10 (red), and a = 30 (green) overlaid on top of the rescaled experimental data (dashed curve) of [18]. Note that model solutions for the “large” dose parameters when a = 1 (black) is also shown here. B and C show the time evolution of model IP3 and PLC, respectively. D shows the unscaled Ca2+ concentration given by the model variable c for the three levels of Aβ and for a = 1 (black). E and F show the time evolution of p and y on the order of hours, respectively. Due to the Aβ decay incorporated in the model, all model solutions will eventually go back to the steady-state values. https://doi.org/10.1371/journal.pone.0246116.g008 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 Ca2+ 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 Ca2+ 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. 4.1 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] Ca2+ 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. Fig 9 shows the mean (solid curves) bounded within one standard deviation (dashed curves) for simulations around the concentration values Aβ = 3, 10, and 30, respectively. Again, we also include the result for Aβ = 1 for comparison and as a lower bound for the large doses range. The results illustrated in Fig 9 show that the model solutions are stable under parameter variation and continue to capture both the changes in amplitude and the peak time. Even if the large doses parameter set given 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 PLC-mediated IP3 production. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Model solutions under variation of parameters. The mean and corresponding standard deviations when the model is simulated for n = 100, 000 stochastically chosen parameter sets. The solid curves corresponds to the mean response and the dashed curves are the standard deviation above and below the mean. A and B illustrate the uncertainty in solutions when parameters are selected from a set with 5% and 10% deviation from the large doses base values given in Table 2, respectively. https://doi.org/10.1371/journal.pone.0246116.g009 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 Ca2+ dynamics [50]. To do this, we characterize the resulting Ca2+ dynamics with two quantities: the peak Ca2+ 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 Ca2+ concentration while Table 4 shows the correlation of the time of peak. The tables list the 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 Ca2+ solutions are the parameters kd, kb, and K1 while the parameters that most increase the amplitude are VQ, kc and ka, as a is increased. Similarly, the parameters that most decrease the time of peak of Ca2+ solutions are the parameters VQ, ka, and kc while the parameters that most increase the time peak are K1, kd and kb. 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the maximum Ca2+ 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. https://doi.org/10.1371/journal.pone.0246116.t003 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the time to peak when maximum Ca2+ 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. https://doi.org/10.1371/journal.pone.0246116.t004 When looking at the PRCC analysis, it appears that the PLC and G-protein rate constants ka, kb, kc, and kd all have a large impact on the solution patterns in terms of solution peak and time to peak. Recall that, ka and kc are the activation rates for PLC and G-proteins, respectively. As the activation rates increase, this will lead to an increase in IP3 production and you will see the peak of the Ca2+ signal occur sooner. On the other hand, kb and kd correspond to the inactivation of PLC and G-proteins, respectively. A higher inactivation rate for both PLC and G-proteins will decrease IP3 production and thus lower the peak amplitude of Ca2+ responses. From a biological perspective this makes sense, once PLC is activated, the production of IP3 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 Ca2+ 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 Ca2+ responses as fewer IP3 are available for binding to the IP3R. 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 IP3 available for IP3R binding. Further analysis on the impact of these parameters is provided below. As noted above, the parameter K1 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 IP3 to one of the IP3R model subunits when no inactivating Ca2+ is present. As such, this parameter helps drive the IP3R dynamics. In the model, an increase in K1 has an inactivating effect on the IP3R since either the unbinding rate of IP3 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 K1 is critical for understanding the Ca2+ patterns of the model. Because of the influence of this parameter on model solutions, this suggests that the IP3R dynamics does contribute to the observed Ca2+ patterns. We analyze the model below to further expand on the influence of K1 on model solutions. As the model suggests that changes in K1 may be dependent on Aβ levels, further investigations on the connection between the IP3R and Aβ may be necessary. 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 Ks appears to be tied with the size of a. More specifically, the PRCC for Ks when a = 3, 10, and 30 are 0.712, 0.514, and −0.147, respectively. This implies that as a increases, altering Ks 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 KPLC and k5p 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. 4.2 Impact of Aβ on IP3R for large doses The impact of Aβ on the IP3 signaling cascade appears to be concentration dependent. Not surprising, the PRCC analysis suggests that the rates ka, kb, kc, and kd 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 VQ. Recall that VQ controls the influence of Aβ on PLC-mediated IP3 production. As such, it is no surprise that VQ also plays a significant role in the solution patterns. Fig 10 shows the impact of altering VQ on model solutions for a = 3, 10, and 30, in A, B, and C, respectively. As the parameter VQ increases from the small doses value of VQ = 7.82 to the large doses value of VQ = 380 we see that the model solutions shift up and to the left. This is highlighted by the curved arrow in each figure. The large doses value of VQ = 380 has been singled out using the solid black trace while 9 other solutions (with various VQ values) are shown as dashed colored traces. The solution for the values VQ = 30, 380, and 480 have been highlighted in each figure for reference. The results of Fig 10 also confirm the PRCC analysis that VQ is positively correlated with the peak amplitude and negatively correlated with the peak time onset. Clearly, altering VQ impacts both the solution amplitude and the time to peak. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. PRCC prediction on solution amplitude and time to peak for model parameter VQ. The impact of VQ is shown as a series of curves for a = 3, a = 10, and a = 30 in A, B, and C, respectively. In each diagram, the curved black arrow tracks the shift in the peak of solutions as VQ takes on various values ranging from VQ = 30 to VQ = 480. The black trace in each diagram represents the baseline VQ value for the large doses parameter set. https://doi.org/10.1371/journal.pone.0246116.g010 Although the PRCC identifies the parameter K1 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 Ca2+ signaling patterns observed experimentally. Varying k−4 has a direct impact on the Ca2+ signal tail and partly controls the decay of the signals, but does not alter the amplitude or time to peak significantly. Both K1 and k−4 are parameters that help control the dynamics of IP3Rs. Shown in Fig 11 are two diagrams that illustrate the impact of Aβ on the IP3R itself through the parameters K1 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 K1 model solutions. Starting with the large doses parameters, we simulate the model by altering K1 from the base small doses value of K1 = 0.13 (bold black trace) and increasing the parameter to the large doses value K1 = 0.13 (smooth red trace). As is suggested by the PRCC analysis, we see that K1 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 Ca2+ 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 IP3R parameters will have various effects on solutions similar to the impact of varying K1. Changes to IP3R parameters seem necessary in order to capture the increasingly rapid decay and suggests that Aβ for large doses may act to desensitize the IP3R. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Impact of large doses of Aβ on IP3R dynamics. The impact of the IP3R model parameters K1 and k−4 are shown in A and B, respectively. The traces shown use the large doses parameters except for the values highlighted in each diagram. The top bold black traces correspond to the model solution when the parameter value for the small doses is used. The red traces are the model solutions for the parameter values of the large doses. The black traces between the bold and red correspond to intermediate parameter values as given in each diagram. https://doi.org/10.1371/journal.pone.0246116.g011 Whether Aβ directly interferes with IP3Rs 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 IP3Rs is affected by Aβ. Of particular interest is the role of the IP3R parameters in capturing the observed rapid decay of Ca2+ signal for large Aβ doses. 4.3 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 IP3 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 Ca2+ 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 Ca2+ 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 Ca2+ regulatory mechanisms will be necessary to describe whole-cell calcium dynamics in a biologically robust way. Our model construction assumes that iAβ42Os (1) act as an agonist for G-protein activation, and (2) affect the maximal rate of PLC mediated IP3 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β42Os 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 Ca2+ 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 VQ ≠ 0. That led us to incorporate the Aβ-dependent term for the maximal rate of PLC mediated IP3 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 IP3 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β42Os. 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 IP3 production cascade. 4.1 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] Ca2+ 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. Fig 9 shows the mean (solid curves) bounded within one standard deviation (dashed curves) for simulations around the concentration values Aβ = 3, 10, and 30, respectively. Again, we also include the result for Aβ = 1 for comparison and as a lower bound for the large doses range. The results illustrated in Fig 9 show that the model solutions are stable under parameter variation and continue to capture both the changes in amplitude and the peak time. Even if the large doses parameter set given 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 PLC-mediated IP3 production. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Model solutions under variation of parameters. The mean and corresponding standard deviations when the model is simulated for n = 100, 000 stochastically chosen parameter sets. The solid curves corresponds to the mean response and the dashed curves are the standard deviation above and below the mean. A and B illustrate the uncertainty in solutions when parameters are selected from a set with 5% and 10% deviation from the large doses base values given in Table 2, respectively. https://doi.org/10.1371/journal.pone.0246116.g009 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 Ca2+ dynamics [50]. To do this, we characterize the resulting Ca2+ dynamics with two quantities: the peak Ca2+ 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 Ca2+ concentration while Table 4 shows the correlation of the time of peak. The tables list the 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 Ca2+ solutions are the parameters kd, kb, and K1 while the parameters that most increase the amplitude are VQ, kc and ka, as a is increased. Similarly, the parameters that most decrease the time of peak of Ca2+ solutions are the parameters VQ, ka, and kc while the parameters that most increase the time peak are K1, kd and kb. 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the maximum Ca2+ 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. https://doi.org/10.1371/journal.pone.0246116.t003 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Partial rank correlation coefficient sensitivity analysis between model parameters (n = 100, 000) and the time to peak when maximum Ca2+ 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. https://doi.org/10.1371/journal.pone.0246116.t004 When looking at the PRCC analysis, it appears that the PLC and G-protein rate constants ka, kb, kc, and kd all have a large impact on the solution patterns in terms of solution peak and time to peak. Recall that, ka and kc are the activation rates for PLC and G-proteins, respectively. As the activation rates increase, this will lead to an increase in IP3 production and you will see the peak of the Ca2+ signal occur sooner. On the other hand, kb and kd correspond to the inactivation of PLC and G-proteins, respectively. A higher inactivation rate for both PLC and G-proteins will decrease IP3 production and thus lower the peak amplitude of Ca2+ responses. From a biological perspective this makes sense, once PLC is activated, the production of IP3 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 Ca2+ 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 Ca2+ responses as fewer IP3 are available for binding to the IP3R. 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 IP3 available for IP3R binding. Further analysis on the impact of these parameters is provided below. As noted above, the parameter K1 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 IP3 to one of the IP3R model subunits when no inactivating Ca2+ is present. As such, this parameter helps drive the IP3R dynamics. In the model, an increase in K1 has an inactivating effect on the IP3R since either the unbinding rate of IP3 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 K1 is critical for understanding the Ca2+ patterns of the model. Because of the influence of this parameter on model solutions, this suggests that the IP3R dynamics does contribute to the observed Ca2+ patterns. We analyze the model below to further expand on the influence of K1 on model solutions. As the model suggests that changes in K1 may be dependent on Aβ levels, further investigations on the connection between the IP3R and Aβ may be necessary. 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 Ks appears to be tied with the size of a. More specifically, the PRCC for Ks when a = 3, 10, and 30 are 0.712, 0.514, and −0.147, respectively. This implies that as a increases, altering Ks 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 KPLC and k5p 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. 4.2 Impact of Aβ on IP3R for large doses The impact of Aβ on the IP3 signaling cascade appears to be concentration dependent. Not surprising, the PRCC analysis suggests that the rates ka, kb, kc, and kd 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 VQ. Recall that VQ controls the influence of Aβ on PLC-mediated IP3 production. As such, it is no surprise that VQ also plays a significant role in the solution patterns. Fig 10 shows the impact of altering VQ on model solutions for a = 3, 10, and 30, in A, B, and C, respectively. As the parameter VQ increases from the small doses value of VQ = 7.82 to the large doses value of VQ = 380 we see that the model solutions shift up and to the left. This is highlighted by the curved arrow in each figure. The large doses value of VQ = 380 has been singled out using the solid black trace while 9 other solutions (with various VQ values) are shown as dashed colored traces. The solution for the values VQ = 30, 380, and 480 have been highlighted in each figure for reference. The results of Fig 10 also confirm the PRCC analysis that VQ is positively correlated with the peak amplitude and negatively correlated with the peak time onset. Clearly, altering VQ impacts both the solution amplitude and the time to peak. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. PRCC prediction on solution amplitude and time to peak for model parameter VQ. The impact of VQ is shown as a series of curves for a = 3, a = 10, and a = 30 in A, B, and C, respectively. In each diagram, the curved black arrow tracks the shift in the peak of solutions as VQ takes on various values ranging from VQ = 30 to VQ = 480. The black trace in each diagram represents the baseline VQ value for the large doses parameter set. https://doi.org/10.1371/journal.pone.0246116.g010 Although the PRCC identifies the parameter K1 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 Ca2+ signaling patterns observed experimentally. Varying k−4 has a direct impact on the Ca2+ signal tail and partly controls the decay of the signals, but does not alter the amplitude or time to peak significantly. Both K1 and k−4 are parameters that help control the dynamics of IP3Rs. Shown in Fig 11 are two diagrams that illustrate the impact of Aβ on the IP3R itself through the parameters K1 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 K1 model solutions. Starting with the large doses parameters, we simulate the model by altering K1 from the base small doses value of K1 = 0.13 (bold black trace) and increasing the parameter to the large doses value K1 = 0.13 (smooth red trace). As is suggested by the PRCC analysis, we see that K1 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 Ca2+ 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 IP3R parameters will have various effects on solutions similar to the impact of varying K1. Changes to IP3R parameters seem necessary in order to capture the increasingly rapid decay and suggests that Aβ for large doses may act to desensitize the IP3R. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Impact of large doses of Aβ on IP3R dynamics. The impact of the IP3R model parameters K1 and k−4 are shown in A and B, respectively. The traces shown use the large doses parameters except for the values highlighted in each diagram. The top bold black traces correspond to the model solution when the parameter value for the small doses is used. The red traces are the model solutions for the parameter values of the large doses. The black traces between the bold and red correspond to intermediate parameter values as given in each diagram. https://doi.org/10.1371/journal.pone.0246116.g011 Whether Aβ directly interferes with IP3Rs 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 IP3Rs is affected by Aβ. Of particular interest is the role of the IP3R parameters in capturing the observed rapid decay of Ca2+ signal for large Aβ doses. 4.3 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 IP3 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 Ca2+ 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 Ca2+ 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 Ca2+ regulatory mechanisms will be necessary to describe whole-cell calcium dynamics in a biologically robust way. Our model construction assumes that iAβ42Os (1) act as an agonist for G-protein activation, and (2) affect the maximal rate of PLC mediated IP3 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β42Os 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 Ca2+ 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 VQ ≠ 0. That led us to incorporate the Aβ-dependent term for the maximal rate of PLC mediated IP3 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 IP3 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β42Os. 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 IP3 production cascade. 5 Discussion Ca2+ 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 Ca2+ 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, IP3-mediated elementary Ca2+ signals, also referred to as puffs, are the building blocks of cellular Ca2+ signaling, and arise through the concerted opening of clustered IP3Rs coordinated via a Ca2+-induced Ca2+-release mechanism [52]. Although the cytosolic Ca2+ dependency of IP3Rs has been well characterized, little is known as to how changes in basal cytosolic [Ca2+] would alter the dynamics of IP3-evoked Ca2+ 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 Ca2+ 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 IP3-signaling pathway has been indicated as a potential mechanism responsible for alteration of the Ca2+ homeostasis in AD neurons. We and others have suggested that a G-protein mediated activation of PLC by iAβ42Os is responsible for the overproduction of IP3 and consequent rise of cytosolic Ca2+ observed in cells exposed to iAβ42Os [14, 18]. Moreover, others have suggested that Aβ may cause cytosolic Ca2+ rise by a mixed mechanisms of PLC-dependent and independent manner [15, 16, 55]. The effect of iAβ42Os on intracellular Ca2+ fluxes have previously been investigated by developing a computational model to study important intracellular Ca2+ pathways in normal and in iAβ42Os affected conditions [27]. However, no upstream IP3 production processes were incorporated in the model. Here, we have illustrated a possible mechanistic way for how iAβ42Os triggers IP3 overproduction with consequent rise in cytosolic Ca2+ by including some mechanisms of upstream IP3 production in the model. Specifically, we pinpoint two main possible sites of action for iAβ42Os to interact in the cascade of events resulting from stimulation of G-protein in the plasma membrane to the release of Ca2+ from the ER. In our previous study [18], we argued that it was unlikely that iAβ42Os act on IP3Rs in the generation of Aβ-induced Ca2+ 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 IP3Rs when large doses are introduced. The analysis illustrated in Fig 11 helps us understand what happens to Ca2+ signaling in the presence of iAβ42Os as changes to IP3Rs occur. The persistent increase of iAβ42Os may alter the sensitivity of IP3Rs to Ca2+ over time. For large doses of iAβ42Os, IP3Rs may become more sensitive to low- or sub-threshold IP3 levels and in turn trigger local and global Ca2+ signaling events. The fact that the parameter k−4 appears to play a major role in the decay of observed Ca2+ signals singles out the potential that iAβ42Os do act on the IP3R itself, at least for large doses. Our model suggests the need for further investigation on the relationship between iAβ42Os and the sensitivity of IP3Rs to IP3 levels. Our approach provides a precise way to incorporate the effects of iAβ42Os on IP3 signaling mechanisms that does not necessarily depend on the choice of the IP3R model. When a saturating binding rate model for the IP3R 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β42Os with G proteins and PLC may be sufficiently captured by the model. Additionally, Toglia et al. [27] have also suggested a relationship between IP3 concentration and iAβ42Os. However, their investigation assume that IP3 concentration levels are impacted by iAβ42Os 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β42Os affects the upstream mechanisms in the IP3 signaling cascade. Although our model considers the impact of iAβ42Os specifically on the IP3 signaling cascade in oocytes, our results could be useful in more complex models of various cells. Existing astrocyte models (such as [34, 56–58]) that incorporate Ca2+ dynamics could be altered to include the effects of Aβ on IP3 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 Ca2+ 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 Ca2+ 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 Ca2+ signals over various timescales. Using the model to better understand what happens to Ca2+ regulation in these simulations can directly influence and suggest how one could control Ca2+ 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 Ca2+ signaling in the presence of increasing concentrations of iAβ42Os. 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β42Os 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. Appendix According to Maravall et al. (2000) [28], changes in Ca2+ concentration are associated with changes in fluorescence through the equation (23) where KD is a dissociation constant, fmax is the intensity of the dye at maximum Ca2+ concentration, Rf = fmax/fmin is the indicator’s dynamic range with fmin being the intensity at minimum Ca2+ concentration, δfmax is the saturation of the Ca2+ indicator, fm = fmax/f0, and c0 is the resting Ca2+ concentration. The values of KD and Rf are associated with attributes of the indicator in a particular cellular environment [28] and as such are independent of cellular properties. Wavelength ratio measurements do not generally depend on dye concentration, optical path length, excitation intensity, or detector efficiency [28]. However, the value of KD and the dynamic range of Rf may vary batch to batch and should be estimated using a similar protocol and cellular cytoplasmic domain [28, 59]. In Eq (23), δfmax is the key parameter needed for establishing the conversion from fluorescence to Ca2+ concentrations. When we fix the initial Ca2+ concentration, c0, we can estimate δfmax using (24) as long as true saturation is attained [28] and KD and Rf are known. In practice, δfmax can be used to estimate the unknown resting Ca2+ concentration by inverting the relationship in Eq (24). With Eqs 23 and 24 in hand, converting fluorescence data to Ca2+ concentrations only requires obtaining values for fm, Rf, and δfmax during experimental procedure. However, these values are often not reported in favor of the traditional δf fluorescence measurements and extracting them from reported changes in fluorescence ratio, or establishing their values a posteriori, can be challenging. As such, in order to complete a conversion for data given in terms of δf we approximate a number of parameters. Since both KD and Rf depend on indicator properties, they can be approximated for a variety of indicators. Based on the experiment in [18], we assume values of KD ≈ 0.2−0.5 μM and that Rf has a dynamic range Rf ≈ 85−100 and note that uncertainties in Rf have minimal affect on Eq (24). We illustrate this in Fig 12 where δfmax is plotted as a function of KD and Rf when c0 = 0.01 μM and c0 = 0.05 μM in A and B, respectively. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Uncertainty in estimation of Rf have minimal effect on data rescaling. Effects of KD and Rf on δfmax under initial Ca2+ concentration c0 = 0.01 in A and c0 = 0.05 in B. Notice that Rf has minimal effect on δfmax while KD alters the value of δfmax. https://doi.org/10.1371/journal.pone.0246116.g012 As can be seen from Fig 12, Rf has little effect on the value of δfmax. This is consistent with the idea that for indicators with a large dynamic range, the exact value of Rf is insignificant [28]. For indicators such as Fluo 4, KD is often assumed to be between 0.25 and 0.5 μM but some studies suggest that KD may have much greater range [60, 61]. Without loss of generality, here we consider a basal Ca2+ concentration of c0 = 0.05 μM and set KD = 0.3 μM and Rf = 100. Because we have no previous knowledge for the value of fm, we consider a range fm ≈ 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 Ca2+ concentrations from the fluorescence data in [18] for various estimates of fm. 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 fm ranging from fm = 1 to fm = 100 (black) with n = 11 (fm = 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Impact of fm on rescaling of Ca2+ data. Changes in scaled experimental data when fm ranges from 1 to 100. In all figures, Rf = 100, KD = 0.3, and c0 = 0.05. Figs A, B, and, C correspond to the scaled data for a = 3, a = 10, and a = 30, respectively. The maximum value of each scaled experimental data set is shown by the open circle. The bold color curve corresponds to fm = 40, the value used throughout the simulations. https://doi.org/10.1371/journal.pone.0246116.g013 To study the impact of the conversion to Ca2+ concentrations, Fig 14A–14C shows the corresponding maximum value of the concentration as a function of fm and the range of δfmax 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 fm or c0, 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. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Peak values as a function of scaling parameters. Peak values of the data conversion as a function of fm and δfmax for a = 3 A, a = 10 B, and a = 30 C. In all simulations, Rf = 100, KD = 0.3, and c0 = 0.05 were fixed. This figure shows that the greatest scaling affects occur when fm and δfmax are large and small, respectively. https://doi.org/10.1371/journal.pone.0246116.g014
TI - Quantifying the dose-dependent impact of intracellular amyloid beta in a mathematical model of calcium regulation in xenopus oocyte
JF - PLoS ONE
DO - 10.1371/journal.pone.0246116
DA - 2021-01-28
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/quantifying-the-dose-dependent-impact-of-intracellular-amyloid-beta-in-9vVKQRpARB
SP - e0246116
VL - 16
IS - 1
DP - DeepDyve
ER -