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

Abstract

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

Fig 1. Dynamics and bifurcation structure for constant IP₃.

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.

Fig 2. Model structure and components.

Modeling assumptions for the location of impact of Aβ on the production of IP₃ with key Ca²⁺ signaling mechanisms included in the closed-cell model. The key model assumptions for how Aβ impacts the IP₃ 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.

Fig 3. Model solutions mimic experimental Ca²⁺ 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 K_s is investigated. A shows an increasing Ca²⁺ signal that settles to an increased steady-state when K_s = 0.15. B shows that oscillations in Ca²⁺ can exhibit increasing amplitudes such as those found in Fig 1D in [18]. C and D show that as the value of K_s decreases, the oscillatory patterns of the model reproduce the spike-like Ca²⁺ signals observed in Fig 1E in [18]. E illustrates an oscillatory solution with an increased steady-state Ca²⁺ homeostasis level when K_s is just above the Hopf point. Both D and E show the traces for c_s, y, and p in blue, red, and black, respectively. F shows a simplification of the scaled bifurcation diagram with the bifurcation parameter K_s. Notice that as K_s decreases from the base value of 0.15, a transition from stable steady-states into periodic oscillations occurs through a Hopf bifurcation around HB3≈0.1242. The dynamics around K_s = 0.09 include multiple Hopf bifurcations and has more complex structure than what is presented here.

Fig 4. Dynamics of Ca²⁺ 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 K_s. The bifurcation diagram when K_s = 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 K_s = 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 K_s 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.

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).

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 Ca²⁺ 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).

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 IP₃ Receptor parameters are kept fixed. Notice that altering the remaining parameters (IP₃ Model, PLC, and G-Protein) cannot capture the observed Ca²⁺ signal. B shows a “best” fit solution when only the IP₃R parameters are kept fixed. Notice that without altering the IP₃R parameters, model solution peaks and decay also do not reproduce the observed experimental behaviors for large doses of Aβ.

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 IP₃ and PLC, respectively. D shows the unscaled Ca²⁺ 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.

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.

Fig 10. PRCC prediction on solution amplitude and time to peak for model parameter V_Q.

The impact of V_Q 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 V_Q takes on various values ranging from V_Q = 30 to V_Q = 480. The black trace in each diagram represents the baseline V_Q value for the large doses parameter set.

Fig 11. Impact of large doses of Aβ on IP₃R dynamics.

The impact of the IP₃R model parameters K₁ and k₋₄ 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.

Fig 12. Uncertainty in estimation of R_f have minimal effect on data rescaling.

Effects of K_D and R_f on δf_max under initial Ca²⁺ concentration c₀ = 0.01 in A and c₀ = 0.05 in B. Notice that R_f has minimal effect on δf_max while K_D alters the value of δf_max.

Fig 13. Impact of f_m on rescaling of Ca²⁺ data.

Changes in scaled experimental data when f_m ranges from 1 to 100. In all figures, R_f = 100, KD = 0.3, and c₀ = 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 f_m = 40, the value used throughout the simulations.

Fig 14. Peak values as a function of scaling parameters.

Peak values of the data conversion as a function of f_m and δf_max for a = 3 A, a = 10 B, and a = 30 C. In all simulations, R_f = 100, KD = 0.3, and c₀ = 0.05 were fixed. This figure shows that the greatest scaling affects occur when f_m and δf_max are large and small, respectively.