Composite mathematical modeling of calcium signaling behind neuronal cell death in Alzheimer’s disease

Abstract

Background

Alzheimer’s disease (AD) is a progressive neurological disorder, recognized as the most common cause of dementia affecting people aged 65 and above. AD is characterized by an increase in amyloid metabolism, and by the misfolding and deposition of β-amyloid oligomers in and around neurons in the brain. These processes remodel the calcium signaling mechanism in neurons, leading to cell death via apoptosis. Despite accumulating knowledge about the biological processes underlying AD, mathematical models to date are restricted to depicting only a small portion of the pathology.

Results

Here, we integrated multiple mathematical models to analyze and understand the relationship among amyloid depositions, calcium signaling and mitochondrial permeability transition pore (PTP) related cell apoptosis in AD. The model was used to simulate calcium dynamics in the absence and presence of AD. In the absence of AD, i.e. without β-amyloid deposition, mitochondrial and cytosolic calcium level remains in the low resting concentration. However, our in silico simulation of the presence of AD with the β-amyloid deposition, shows an increase in the entry of calcium ions into the cell and dysregulation of Ca ²⁺ channel receptors on the Endoplasmic Reticulum. This composite model enabled us to make simulation that is not possible to measure experimentally.

Conclusions

Our mathematical model depicting the mechanisms affecting calcium signaling in neurons can help understand AD at the systems level and has potential for diagnostic and therapeutic applications.

Background

Alzheimer’s disease (AD) is characterized by the deposition of β-amyloid (A β) oligomers in and around neurons in the brain accompanied by dysfunctional neuronal calcium homeostasis. Autophagy is generally an efficient mechanism for removing amyloids. During the onset of AD, autophagy is increased but the transfer of autophagic vesicles to the lysosomes is blocked [1]. This may contribute to the accumulation of amyloids. There is increasing evidence to support the hypothesis that A β induces an up-regulation of intracellular Ca ²⁺ and leads to AD. Multiple studies on AD mouse models have shown that Ca ²⁺ dysregulation leads to increased Ca ²⁺ entry into the cytoplasm resulting in neuronal cell death and AD [2, 3].

The observations about the effect of A- β oligomers on neuronal calcium signaling led to the formulation of the calcium hypothesis of AD [4]. The basic argument behind the hypothesis is that the activation of the amyloidogenic pathway results in a remodeling of the neuronal calcium signaling pathway. The up-regulation of Ca ²⁺ distorts the normal neuronal Ca ²⁺ signaling by increasing the amount of Ca ²⁺ being taken up by the mitochondria. A sustained increase in the mitochondrial Ca ²⁺ may activate the mitochondria to initiate the intrinsic pathway of Ca ²⁺-induced apoptosis, as described in the calcium hypothesis of Alzheimer’s disease [4].

According to Berridge, the increased output of Ca²⁺ due to the hypersensitivity of the Ca²⁺ signaling system may activate the mitochondria to initiate the intrinsic pathway of Ca²⁺-induced apoptosis by opening up the mitochondrial permeability transition pore (PTP), causing collapse of the mitochondrial membrane potential and releasing cytochrome c and other factors that activate the caspase cascade responsible for apoptosis.

Ichas and Mazat [5] demonstrated that the mitochondrial PTP operates at the crossroads of 2 distinct physiological pathways i.e. the Ca²⁺ signaling network during the life of the cell and the effector phase of the apoptotic cascade during Ca²⁺-dependent cell death. It has 2 open conformations correspondingly. The low-conductance state, which allows the diffusion of small ions like Ca²⁺, is pH-operated, promoting spontaneous closure of the channel. A high-conductance state, which allows the unselective diffusion of big molecules, stabilizes the channel in open conformation [5].

Mitochondria in open high-conductance state can no longer maintain a proton gradient, and thus cannot sustain oxidative phosphorylation, resulting in an arrest of aerobic ATP synthesis (necrotic cell death). This also results in an oncotic imbalance in the mitochondria causing it to swell up. The cristae formed by the inner membrane unfold, leading to rupture of the outer membrane that brings into direct communication the former intermembrane space and the extra-mitochondrial medium. Soluble components like cytochrome c and Apoptosis Inducing Factor (AIF), which are normally trapped in the intermembrane space, are released into the cytosol, thereby inducing cell apoptosis [5].

