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

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 2+ 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 2+ and leads to AD. Multiple studies on AD mouse models have shown that Ca 2+ dysregulation leads to increased Ca 2+ 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 2+ distorts the normal neuronal Ca 2+ signaling by increasing the amount of Ca 2+ being taken up by the mitochondria. A sustained increase in the mitochondrial Ca 2+ may activate the mitochondria to initiate the intrinsic pathway of Ca 2+-induced apoptosis, as described in the calcium hypothesis of Alzheimer’s disease [4].

According to Berridge, the increased output of Ca2+ due to the hypersensitivity of the Ca2+ signaling system may activate the mitochondria to initiate the intrinsic pathway of Ca2+-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 Ca2+ signaling network during the life of the cell and the effector phase of the apoptotic cascade during Ca2+-dependent cell death. It has 2 open conformations correspondingly. The low-conductance state, which allows the diffusion of small ions like Ca2+, 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 Ca2+ oscillations, assuming intracellular Ca2+ 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:

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.

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 2+ (CAC), Mitochondrial Ca 2+ (CAM), Endoplasmic Reticulum Ca 2+ (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 2+ concentration (CAC), Mitochondrial Ca 2+ concentration (CAM), Endoplasmic Reticulum Ca 2+ 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 2+ 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 2+ signaling.

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

The inclusion of the De Young-Keizer model for inositol-1,4,5-triphosphate (IP3)-mediated Ca 2+ 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 2+ signaling in a cell. The basic mitochondrial function has been taken from the Magnus-Keizer model, with a modified mechanism for Ca 2+ 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.

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 2+. 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 2+ load.

The Amyloid metabolism model

It is widely known that amyloids perturb Ca 2+ homeostasis, and β-amyloids perturb the balance between Ca 2+ entry in and extrusion out of the cytoplasm. In healthy neurons, these processes equilibrate, leading to a basal Ca 2+ level in the range of 50-100 nM [16]. Studies on the cortical neurons of AD stricken animals found a basal Ca 2+ 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 2+ levels observed in cortical neurons with the accumulation of amyloids. Furthermore, research conducted on the production of the toxic oligomers found that Ca 2+ 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].

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

Table 1

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

Parameters	Value	Biological significance
V	1 ml	Total volume
uMmM	1000	Converts μM to mM
τ min	60	Converts minutes to seconds
p cytosol	0.5	Proportion of volume occupied my cytosol
d cytosol	75 mg/ml	Density of cytosolic protein mg/ml
p mito	0.05	Proportion of volume occupied my mitochondria
d mito	1000 mg/ml	Density of mitochondrial protein mg/ml
p er	0.1	Proportion of volume occupied my ER
d er	1000 mg/ml	Density of ER protein mg/ml
c mito	0.0725 nmolmV∗mg	Mitochondrial calcium concentration
ρ uni	300 nmolmg∗min	Maximum rate of transport through mitochondrial uniporter
ρNaCa	3 nmolmg∗min	Maximum rate of transport through Na+/Ca2+ exchanger
ρ res	0.4	Mitochondrial respiration co-efficient
ρF1	0.7 nmolmg∗min	Mitochondria - Fo/F1 ATPase
[Pi]m	20 mM	Concentration of free phosphates
ρ leak	0.2 nmolmg∗min	Mitochondrial membrane proton leak
J red,basal	20 nmolmg∗min	NADH reduction rate
J max,ANT	900 nmolmg∗min	ATP/ADP antiport flux
glc	1 mM	Glucose concentration in cytosol
J hyd,max	30.1 mM	Cytosol hydrolysis of ATP
V IP3	3000 μM	IP3 receptor volume
V leak	0.1	proportion of leakage from IP3 receptor
J leak	0.1	ER leak
d IP3	0.25 μM	IP3 receptor sensitivity
d ACT	1 μM	IP3 receptor Ca2+ activation constant
d INH	1.4 μM	IP3 receptor Ca2+ inhibition constant
τ	4 s	IP3 receptor inhibitory time constant
V serca	110 nmolmg∗min	SERCA pump flux
k serca	0.4 μM	SERCA pump Ca2+ sensitivity
V 1	0.0065 nM/s	Constant rate of β-amyloid synthesis
V α	0.05 nM/s	Maximal rate of β-amyloid synthesis
K α	120 nM	Half-saturation constant
k β	0.2 nM3/s	Rate constant of increased Ca2+ entry
K 1	0.01 s−1	Rate constant of β-amyloid elimination
n	2	Hill coefficient for activation of β-amyloid synthesis
m	4	Cooperativity coefficient
f m	0.0003	Mitochondrial Ca2+ buffering coefficient
f i	0.01	Cytosolic Ca2+ buffering coefficient