Thul [6] described the use of Ordinary Differential Equations (ODEs) to model intracellular Ca²⁺ oscillations, assuming intracellular Ca²⁺ concentration to be spatially homogeneous. The use of ODEs is widespread among modelers because: (a) the study of ODEs is computationally well-supported, with a large body of techniques available to investigate ODEs in great detail, and (b) lack of sufficient experimental data to develop a spatially extended model.

In this paper the processes mentioned above have been modeled mathematically using ODEs to allow for quantitative understanding of the dynamics of neuronal cell death in AD. We achieved this by integrating three models: Fall-Keizer Model [7], Mitochondrial PTP Model [8], and Amyloid Metabolism Model [9]. These models are explained in detail in the next section. The results obtained are qualitatively consistent with all the three papers. To demonstrate the validity of our composite model, we tested two hypotheses proposed by these individual models:

• (1) When there is no abnormality in β-amyloid folding and deposition, the initiation of an action potential does not lead to long-term sustained oscillations (an extension of the result detailed in [7]). This will be explained further in the Results section.

• (2) When β-amyloid misfolding affects calcium signaling within the neuron, the action potential is prevented from dying down immediately (a more complex expression of the qualitative trend represented by [9]). As β-amyloid deposition increases with time, the rate of entry of calcium ions into the cell increases, thus allowing the calcium ion oscillations to continue by maintaining high calcium ion levels in the cytosol. As this rate of entry continues to increase, the cytosol and mitochondria attempt to release excess calcium ions to one another more frequently and in smaller amounts, thus resulting in smaller but more rapid oscillations.

We used time-course simulation to relate the occurrence of events to biological processes, thereby verifying our model. We also quantified certain findings from the time-course simulations, with special emphasis on the time taken for the PTP to open in high conductance state.

Methods

The process of building a composite model from individual models is depicted in the flow chart in Fig. 1. Our objective is that the composite model should not only satisfy the properties of the individual models but also help lead biologists using this model to additional insights regarding the processes described here. Specifically, for an in-depth study of calcium signaling we integrate key components from the three models to investigate the dynamics of calcium oscillations to neuronal cell death. The three models are explained in the following subsections.

Fig. 1

Flow chart of analyzing existing models and constructing one composite model

As depicted in the model schematic diagram (Fig. 2), the composite model consists of four main molecular species, namely Cytosolic Ca ²⁺ (CAC), Mitochondrial Ca ²⁺ (CAM), Endoplasmic Reticulum Ca ²⁺ (CAER) and beta-amyloid (A β). In addition, a node is included to represent the mitochondria in high-conductance state (PTPh).

Fig. 2

Composite Model Schematic relating Cytosolic Ca ²⁺ concentration (CAC), Mitochondrial Ca ²⁺ concentration (CAM), Endoplasmic Reticulum Ca ²⁺ concentration (CAER), β-amyloid (a) concentration and PTP high-conductance activation state (PTPh)

The Fall-Keizer model

The Fall-Keizer model is an integrated model depicting mitochondrial Ca ²⁺ handling and metabolic function. It integrates the Magnus-Keizer model [10, 11] and the De Young-Keizer model [12]. The Magnus-Keizer model is a comprehensive mitochondrial model with six proton transfer mechanisms that affect Ca ²⁺ signaling.

A key motivation for using the Fall-Keizer model was to improve the original Magnus-Keizer model, by modifications to the Ca ²⁺ uniporter so that the prediction of Ca ²⁺ signaling can be more accurate.

The inclusion of the De Young-Keizer model for inositol-1,4,5-triphosphate (IP3)-mediated Ca ²⁺ release along with appropriate scaling and provisions for accommodating different cell types allows a modeler to easily shape this model to his or her objective. Furthermore, the comprehensive nature and the modularity of the model that have been maintained by Fall and Keizer make this model an ideal choice for the purposes of this paper.

The model has over 12 variables interacting with each other to depict the Ca ²⁺ signaling in a cell. The basic mitochondrial function has been taken from the Magnus-Keizer model, with a modified mechanism for Ca ²⁺ uptake introduced by Fall and Keizer [7]. Furthermore, the modular nature of this model allows for addition of the PTP characteristics which were omitted in the Fall-Keizer model.

$$ {{\begin{aligned} \frac{dCAC}{dt}= \frac{f_{i}}{V_{c}\tau_{min}}\left[M\left(J_{Na,ex}^{Ca}-J_{uni}^{Ca}-J_{PTP}^{Ca}\right)-E\left(J_{serca}-J_{er,out}\right)\right] \end{aligned}}} $$

(1)

Equation (1) was obtained from the Fall-Keizer model [7]. CAC is directly proportional to the total volume of mitochondria in the cell (M) multiplied by the rate of transfer of calcium ions into the cytosol from the mitochondria, which is dependent on the sodium-calcium ion exchanger, the calcium uniporter and the PTP. CAC is also directly proportional to the total volume of ER in the cell (E) multiplied by the rate of transfer of calcium ions into the cytosol from the ER, which is affected by the SERCA pump and the leakage of calcium ions from the ER.

Since our focus is on mitochondria-induced cell apoptosis, this model provides us with the appropriate foundation to build upon.

The Mitochondrial PTP model

The model proposed by Oster et al. [8] provides a representation of the mitochondrial permeability transition pore (PTP) behavior. In our model, we have focused on the high conductance state of PTP. In its high-conductance conformation, PTP opening induces unselective solute fluxes that dissipate the concentration gradients of relatively big molecules. However, since most proteins remain trapped in the matrix, the resulting oncotic imbalance (at least in vitro) causes high amplitude swelling of the organelle. The subsequent unfolding of the inner membrane causes rupture of the outer membrane, which results in the release of soluble components (mainly cytochrome c and Apoptosis Inducing Factor (AIF)) that are normally located in the intermembrane space [13–15].

The transition of the pore to a high-conductance state requires prolonged levels of high mitochondrial Ca ²⁺. Once the pore opens to this state, it remains open, leading to cell death. The model proposed by Oster et al. [8] assumes that whether or not the pore enters the high-conductance state depends on a secondary slow process, which in turn depends on the overall mitochondrial Ca ²⁺ load.

The Amyloid metabolism model

It is widely known that amyloids perturb Ca ²⁺ homeostasis, and β-amyloids perturb the balance between Ca ²⁺ entry in and extrusion out of the cytoplasm. In healthy neurons, these processes equilibrate, leading to a basal Ca ²⁺ level in the range of 50-100 nM [16]. Studies on the cortical neurons of AD stricken animals found a basal Ca ²⁺ level of around 250 nM, i.e. around twice that found in controls [17].

Based on these observations, we believe that the model proposed by Caluwe and Dupont [9] provides an accurate representation of the relationship between β-amyloid protein concentration and cytoplasmic calcium concentration. Their depiction of a positive feedback loop supports the increased basal Ca ²⁺ levels observed in cortical neurons with the accumulation of amyloids. Furthermore, research conducted on the production of the toxic oligomers found that Ca ²⁺ ions actually promote the synthesis of β-amyloids [18]. These results support the existence of a bistable state switch in the neurons due to the positive feedback loop between the amyloid concentration and the calcium ion concentration. Equation (2), describing the change in intracellular calcium concentration, has been taken from the Caluwe and Dupont model [9].

$$ \frac{dCAC}{dt}=V_{2}+K_{\beta}a^{m} - k_{2} * CAC $$

(2)

This equation shows the effect of amyloids on intracellular calcium concentration. It was formulated based on the assumption that amyloids increase intracellular calcium concentration by increasing the permeability of plasma membranes. This assumption is also made in our model. Here, the rate at which calcium enters the cytoplasm and the rate of eliminating calcium ions from the cytoplasm are assumed to be constant. Furthermore, calcium ion concentration in the cytoplasm is assumed to have first order kinetics.

Composite model

We used XPPAUT [19], an.ode model file as in Additional file 1 and MATLAB to plot our time-course simulations. Details of the composite model constructed based on the afore-mentioned three models are given below. A schematic diagram of the composite model is illustrated in Fig. 2. The parameters used in our composite model were obtained from the individual models.