Table 2

List of parameters used in IP3 step function for IP3-mediated calcium oscillations [7–9]

Parameters	Value	Biological significance
b a s e l i n e	0.3 μM	Base concentration
a m p l i t u d e	0.3 μM	Oscillation amplitude
i n i t	10 ms	Initial time
d u r a t i o n	100 ms	Duration of oscillations

Table 3

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

Parameters	Value	Biological significance
C A M ∗	4 μM	Threshold Mitochondrial Ca2+
y ∗	0.8	Secondary process threshold
f_H_M	1.28 x 10 −6	Fast buffering constant for protons in mitochondria
p 1	0.022	Parameter p1
p 2	0.0001	Parameter p2
p 3	0.0231	Parameter p3
p 4	0.0001	Parameter p4
a m p τ	26000	Amplitude for time constant
p 6	0.001	Parameter p6
permlH	3.0	PTP Permeability to protons
p e r m Ca	0.4	PTP Permeability to calcium ions
postptp	2	PTP opening indication constant

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 Ca2+ entry multiplied by the amyloid concentration raised to its cooperativity coefficient.

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 2+ ion exchanger.

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

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.

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]

Equation	Biological significance
dNADHm/dt=(Jred−Jo)∗M/(uMmM∗Vm∗minute)	ODE for mitochondrial NADH concentration change (mM/s)
dADPm/dt=(JANT−Jp,TCA−Jp,F1)∗M/(uMmM∗Vm∗minute)	ODE for mitochondrial ADP concentration change (mM/s)
dADPi/dt=(−JANT∗M+(Jhyd−Jp,gly)∗C)/(uMmM∗Vc∗minute)	ODE for cytosolic ADP concentration change (mM/s)
dPSI/dt=−−JresH+J_H,F1+J_ANT+JPTPH+JLH+2∗J_uni+2∗JPTPCa∗M/c_mito∗minute	ODE for mitochondrial inner membrane voltage
dh/dt=(dINH−(CAC+dINH)∗h)/τ	ODE for change in percentage of closed channels

Table 5

List of initial conditions [7–9]

Parameters	Value	Biological significance
PSI(0)	164 mV	Base potential at time 0
CAM(0)	0.05 μM	Mitochondrial Ca2+ concentration at time 0
CAC(0)	0.05 μM	Cytosolic Ca2+ concentration at time 0
CAER(0)	11 μM	ER Ca2+ concentration at time 0
ADPm(0)	4.46 mM	Mitochondrial ADP concentration at time 0
ADPi(0)	0.028 mM	Cytosolic ADP concentration at time 0
NADHm(0)	0.16 mM	Mitochondrial NADH concentration at time 0
h(0)	95%	Percentage of closed channels at time 0
a(0)	0 mM	β-amyloid concentration at time 0
y(0)	0	Secondary slow process involved in PTP opening at time 0
PTPh(0)	0	PTP closed at time 0
PTPl(0)	0	PTP closed at time 0

Table 6

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

Equation	Biological significance
M=V∗pmito∗dmito	Calculation of mitochondrial protein amount
C=V∗pcytosol∗dcytosol	Calculation of cytosolic protein amount
E=V∗per∗der	Calculation of ER protein amount

Table 7

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

Equation	Biological significance
Vm=(V∗pmito)	Calculation of mitochondrial compartment volumes
Vc=(V∗pcytosol)	Calculation of cytosolic compartment volumes
Ve=(V∗per)	Calculation of ER compartment volumes

Table 8

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

Equation	Biological significance
ATPm=(12∗dmito/uMmM)−ADPm	Mitochondrial ATP concentration
NAD=(8∗dmito/uMmM)−NADHM	Mitochondrial NAD concentration

Table 9

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