The model equations for these molecular species are given below and the model parameters are given in Table 1. For the meaning of the rate constants and other parameters readers are to refer to Tables 1, 2 and 3. In an attempt to extend the models described previously and make them more comprehensive, we have combined components and effects that act on the same ion concentration. For instance, here, we used the components affecting CAC in Eqs. (1) and (2) to create the composite equation shown in Eq. (3):

$$ {{\begin{aligned} \frac{dCAC}{dt}\,=\, \frac{f_{i}}{V_{c}\tau_{min}}\!\left[\!M\left(J_{Na,ex}^{Ca}\,-\,J_{uni}^{Ca}-J_{PTP}^{Ca}\right)\,-\,E\left(J_{serca}\,-\,J_{er,out}\right)\right] \,+\, k_{\beta}a^{m} \end{aligned}}} $$

(3)

Table 1 List of parameters used in our model with their values and biological significance [7–9]

Table 2 List of parameters used in IP₃ step function for IP₃-mediated calcium oscillations [7–9]

Table 3 List of parameters used in PTP Integration for PTPh and PTPl states [7–9]

This equation shows the effect of amyloid concentration inducing calcium ion entry into the cell by taking into account its effect of increasing plasma membrane permeability. The effect of amyloid deposition on cytosolic calcium concentration is quantified as the rate of increase of the Ca²⁺ entry multiplied by the amyloid concentration raised to its cooperativity coefficient.

$$ \frac{dCAM}{dt}=\frac{f_{m}M}{V_{m}\tau_{min}}\left(J_{uni}^{Ca}-J_{Na,ex}^{Ca}+J_{PTP}^{Ca}\right) $$

(4)

The dynamics of mitochondrial calcium ion concentration is shown in Eq. (4), obtained from Oster et al. [8]. It increases with the calcium ion influx through the uniporter and PTP, and it decreases with the calcium ion efflux through the Na ⁺-Ca ²⁺ ion exchanger.

$$ \frac{dCAER}{dt}=\frac{f_{i}E}{V_{e}\tau_{min}}\left(J_{serca}-J_{er,out}\right) $$

(5)

ER calcium ion concentration varies with influx through the SERCA (Sarcoplasmic/Endoplasmic Reticulum Calcium ATPase) pump and efflux by leakage of calcium ions, obtained from Fall and Keizer (2001) [7].

$$ \frac{da}{dt}=V_{1}+V\alpha\left(\frac{1}{1+\frac{K_{\alpha}^{n}}{CAC^{n}}}\right)-K_{1}a $$

(6)

Equation (6) focuses on the rate of change of amyloid concentration (a) in the cell. Like Caluwe and Dupont [9], we picked amyloid concentration as a generic quantity, encompassing intracellular and extracellular amyloid concentrations, along with amyloid compounds of different lengths and oligomerization states. The rate of amyloid synthesis is assumed to be a constant. The increase of the cytosolic calcium concentration increases the amyloid concentration as well, according to Caluwe and Dupont [9]. Moreover, the increased amyloid concentration negatively affects its own increase in gradient, allowing some form of regulation of amyloids in spite of cytosolic calcium ion concentration.

$$ \frac{dPTP_{h}}{dt}=\frac{PTP_{h}^{\infty}-PTP_{h}}{\tau_{h}} $$

(7)

Equation (7) is derived from Oster et al. [8]. It depicts the mitochondrial PTP in its high conductance state as a binary state controlled by a Heaviside function. This depends on a secondary slow process, which is also a Heaviside function function that depends on a mitochondrial calcium ion concentration crossing a particular threshold. This is the time constant used to provide a time lag in simulating the actual switching of the PTP state. Tables 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 and 17 contain the equations, parameters and initial conditions from other models that have been used in this model.

Table 4 List of model ODEs from [7–9]

Table 5 List of initial conditions [7–9]

Table 6 List of model equations used in calculation of protein amounts [7–9]

Table 7 List of model equations used in calculation of compartment volumes [7–9]

Table 8 List of model equations used in nucleotide conversion/conservation relations (obtained from [10, 11])

Table 9 List of model equations used in calculation of proportion of free nucleotides (obtained from [11])

Table 10 List of model equations used in calculation of proportion of charged, free nucleotides (obtained from [11])

Table 11 List of model equations used in mitochondrial Ca²⁺ handling (obtained from [10])

Table 12 List of model equations used in calculation of mitochondrial respiration equations (obtained from [11])

Table 13 List of model equations used in calculation of mitochondrial Fo/F1-ATPase equations (obtained from [11])

Table 14 List of model equations used in calculation of cytosolic components (obtained from [11])

Table 15 List of model equations used in calculation of ER Ca²⁺ handling (modified from [21])

Table 16 List of model equations used in IP₃ step function (modified from [21])

Table 17 List of model equations used in PTP Integration (modified from [21])

Results

In the absence of pathology

Figures 3 and 4 show the time-course simulation of CAC, CAM and PTPh in a neuron stimulated by an action potential but devoid of deposition of misfolded β-amyloid oligomers. The β-amyloid concentration has been mathematically made null to depict the normal response of these parameters in a neuron to such a stimulus. In this scenario, there is an initial spike in cytosolic calcium level due to the action potential, which is subsequently followed by oscillations of lower amplitude in the ER, cytosol and mitochondria. The oscillations during such a stimulus are consistent with those observed in the Fall-Keizer model [7]. On adding mitochondria to IP₃-mediated calcium ion release in isolation, long-term oscillations are observed (as detailed in [7]). However in our model, on incorporation of uniporter and PTP transport, there is transport of calcium ions outside this oscillatory mechanism, because of which the oscillations cannot be maintained and will stop. These oscillations are temporary, and die down after around 110 ms, restoring the calcium level in the ER, mitochondria and cytosol to their original values. There is naturally no stimulation of the mitochondrial PTP to open in the high conductance state due to the low resting concentration of Ca ²⁺ in the mitochondria.

Fig. 3

Overall time-course simulation showing cytosolic (CAC) and mitochondrial (CAM) Ca ²⁺ ion concentration along with tendency for opening of PTP in high-conductance state (PTPh) over a time of 500ms in the absence of β-amyloid deposition (by setting V₁ = 0, V _α = 0 and K₁ = 0)

Fig. 4

Addition of ER Ca ²⁺ ion concentration (CAER) to Fig. 3 and a visual representation of the scale of difference between the resting ER Ca ²⁺ level and cytosolic and mitochondrial resting Ca ²⁺ levels

In this case, the neuron responds as it would normally do to a stimulus, with the Ca ²⁺ ions being released rapidly from the ER via the IP₃ channels to maintain calcium homeostasis across the cell. This sudden spurt results in the onset of oscillations of Ca ²⁺ ion concentration in the mitochondria, cytosol and ER. There is no further activity, and the concentration of Ca ²⁺ is sustained at a low equilibrium value.

In the presence of pathology

In Fig. 5, the neuron is affected by β-amyloid deposition, causing an increase in the entry of calcium ions into the cell and dysregulation of Ca ²⁺ channel receptors on the ER. At the start of the time-course simulation, an action potential is simulated causing the calcium levels in the cytosol to jump to over 1.4 μM. Then, owing to Calcium-Induced Calcium Release (CICR), the process of sequestering and subsequent release of Ca ²⁺ ions, with the ER and mitochondria, the calcium level continues to oscillate for a while until it dies down and returns to its original level of around 0.05-0.1 μM. As the amyloid concentration grows in the process, there is a sudden jump in the cytosolic and correspondingly mitochondrial calcium ion concentration due to the sudden release of a large volume of calcium ions from the ER (not shown in this graph). Amyloid metabolism affects the ryanodine and IP₃ receptors located on the surface of the ER that regulate calcium release from the ER [20]. The sudden rise in calcium levels can be attributed to the ER sequestering a large proportion of the calcium that was initially in the cytosol and mitochondria.

Fig. 5

Graph depicting cytosolic (CAC) and mitochondrial (CAM) Ca ²⁺ ion concentration, β-amyloid concentration (a) and tendency for opening of PTP in high-conductance state (PTPh) over a time of 1000ms - in the presence of pathology (using the provided equations, parameter values and initial conditions in the XPP.ode file)