Equation	Biological significance
ADPmf=0.8∗ADPm	Unbound mitochondrial ADP concentration
ADPif=0.3∗ADPi	Unbound cytosolic ADP concentration

Table 10

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

Equation	Biological significance
[ADP3−]m=0.45∗ADPmf	Unbound, 3- charged mitochondrial ADP concentration
[ADP3−]i=0.45∗ADPif	Unbound, 3- charged cytosolic ADP concentration
[MgADP−]i=0.55∗ADPif	Unbound, 1- charged cytosolic MgADP− concentration
[ATP4−]i=0.05∗ATPi	Unbound, 4- charged cytosolic ATP concentration
[ATP4−]m=0.05∗ATPm	Unbound, 4- charged mitochondrial ATP concentration

Table 11

List of model equations used in mitochondrial Ca2+ handling (obtained from [10])

Equation	Biological significance
Mitochondrial uniporter equations
MWCnum=(CAC/6)∗((1+(CAC/6))3)	MWC numerator
MWCdenom=((1+(CAC/6))4)+(50/((1+(CAC/0.38))2.8))	MWC denominator
MWC=MWCnum/MWCdenom	MWC fraction value
VuniD=PSI−91/13.35	Uniporter potential exponent
J_uni=ρ_uni∗VuniD∗MWC−CAM∗exp−VuniD/1−exp−VuniD∗1−PTP_h	Rate of transport through uniporter considering PTP in high conductance state
Na+/Ca2+ exchanger equations
VncD=expPSI−91/53.4	Na+/Ca2+ exchanger potential generated
J_nc=ρ_nc∗VncD∗1/1+9.4/30∗∗2∗1/1+0.003∗Dmito/CAM∗1−PTP_h	Rate of Na+/Ca2+ exchange

Table 12

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

Equation	Biological significance
Ares=(1.35e18)∗NADHM0.5/(NAD)0.5	Ares = affinity bracketed expression
VresD=exp0.191∗PSI	Respiration potential generated
Proton pump equations
r1=7e−7	Variable r1
r2=(2.54e−3)∗Ares	Variable r2
r_3=0.639∗VresD	Variable r3
r4=7.58e13+(1.57e−4)∗Ares	Variable r4
r_5=1.73+A_res∗1.06e−17∗VresD	Variable r5
JresH=360∗ρ_res∗r1+r2−r3/r4+r5	Rate of transport through proton pump during respiration
Oxygen consumption rate equations
o1=Ares∗2.55e−3	Variable o1
o2=Ares∗2.00e−5	Variable o2
o_3=0.639∗VresD	Variable o3
o_4=VresD∗A_res∗8.63e−18	Variable o4
o5=(1+Ares∗2.08e−18)∗7.54e13	Variable o5
o_6=1.73+1.06e−17∗A_res∗VresD	Variable o6
Jo=30∗ρres∗(o1+o2−o3+o4)/(o5+o6)	Rate of oxygen consumption during respiration

Table 13

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

Equation	Biological significance
AF1=(1.71e9)∗(ATPm)/(ADPmf∗pim)	AF1 = affinity bracketed expression
VF1D=exp(0.112∗PSI)	ATPase potential generated
F0/F1 ATPase phosphorylation of ADPm
f1=10.5∗AF1	Variable f1
f_2=166∗VF1D	Variable f2
f_3=(4.85e−12)∗A_F1∗VF1D	Variable f3
f4=(1e7+0.135∗AF1)∗275	Variable f4
f_5=(7.74+(6.65e−8)∗A_F1)∗VF1D	Variable f5
Jp,F1=−60∗ρF1∗((f1−f2+f3)/(f4+f5))	Rate of F0/F1 ATPase phosphorylation
J_H,F1=−180∗ρ_F1∗0.213+f_1−169∗VF1D/f_4+f_5	Proton flux due to ATPase
JH,leak=ρleak∗(PSI+24.6)	Mitochondrial membrane proton leak
fPDH=1/(1+(1.1∗(1+(15/(1+(CAM/0.05))2))))	Fraction of activated pyruvate
Jred=Jred,basal+6.3944∗fPDH∗Jgly,total	NADH reduction rate
ATP/ADP antiport flux
ant1=([ATP4−]i/[ADP3−]i)∗([ADP3−]m/[ATP4−]m)∗exp(−PSI/26.7)	Variable ant1
ant2=1+([ATP4−]i/[ADP3−]i)∗exp(−PSI/53.4)	Variable ant2
ant3=1+([ADP3−]m/[ATP4−]m)	Variable ant3
JANT=Jmax,ANT∗((1−ant1)/(ant2∗ant3))	Rate of Adenine Nucleotide Translocator (ANT) activity
Phosphorylation of ADPm from TCA cycle
Jp,TCA=(Jred,basal/3)+0.84∗fPDH∗Jgly,total	Unbound, 3- charged mitochondrial ADP concentration