As this continues, oscillation of calcium levels resumes, except this time the mitochondrial calcium ion concentration continuously increases, indicating an increase in sequestration of calcium ions by the mitochondria. As this constant rise is maintained for a period of time, as detailed in Oster et al. [8], a slower secondary process is activated, which results in the PTP opening in the high-conductance state. The opening of PTP in high-conductance state detailed by Oster et al. [8] follows exactly the same pattern, with the high, sustained levels of mitochondrial calcium ion concentration resulting in the activation of a slower, secondary process. This secondary slow process, on completion, results in the PTP opening in high conductance state.

The β-amyloid concentration rises slowly with the growth in cytosolic calcium, and this observation can be attributed to the relatively slow nature of amyloid processing to produce β-amyloids as compared to transfer of Ca ²⁺ ions (Fig. 6).

Fig. 6

Zooming into 0-250 ms of Fig. 5 to visualise the oscillations and growth of amyloid concentration

Figure 7 shows that the PTPh starts opening at around 742 ms and takes 123 ms to reach the high-conductance state at around 865 ms. At the point of opening of the PTP in high conductance state (near 741 ms), mitochondrial Ca ²⁺ rushes out through the PTP and into the cytosol, which in turn is sequestered by the ER. The spurted release results in an oscillation of Ca ²⁺ ions between the ER and the cytosol, which dies down as the mitochondrial calcium level stagnates, indicating that the mitochondria is no longer functional.

Fig. 7

Zooming into 550-1000 ms of Fig. 5 to observe changes on opening of the PTP

When the PTPh opens in the high conductance state, calcium ions rush out of the mitochondria causing the CAM to level off. As the mitochondria stop functioning due to the PTP opening in high conductance, its calcium ion concentration is frozen, with CAC and CAER skyrocketing due to uncontrollable amounts of calcium ions flowing into the cell. This may be because the amyloid concentration has also increased over time (although in a much more regulated fashion), causing the plasma membrane to become more permeable and allowing more calcium ions into the neuron. Consequently, a jump occurs in cytosolic calcium ion concentration, which oscillates for a while and once stabilized begins to increase at an accelerated rate.

Conclusion and discussion

This paper presents a mathematical model of the biological processes involved in the deposition of β-amyloids in and around the neurons and its effect on neuronal calcium signaling homeostasis. Using a detailed mitochondrial model of calcium signaling [7], the paper relates this change in calcium signaling in the cytosol to the calcium level in the mitochondrial matrix. Moreover, the introduction of the permeability transition pore and its characteristics allowed the model to depict the irreversible onset of apoptosis. By combining the models and features regarding β-amyloid deposition, mitochondrial calcium signaling and permeability transition pore activity, we simulated the opening of the permeability transition pore using amyloid deposition as the trigger. We found that, using the parameters and model equations above, high-conductance state is reached at around 865 ms after amyloid deposition begins.

The lack of comprehensive models of subcellular dynamics in Alzheimer’s disease poses a challenge for computational scientists to explore. We envisage a need for integrated modeling, i.e. selecting individual and specialized models and finding ways to combine them, to formulate a composite model of a neuron’s cell fate in AD.

In this paper, we have modeled only a single neuron undergoing mitochondrial PTP-related apoptosis, which is one of the mechanisms of cell death in AD. We have not studied necrosis due to Ca ²⁺ excitotoxicity, or looked at the effect of β-amyloid deposition and misfolding on cognitive functions, both of which can substantially extend the scope of this research. Furthermore, the effect of the neuron on its neighboring neurons in such a scenario has not been considered in this paper, as we assumed the neuron to exist in isolation. These issues are pertinent, and extending our composite model to address them could significantly improve our understanding of the disease.

In summary, we constructed a composite model by integrating three individual models to recapitulate a sequence of events and their repercussions that eventlead to neuronal death. Through this model, we are able to shed new light on a single sequence of processes starting from the deposition and misfolding of β-amyloid in and around the neuron all the way to neuronal apoptosis. Compared with existing models, this model provides a more comprehensive view of these molecular processes occurring in Alzheimer’s disease. It represents a step towards building more realistic models to facilitate the diagnosis and treatment of this aging-related disease.