Table 14

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

Equation	Biological significance
Glycolytic rate based on hexokinase
glynum=(123.3∗(1+1.66∗glc)∗(glc∗ATPi))∗0.0249	Glycolytic rate numerator
glydenom=1+(4∗ATPi)+((1+2.83∗ATPi)∗1.3∗glc)+((1+2.66∗ATPi)∗0.16∗glc2)	Glycolytic rate denominator
Jgly,total=glynum/glydenom	Glycolytic rate
Jp,gly=2∗Jgly,total	Phosphorylation of ADPi from glycolysis
Jhyd=41∗(ATPi)+Jhyd,max/(1+(8.7/glc)2.7)	Cytosolic hydrolysis of ATP

Table 15

List of model equations used in calculation of ER Ca2+ handling (modified from [21])

Equation	Biological significance
Jer,out=(VIP3∗((IP3/(IP3+dIP3))3)∗((CAC/(CAC+dACT))3)∗(h3)+vleak)∗(CAER−CAC)	IP3 receptor and leak
J_serca=V_serca∗CAC^2/kserca2+CAC^2	SERCA pump

Table 16

List of model equations used in IP3 step function (modified from [21])

Equation	Biological significance
stepupf=heav(t−init)	Heaviside step up function
stepdownf=heav(t−(init+duration))	Heaviside step down function
IP3=baseline+amplitude∗(stepupf−stepdownf)	IP3 step function

Table 17

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

Equation	Biological significance
τy=1000∗((1000/cosh(CAM/0.1))+0.1)	Time constant for secondary slow process
τh=τy/8	Time constant for PTP high conductance state
PTPh∞=heavy−y^∗	Heaviside step function for PTP max value
y∞=heav(CAM−CAM∗)	Heaviside step function for y threshold value
dy/dt=(y∞−y)/τy	Secondary slow process involved in opening of PTP high conductance state
dPTP_h/dt=PTPh∞−PTP_h/τ_h	PTP high conductance state dynamics
JPTPH=permlH∗PTP_l∗PSI∗H_M−0.0000000398∗exp−37.434∗PSI/1−exp−37.434∗PSI	Proton flux through PTP in high conductance state
JPTPCa=perm_Ca∗PTP_l∗J_uni∗1−postptp∗PTP_h	Rate of Ca2+ ion transport across PTP
dH_M/dt=f_H_M/τ_h∗JLH+JF1H−JresH+JPTPH	Change in mitochondrial proton concentration
τl=p6+ampτ/cosh((HM−p3)/p4)	Time constant for PTP low conductance state
PTPl∞=0.5∗1+tanhp1−H_M/p_2	Rate of change of polling function
dPTP_l/dt=PTPl∞−PTP_l/τ_l	PTP low conductance state dynamics

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 IP3-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 2+ in the mitochondria.

Fig. 3

Overall time-course simulation showing cytosolic (CAC) and mitochondrial (CAM) Ca 2+ 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 V1 = 0, V = 0 and K1 = 0)

Fig. 4

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

In this case, the neuron responds as it would normally do to a stimulus, with the Ca 2+ ions being released rapidly from the ER via the IP3 channels to maintain calcium homeostasis across the cell. This sudden spurt results in the onset of oscillations of Ca 2+ ion concentration in the mitochondria, cytosol and ER. There is no further activity, and the concentration of Ca 2+ 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 2+ 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 2+ 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 IP3 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 2+ 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 2+ 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 2+ 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 2+ 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 2+ 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.

This file contains the parameters, initial conditions and equations used to generate the results in this paper. It can be run using XPP [19]. (ODE 8.18 kb)