A Multi-Scale Computational Model of Excitotoxic Loss of Dopaminergic Cells in Parkinson's Disease

Parkinson's disease (PD) is a neurodegenerative disorder caused by loss of dopaminergic neurons in substantia nigra pars compacta (SNc). Although the exact cause of cell death is not clear, the hypothesis that metabolic deficiency is a key factor has been gaining attention in recent years. In the present study, we investigated this hypothesis using a multi-scale computational model of the subsystem of the basal ganglia comprising the subthalamic nucleus (STN), globus pallidus externa (GPe), and SNc. The proposed model is a multiscale model in that interaction among the three nuclei are simulated using more abstract Izhikevich neuron models, while the molecular pathways involved in cell death of SNc neurons are simulated in terms of detailed chemical kinetics. Simulation results obtained from the proposed model showed that energy deficiencies occurring at cellular and network levels could precipitate the excitotoxic loss of SNc neurons in PD. At the subcellular level, the models show how calcium elevation leads to apoptosis of SNc neurons. The therapeutic effects of several neuroprotective interventions are also simulated in the model. From neuroprotective studies, it was clear that glutamate inhibition and apoptotic signal blocker therapies were able to halt the progression of SNc cell loss when compared to other therapeutic interventions, which only slowed down the progression of SNc cell loss.


INTRODUCTION
Parkinson's disease (PD) is predominantly considered as a motor disorder, which affects more than 6 million people around the world (Chaudhuri et al., 2006). It is caused by the loss of dopaminergic neurons in substantia nigra pars compacta (SNc) situated in the midbrain region (Fu et al., 2018). The cardinal symptoms of PD, such as tremor, rigidity, bradykinesia, and postural instability (Goldman and Postuma, 2014) are thought to be considered as the first sign of PD pathogenesis. However, other symptoms such as anosmia (loss of smell) (Omori and Okutani, 2020), constipation (Lubomski et al., 2020), sleep disorders (specifically rapid eye movement behavior sleep disorder) (Postuma et al., 2019), and depression (Bayram et al., 2020) also emerge well before motor impairments. Recently, several studies pointed to the fact that even before the emergence of motor and non-motor symptoms, pathogenesis begins with metabolic abnormalities occurring at different levels of neural hierarchy: subcellular, cellular, and network levels (Bolam and Pissadaki, 2012;Pissadaki and Bolam, 2013;Pacelli et al., 2015;Limphaibool et al., 2018;Nam et al., 2018;Muddapu et al., 2019. With the help of a computational model, Muddapu et al. (2019) have recently suggested that the excitotoxic loss of SNc cells might be due to energy deficiencies occurring at different levels in the neural hierarchy.
If metabolic factors are indeed the deep underlying reasons behind PD pathogenesis, it is a hypothesis that deserves closer attention because any positive proof regarding the role of metabolic factors puts an entirely new spin on PD research. Unlike current therapeutic approaches that manage the symptoms rather than provide a cure, the new approach can, in principle, point to a more lasting solution. If inefficient energy delivery or energy transformation mechanisms are the reason behind degenerative cell death in PD, relieving the metabolic load on the vulnerable neuronal clusters, by intervening through surgical and/or pharmacological approaches could prove to be a decisive treatment for PD, a disease that previously proved itself to be intractable to current therapeutic approaches.
In this work, using computational models, we investigated the hypothesis that the major cause of SNc cell loss in PD might be due to energy deficiency occurring in SNc neurons. In the proposed modeling study, we focused on excitotoxicity in SNc caused by the subthalamic nucleus (STN), which is precipitated by energy deficiency and on understanding the mechanism behind neurodegeneration during excitotoxicity. Moreover, it aims to suggest therapeutic interventions that can reduce the metabolic burden on the SNc neurons, which in turn can delay the progression of SNc cell loss in PD. In Muddapu and Chakravarthy (2017), we proposed a preliminary computational spiking network model of STN-mediated excitotoxicity in SNc with an abstract treatment of apoptosis. Building on the previous version of the model in Muddapu et al. (2019), we have improved the excitotoxicity model by incorporating more biologically plausible dopamine plasticity and also explored the therapeutic strategies to slow down or halt SNc cell loss. In , we proposed a comprehensive computational model of SNc cell which helped us to understand the pathophysiology of neurodegeneration at a subcellular level in PD. The previous excitotoxicity models (Muddapu and Chakravarthy, 2017;Muddapu et al., 2019) use a slightly abstract single neuron model of SNc with a simplified representation of programmed cell death and explored the network level interactions that possibly lead to cell death. In the present paper, we used a biophysically elaborate model of SNc neurons and studied the link between energy deficiency and molecular level changes that occur inside SNc neurons under PD conditions. Specifically, we considered intracellular molecular mechanisms such as energy metabolism, dopamine turnover processes, calcium buffering mechanisms, and apoptotic signal pathways. Using such a hybrid model, we were able to show excitotoxic loss of SNc neurons, which was precipitated by energy deficiency and postulated the possible mechanism behind excitotoxic neurodegeneration. Moreover, we also explored various therapeutic interventions to halt or slow down the progression of SNc cell loss.

METHODS
The proposed hybrid excitotoxicity model (HEM) in PD consisted of three nuclei from basal ganglia, namely STN, SNc, and globus pallidus externa (GPe). The SNc neuron was modeled as a biophysical neuronal model, and STN and GPe neurons (Muddapu et al., 2019) were modeled using the Izhikevich neuron model (Izhikevich, 2007). Neurons in each nucleus were arranged as a two-dimensional lattice (Figure 1). All the simulations were carried out on the MATLAB (RRID:SCR_001622) platform, where all models were numerically integrated with a time step of 0.1 ms.

Izhikevich (Spiking) Neuronal Model (STN, GPe)
The Izhikevich neuronal models are capable of exhibiting biologically realistic firing patterns at relatively low computational expense (Izhikevich, 2003). The proposed model of HEM consisted of GPe and STN neurons modeled as Izhikevich spiking neuron models, where the Izhikevich neuronal model parameters were adopted from the literature (Michmizos and Nikita, 2011;Mandali et al., 2015). The neuronal population sizes in the model were selected based on the anatomical data of rat basal ganglia (Oorschot, 1996;Arbuthnott and Wickens, 2007). The external bias current (I x ) was adjusted to match the firing rate of nuclei with published data (Tripathy et al., 2015). The Izhikevich neuron model of GPe and STN consisted of two variables, namely membrane potential (v x ), and the membrane recovery variable (u x ): Resetting: if where, v x ij , and u x ij are the membrane potential, and the membrane recovery variables of a neuronal type x at the location i, j , respectively, I x ij , and I syn ij are the external bias, and total synaptic currents received to a neuronal type x at the location i, j , respectively, C x is the membrane capacitance of a neuronal type x, a, b, c, d are Izhikevich neuronal model parameters, v x peak is the maximum (peak) membrane voltage set to a neuronal type x (where x = STN or GPe).

Biophysical (Conductance-Based) Neuron Model (SNc)
The biophysical neuron model of SNc in the proposed HEM was adopted from  cellular and molecular processes such as ion channels (active ion pumps, ion exchangers), calcium buffering mechanisms (calcium-binding proteins, organelles sequestration of calcium), energy metabolism pathways (glycolysis and oxidative phosphorylation), dopamine turnover processes (synthesis, storage, release, reuptake, and metabolism), and apoptosis (endoplasmic reticulum-stress and mitochondrial-induced apoptosis). The dynamics of SNc membrane potential v SNc is given as, where, F is the Faraday's constant, C SNc is the SNc membrane capacitance, v cyt is the cytosolic volume, A pmu is the cytosolic area, J m,Na is the sodium membrane ion flux, J m,Ca is the calcium membrane ion flux, J m,K is the potassium membrane ion flux, and J inp is the overall input current flux. A more detailed description of the SNc neuron model was provided in . In the proposed model, intracellular calcium concentration in the SNc neuron was dependent on calcium-binding proteins, mitochondria (MT), and endoplasmic reticulum (ER) (Zündorf and Reiser, 2011;Zaichick et al., 2017 where, J m,Ca is the flux of calcium ion channels, J calb is the calcium buffering flux by calbindin, J cam is the calcium buffering flux by calmodulin, J serca,er is the calcium buffering flux by ER uptake of calcium through sarco/endoplasmic reticulum calcium-ATPase, J ch,er is the calcium efflux from ER by calcium-induced calcium release mechanism, J leak,er is the calcium leak flux from ER, J mcu,mt is the calcium buffering flux by MT uptake of calcium through mitochondrial calcium uniporters (MCUs), and J out,mt is the calcium efflux from MT through sodium-calcium exchangers, mitochondrial permeability transition pores (PTPs), and nonspecific leak flux.

Synaptic Connections
The synaptic connectivity among different neuronal populations was modeled as a standard single exponential model of postsynaptic currents (Humphries et al., 2009) as follows: The N-Methyl-D-aspartic Acid (NMDA) current was regulated by voltage-dependent magnesium channels which were modeled as, where, h x→y ij is the gating variable for the synaptic current from neuronal type x to neuronal type y (where x → y = {STN → STN, STN → SNc, STN → GPe, GPe → STN, GPe → GPe), τ Recep is the decay constant for the synaptic receptor, S x ij is the spiking activity of a neuronal type x at time t, W x→y is the synaptic weight from neuron x to y, v y ij is the membrane potential of the neuronal type y at the location (i, j), E Recep is the receptor-associated synaptic potential (Recep = GABA/AMPA/NMDA), and Mg 2+ is the magnesium ion concentration. The time constants of Gamma-Amino Butyric Acid (GABA), Alpha-amino-3-hydroxy-5-Methyl-4-isoxazole Propionic Acid (AMPA), and NMDA in STN and GPe were chosen from Götz et al. (1997) are given in Table 1.

Lateral Connections
The lateral connections in SNc, STN, and GPe, were modeled as Gaussian neighborhoods (Muddapu et al., 2019), where, w m→m ij,pq is the weight of lateral connection strength of a neuronal type m at the location (i, j), d ij,pq is the distance of neuron at the location (i, j) from the center neuron (p, q), R m is the standard deviation of Gaussian, and A m is the amplitude of lateral synaptic strength (where m = GPe or STN or SNc).
The lateral connections within SNc and GPe populations were considered as inhibitory and within STN as excitatory (Muddapu et al., 2019) (Figure 1). The lateral currents in the STN and GPe were modeled similar to Equations (6-8) and in the case of SNc which was modeled as, where, I x→y ij is the synaptic current from neuronal type x to neuronal type y, W x→y is the weight of synaptic strength from neuronal type x to neuronal type y, v x ij , and v y ij are the membrane potential of the neuronal type x and y, respectively, at the location (i, j), E GABA is the GABAergic receptor reversal potential, and s x→y ij is the synaptic gating variable from neuronal type x and y at the location (i, j) (where x → y = {SNc → SNc} only). The parametric values of α, β, θ g , θ H g , σ H g are adopted from Rubin and Terman (2004) and given in Table 2.

Effect of Dopamine on Synaptic Plasticity
Dopamine (DA) modulated lateral connection strength in STN, SNc, and GPe populations. As DA levels increased, the lateral connection strength in GPe and SNc increased, whereas, in the case of STN, it decreased. DA-modulation of lateral connection strength is modeled as, A GPe = s GPe min * e (cdgpe * DA s (t)) where, s STN max , s GPe min , and s SNc min are lateral connection strengths of STN, GPe, and SNc, respectively, where dopamine effect on the basal spontaneous activity of the neuronal population was minimal, cd stn , cd gpe , and cd snc are the parameter values which modulated the influence of dopamine on the lateral connections in STN, GPe, and SNc populations, respectively, DA s (t) is the instantaneous DA level. The instantaneous DA was estimated by the spatial average DA concentration of all the terminals at a given instant. All parameter values are given in Table 3.
The post-synaptic effects of DA in SNc, STN, and GPe were modeled similar to Muddapu et al. (2019), Frontiers in Neuroinformatics | www.frontiersin.org is the parameter that affects the post-synaptic current, and DA s (t) is the instantaneous dopamine level.

SNc
The total synaptic current received by an SNc neuron at the lattice position i, j is the summation of the glutamatergic input from the STN neurons, considering both NMDA and AMPA receptor activation, and lateral GABAergic current from other SNc neurons. is the lateral GABAergic current from other SNc neurons; F STN→SNc is the scaling factor for the glutamatergic current from STN neuron.

GPe
The total synaptic current received by a GPe neuron at the lattice position i, j is the summation of the glutamatergic input from the STN neurons considering both NMDA and AMPA receptors activation and the lateral GABAergic current from other GPe neurons.

Calcium-Induced Neurodegeneration in SNc
Calcium plays a dual role in living organisms-as a survival factor or a ruthless killer (Orrenius et al., 2003). For the survival of neurons, minimal (physiological) levels of glutamate stimulation are required. Under normal conditions, calcium concentration within a cell is tightly regulated by pumps, transporters, calciumbinding proteins, ER, and MT (Wojda et al., 2008;Surmeier et al., 2011). Due to prolonged calcium influx driven by excitotoxicity, calcium homeostasis within the cell is disrupted, which results in cellular imbalance, leading to activation of apoptotic pathways (Bano and Ankarcrona, 2018). The SNc soma undergoes degeneration when there is an abnormal calcium build up inside the cell resulting in calcium loading inside ER and MT, which leads to ER stress-induced and MT-induced apoptosis, respectively (Malhotra and Kaufman, 2011). The proposed model will also include the apoptotic processes inside SNc neurons, which get activated when calcium levels in the neuron cross a certain threshold as a result of overexcitation and/or metabolic deficiency. In the proposed HEM model, we incorporated a mechanism of programmed cell death, whereby an SNc neuron under high stress (high calcium levels) kills itself. The stress in a given SNc neuron is observed by monitoring the intracellular calcium concentrations in the cytoplasm, ER, and MT.
The SNc neuron undergoes ER-stress-induced apoptosis when calcium levels in ER cross a certain threshold (ER thres ). Under such conditions, the particular SNc neuron gets eliminated as follows, where, Ca ER ij is the calcium concentration in the ER, ER thres is the calcium concentration threshold after which ER-stress induced apoptosis is initiated ER thres = 2.15 × 10 −3 mM , v SNc ij is the membrane voltage of neuron at the lattice position i, j .
The ER calcium concentration ([Ca er ]) dynamics is given by, where, β er is the ratio of free calcium to total calcium concentration in the ER, ρ er is the volume ratio between the ER and cytosol, J serca,er is the calcium buffering flux by ER uptake of calcium through SERCA, J ch,er is the calcium efflux from ER by CICR mechanism, and J leak,er is the calcium leak flux from ER. The SNc neuron undergoes mitochondria-induced apoptosis when calcium levels in mitochondria cross a certain threshold (MT thres ). Then that particular SNc neuron will be eliminated as follows, where, Ca MT ij is the calcium concentration in mitochondria, MT thres is the calcium concentration threshold after which mitochondria-induced apoptosis gets initiated (MT thres = 0.0215 mM), and v SNc ij is the membrane voltage of neuron at the lattice position i, j .
The MT calcium concentration ([Ca mt ]) dynamics is given by, where, β mt is the ratio of free calcium to total calcium concentration in the ER, ρ mt is the volume ratio between the MT and cytosol, J mcu,mt is the calcium buffering flux by MT uptake of calcium through MCUs, and J out,mt is the calcium efflux from MT through sodium-calcium exchangers, mPTPs and non-specific leak flux.

Neuroprotective Strategies
Any therapeutic interventions which result in the slowdown of SNc cell loss can be considered as neuroprotective in nature. It can be achieved either by a blockage of glutamatergic receptors on SNc or by attenuating the pathological oscillations in STN pharmacologically or surgically (Rodriguez et al., 1998). Just as in the previous model (Muddapu et al., 2019), glutamate inhibition, dopamine restoration, subthalamotomy, and deep brain stimulation therapies were implemented in the proposed HEM model. In addition to these therapeutics, since the present model captured several molecular mechanisms at the subcellular level, we simulated therapeutic interventions at the cellular level such as calcium channel blockers, enhancement of calcium-binding proteins (CBPs) expression, and apoptotic signal blockers.

Glutamate Inhibition Therapy
In glutamate inhibition therapy (such as MK-801, NBQX, LY-404187), the excitatory drive from STN neurons to SNc neurons is reduced by altering synaptic weight from STN to SNc (W STN→SNc ) (Johnson et al., 2009;Zhang et al., 2019). In the proposed excitotoxicity model, the glutamate inhibition therapy was implemented by the following criterion, where, W STN→SNc (N sc , t) is the change in synaptic weight of STN to SNc connection based on the number of SNc cells surviving (N sc (t)) at a particular time (t), W 0 STN→SNc is the initial synaptic weight of STN to SNc connection, N sc (t) is the number of SNc neurons surviving at a particular time (t), δ GII is the extent of glutamate inhibition, pcl is the percentage of SNc cell loss (25 %) at which intervention of a particular treatment was employed pcl = 0.25 , P snc is the population size of SNc neurons, and T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed. In the present study, the therapeutic intervention was given at 25% SNc cell loss.

Dopamine Restoration Therapy
In dopamine restoration therapy (such as levodopa, DA agonists), the dopamine drive to STN neurons is restored by the increased administration of external dopamine (δ DAA ) (Vaarmann et al., 2013). In the proposed excitotoxicity model, the dopamine restoration therapy was implemented by the following criterion, where, DA (N sc , t) is the change in dopamine level based on the number of SNc cells surviving at a particular time (t) (N sc (t)), DA s (t) is the instantaneous DA level, N sc (t) is the number of SNc neurons surviving at a particular time (t), δ DAA is the extent of dopamine content restoration, and T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed.

Subthalamotomy
In subthalamotomy therapy, the excitatory drive from STN neurons to SNc neurons is reduced by silencing or removing STN neurons (δ les ) (Wallace et al., 2007;Jourdain et al., 2014). In the proposed excitotoxicity model, the subthalamotomy therapy was implemented by the following criterion, where, δ LES is the percentage of STN lesioned which is in the following range: {5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}, N sc (t) is the number of SNc neurons surviving at a particular time (t), and T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed.

Deep Brain Stimulation (DBS) in STN
In our previous study (Muddapu et al., 2019), we proposed that the neuroprotective effect of DBS was due to increased axonal and synaptic failures in the stimulation site (Rosenbaum et al., 2014). Besides, we also suggested that a biphasic current with four contact point (FCP) stimulation configuration showed a maximal neuroprotective effect. The DBS parameters such as amplitude (A DBS ), frequency f DBS = 1 T DBS , and pulse width (δ DBS ) were adjusted by using clinical settings as a constraint (Moro et al., 2002;Garcia et al., 2005). The biphasic current waveform (P BW ) was generated as the following, Frontiers in Neuroinformatics | www.frontiersin.org where, t k are the onset times of the current pulses, A DBS is the amplitude of the current pulse, and δ DBS is the current pulse width. DBS stimulation with biphasic current waveform and four contact point configuration is given as, In the proposed excitotoxicity model, the DBS therapy was implemented by the following criterion, where, I DBS−STN ij (N sc , t) is the stimulation current to STN neuron at position i, j based on the number of SNc cells surviving at a particular time (t), N sc (t) is the number of SNc neurons surviving at a particular time (t), and T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed.
DBS therapy in STN based on the disruptive hypothesis proposed earlier (Rosenbaum et al., 2014;Muddapu et al., 2019) was implemented by the following criterion, where, W STN→SNc (S DBS , t) is the change in synaptic weight of STN to SNc based S DBS = {ON, OFF} at a particular time (t), S DBS is indicator function which modulates the activation and inactivation of DBS, W ASF is the weight matrix based on the percentage of axonal and synaptic failures (δ ASF = 5, 10,20,30,40,50,60,70,80,90,100). The DBS parameters are given in Table 4.

Calcium Channel Blockers
It was reported that calcium channels in the SNc neuron contribute to neurodegeneration in PD (Benkert et al., 2019). In calcium channel blocker therapy (involving calcium blockers such as dihydropyridine, amlodipine), the excess calcium influx into SNc neurons is blocked by reducing the flux through calcium channels in the proposed model (Ritz et al., 2010;Liss and Striessnig, 2019). In the proposed excitotoxicity model, the calcium channel blocker therapy was implemented by the following criterion, where, I CaL ij (N sc , t) is the instantaneous calcium current of SNc neuron at position i, j based on the number of SNc cells surviving at a particular time (t), N sc (t) is the number of SNc neurons surviving at a particular time (t), δ CB is the extent of calcium channel blockers, and T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed.

Enhancement of CBP Expression
It was reported that the expression of CBPs, in general, is very low in the case of SNc neurons (Chung et al., 2005;Greene et al., 2005;Mendez et al., 2005) and reduces even further under PD (Yamada et al., 1990;Hurley et al., 2013;Zaichick et al., 2017;Fairless et al., 2019). Overexpression of CBPs was found to be neuroprotective in PD (Yuan et al., 2013;Inoue et al., 2019;McLeary et al., 2019). In the enhancement of CBPs expression, the calciumbinding proteins such as calbindin (Calb) and calmodulin (Cam) concentration was increased by the administration of external CBP (δ ECP ). In the proposed excitotoxicity model, the enhancement of CBPs expression therapy was implemented by the following criterion, where, CBP x ij (N sc , t) is the CBP concentration based on the number of SNc cells surviving at a particular time (t) (N sc (t)), CBP x ij (t) is the instantaneous CBP concentration where x = {Calb, Cam}, δ ECP is the extent of CBP enhancement, T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed, and N sc (t) is the number of SNc neurons surviving at a particular time (t ).

Apoptotic Signal Blockers
It is known that neurons which are under stress (increased calcium levels) undergo neurodegeneration through apoptosis (Michel et al., 2016). In apoptotic signal blocker therapy (such as azilsartan, TCH346, and CEP-1347), the apoptotic signaling pathways are blocked by inhibiting the activation of caspase 9 and caspase 3 (Yacoubian and Standaert, 2009;Perier et al., 2012;Gao et al., 2017). In the proposed excitotoxicity model, the apoptotic signal blocker therapy was implemented by the following criterion, where, APOP ij (N sc , t) is the apoptotic signal based on the number of SNc cells surviving at a particular time (t) (N sc (t)), APOP ij (t) is the instantaneous apoptotic signal, δ AB is the extent of apoptotic signal blockage, T l represents the number of SNc cells surviving at which intervention of a particular treatment was employed, and N sc (t) is the number of SNc neurons surviving at a particular time (t ).

RESULTS
We have investigated the neuronal model of SNc for their characteristic firing patterns (Figure 2) and studied the effect of energy deficiency on the model response (Figure 3). We then extensively studied the effect of energy deficiency and overstimulation by STN on the survivability of SNc neurons (Figure 4). Finally, we have explored various therapeutics such as glutamate inhibition, dopamine restoration, subthalamotomy, deep brain stimulation, calcium channel blockers, enhancement of CBPs, and apoptotic signal blockers (Figures 5, 6) which slows down SNc cell loss resulting in neuroprotection.

Characteristic Firing of SNc Neurons
SNc neurons exhibit two distinct firing patterns: low-frequency irregular tonic or background firing (3 − 8 Hz) and highfrequency regular phasic or bursting (∼ 20 Hz) (Grace and Bunney, 1984a,b). In the proposed HEM model, the SNc neurons showed spontaneous (tonic) spiking with a firing rate of ∼  4 Hz (Figure 2A) (Grace and Bunney, 1984b). The intracellular calcium concentration of the SNc neuron during resting state was ∼ 1 × 10 −4 mM, and it rose to higher than 1 × 10 −3 mM upon arrival of the action potential ( Figure 2D) (Dedman and Kaetzel, 1997;Ben-Jonathan and Hnasko, 2001;Wojda et al., 2008). Dopamine released by the SNc neuron during tonic spiking was peaked at ∼ 45 × 10 −6 mM, which was in the range of (34 − 48) × 10 −6 mM observed experimentally ( Figure 2G) (Garris et al., 1997). Under energy deficiency conditions (where intracellular ATP concentration was clamped at 0.411), SNc neurons exhibited tonic bursting with two spikes per burst in the absence of external current ( Figure 2B) (Grace and Bunney, 1984a). The intracellular calcium concentration of the SNc neuron during resting state was below 1 × 10 −4 mM, and rose up to 1 × 10 −3 mM upon arrival of the action potential similar to tonic spiking ( Figure 2E). Dopamine released by the SNc neuron during tonic bursting, at its peak was as low as ∼ 10 × 10 −6 mM ( Figure 2H), which was in the range of (7 − 30) × 10 −6 mM  observed experimentally (Chen, 2005;Koshkina, 2006). When an external current was applied (continuous monophasic current with duration 200 ms and amplitude of 100 × 10 −6 pA), the SNc neuron exhibited phasic bursting with more than two spikes per burst ( Figure 2C). The intracellular calcium concentration of the SNc neuron during resting state was ∼ 1×10 −4 mM, and rose to a higher value than 1 × 10 −3 mM upon stimulation ( Figure 2F). Dopamine released by the SNc neuron during phasic bursting peaked as high as ∼ 150 × 10 −6 mM (Figure 2I), which is in the range of (150 − 400) × 10 −6 mM observed experimentally (Schultz, 1998).

The Effect of Energy Deficiency on SNc Neurons
Under energy-deficient conditions, the SNc neuron exhibited bursting with a decreased firing rate (Figures 3A,B). During bursting, the average intracellular calcium concentration of the SNc neuron rose to a peak value of 2.5 × 10 −3 mM ( Figure 3C). Under energy-deficient conditions, average extracellular dopamine concentration released by the SNc neuron was very low (∼ 4 × 10 −6 mM) compared to normal conditions (∼ 14 × 10 −6 mM) ( Figure 3D).

The Effect of Energy Deficiency on Excitotoxicity in SNc
To understand the effect of energy deficiency on excitotoxicity in SNc, glucose and oxygen levels were reduced to very low values.
As the percentage of SNc cells in energy deficiency increased, the percentage of surviving SNc cells decreased in a sigmoidal manner ( Figure 4A). When there were more than 40% of SNc cells under energy-deficient conditions, there was a significant decrease in the percentage of surviving SNc cells.

The Effect of Energy Deficit and STN on the SNc
In order to understand the influence of STN in SNc excitotoxicity under energy deficiency, glucose and oxygen levels were changed to very low values along with varying synaptic weight from STN to SNc (

Neuroprotective Strategies
The proposed hybrid excitotoxicity model was extended to study the effect of different therapeutic interventions on the progression of SNc cell loss. The following therapeutic interventions were simulated: glutamate inhibition, dopamine restoration, subthalamotomy, deep brain stimulation, calcium channel blockers, enhancement of calcium-binding proteins, and apoptotic signal blockers. During the implementation of neuroprotective strategies, the percentage of SNc cells under energy deficiency was fixed at 80%, and synaptic weight from STN to SNc was fixed at 0.3.

Glutamate Inhibition Therapy
According to the protocol mentioned in the methods section, the neuroprotective effect of glutamate agonists and antagonists on the percentage loss of SNc cells was implemented. The glutamate inhibition therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 5A) for varying extent of glutamate inhibition (δ GII ). As the extent of glutamate inhibition increased, the percentage of surviving SNc cells increased in a linear fashion (Figure 5A). At lower values of δ GII , the percentage of surviving SNc cells was below 50%, and at higher values, it was around 75%, which was same as when the therapeutic intervention was initiated.

Dopamine Restoration Therapy
According to the protocol mentioned in the methods section, the neuroprotective effect of dopamine agonists on the percentage loss of SNc cells was implemented. The dopamine restoration therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 5B) for varying extent of dopamine restoration (δ DAA ). As the extent of dopamine restoration increased, the percentage of surviving SNc cells increased in a linear fashion with noisy response (Figure 5B). At lower values of δ DAA , the percentage of surviving SNc cells was below 50%. At higher values of δ DAA (> 1.1), the percentage of surviving SNc cells was around 60%. However, at intermediate levels of δ DAA (0.9 < δ DAA ≤ 1.1), the percentage of surviving SNc cells peaked close to 75%.

Subthalamotomy
According to the protocol mentioned in the methods section, the neuroprotective effect of subthalamotomy on the percentage loss of SNc cells was implemented. The subthalamotomy therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 5C) for varying extent of STN lesioning (δ LES ). As the extent of subthalamotomy increased, the percentage of surviving SNc cells increased in a hyperbolic fashion (Figure 5C). At lower values of δ LES , the percentage of surviving SNc cells was below 50%. With just δ LES value of 20%, the percentage of surviving SNc cells was increased to > 60%. For δ LES values above 20%, the percentage of surviving SNc cells was steadily increased to a maximum value of 73% for δ LES value of 100%.

Deep Brain Stimulation (DBS) in STN
According to the protocol mentioned in the methods section, the neuroprotective effect of DBS on the percentage loss of SNc cells was implemented. The DBS therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 5D) for varying extent of axonal and synaptic failures (δ ASF ). As the extent of axonal and synaptic failures increased, the percentage of surviving SNc cells increased in a linear fashion (Figure 5D). At lower values of δ ASF , the percentage of surviving SNc cells was below 50%, and at higher values, it peaked at 73% for δ ASF value of 100%.

Calcium Channel Blockers
According to the protocol mentioned in the methods section, the neuroprotective effect of calcium channel blockers on the percentage loss of SNc cells was implemented. The calcium channel blocker therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 6A) for varying extent of calcium channel blockers (δ CCB ). As the extent of calcium channel blockers increased, the percentage of surviving SNc cells increased in a sigmoidal fashion (Figure 6A). At lower values of δ CCB , the percentage of surviving SNc cells was below 50%. At δ CCB values of higher than 0.5, the percentage of surviving SNc cells was increased to > 70%.

Enhancement of CBP Expression
According to the protocol mentioned in the methods section, the neuroprotective effect of CBP enhancement on the percentage loss of SNc cells was implemented. The enhancing CBP therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 6B) for varying extent of enhancement of CBP (δ ECP ). As the extent of enhancement of CBP increased, the percentage of surviving SNc cells increased in a sigmoidal fashion (Figure 6B).
At lower values of δ ECP , the percentage of surviving SNc cells was below 50%. At δ CCB values of > 0.7, the percentage of surviving SNc cells suddenly increased to > 60%.

Apoptotic Signal Blockers
According to the protocol mentioned in the methods section, the neuroprotective effect of apoptotic signal blockers on the percentage loss of SNc cells was implemented. The apoptotic signal blockers therapy was initiated at 25% of SNc cell loss (indicated by red trace in Figure 6C) for varying extent of apoptotic signal blockers (δ ASB ). As the extent of apoptotic signal blockers increased, the percentage of surviving SNc cells increased in a sigmoidal fashion (Figure 6C). At lower values of δ ASB , the percentage of surviving SNc cells was around 50%. At δ ASB values of > 0.5, the percentage of surviving SNc cells seems to be same as when the therapeutic intervention was initiated. In other words, there was no SNc cell loss after the therapeutic intervention was initiated.

DISCUSSION
The goal of this work was to develop a multi-scale model of excitotoxicity in SNc, which would help us in understanding the mechanism behind neurodegeneration due to excitotoxicity under energy deficiency conditions. The present model was the integration of the comprehensive SNc model, which we had developed earlier ) and the previous model of excitotoxicity (Muddapu et al., 2019). From the simulation results, it suggests that compensatory mechanisms exist that try to maintain calcium homeostasis even under energy deficiency conditions. We also suggested therapeutic intervention, which can be neuroprotective in nature. Moreover, the model also gave us scope to explore more novel approaches to halt or slow down the SNc cell loss in PD. From our previous computational studies (Chander and Chakravarthy, 2012;Chhabria and Chakravarthy, 2016), we were able to show that cortical neurons at lower ATP levels exhibited a bursting type of firing patterns. In the present model, we were able to show that SNc neurons exhibited a bursting type of firing pattern under low glucose and oxygen levels ( Figure 3B). The firing rate of the SNc neuron decreased under energy deficiency as a result of bursting, where bursting can be considered as a compensatory mechanism and also an indicator of energy imbalance (Sverrisdóttir et al., 1998(Sverrisdóttir et al., , 2000de Kloet et al., 2005). The intracellular calcium accumulation was a result of failed efflux of calcium by ATP-dependent calcium pumps due to energy deficiency. In the present model of SNc neurons, we have considered a calcium-induced dopamine release mechanism (Oheim et al., 2006;Lee et al., 2009). Even with higher levels of calcium, extracellular dopamine released by SNc neurons was low as a result of failed packing of dopamine into vesicles due to energy deficiency ( Figure 3D) (Blakely and Edwards, 2012;Hnasko and Edwards, 2012).
The relationship between the percentage of SNc cells under energy deficiency and the percentage of surviving SNc cells is threshold-like, where there was a significant drop in the percentage of surviving SNc cells after more than 40% of SNc cells were under energy-deficient conditions. It suggests that there might be mechanisms which compensate at low energy levels in order to enable normal functioning of the SNc neuron (Navntoft and Dreyer, 2016).
Calcium plays an essential role in the normal functioning of a neuron; any imbalance in calcium homeostasis leads to pathogenesis (Orrenius et al., 2003). In the case of SNc neurons, these compensatory mechanisms may play an important role in maintaining calcium homeostasis. In the present SNc model, calcium fluctuations were monitored by three compensatory mechanisms: excess calcium binding to calciumbinding proteins (calbindin and calmodulin), excess calcium sequestered into ER by active pumps, and excess calcium taken up into MT by mitochondrial calcium uniports. Under energy deficiency, calcium homeostasis was disrupted, resulting in calcium accumulation due to the inactivation of ATP-dependent calcium pumps.
As the first line of defense to maintain calcium homeostasis, excess calcium binds to locally available calcium-binding proteins in the cytoplasm. If calcium-binding proteins could not handle the accumulated calcium, as the next line of defense, excess calcium is sequestered into ER by sarco/endoplasmic reticulum calcium-ATPase. If accumulated calcium exceeds the capacity of ER to sequester it, as a final line of defense, excess calcium is taken up by MT by mitochondrial calcium uniports. If accumulated calcium goes beyond the capacity of these three compensatory mechanisms, calcium builds up in both ER and MT, leading to neurodegeneration by ER-stress and mitochondrial-induced apoptotic mechanisms, respectively (Figure 7) (Bano and Ankarcrona, 2018). Thus, the thresholdlike relationship between the percentage of SNc cells under energy deficiency and the percentage of surviving SNc cells can be explained.

Neuroprotective Strategies
The reason behind the lower percentage of surviving SNc cells at higher values of δ GII might be due to a mechanism known as "weak excitotoxicity" (Albin and Greenamyre, 1992). Due to this mechanism, lower levels of available energy lead to calcium accumulation in SNc neurons, which undergo degeneration due to excitotoxicity. However, calcium accumulation in SNc neurons at lower values of δ GII was low as a result of reduced excitation of STN on SNc even though they were exposed to energy deficiency. Thus, glutamate inhibition therapy was proven to be neuroprotective to SNc neurons under energy deficiency conditions (Chan et al., 2010;Masilamoni et al., 2011;Betts et al., 2012;Ferrigno et al., 2015); moreover, it also seems to halt the SNc cell loss.
At intermediate levels of δ DAA (0.9 < δ DAA ≤ 1.1), the percentage of surviving SNc cells was maximum as a result of the total restoration of the dopaminergic drive to the STN neuronal population. However, beyond 1.1, δ DAA results were in the lower percentage of surviving SNc cells when compared to FIGURE 7 | Proposed putative mechanism of excitotoxicity in SNc. Normal-healthy SNc neuron (normal calcium homeostasis); Vulnerable-vulnerable SNc neuron (compromise calcium homeostasis); First line-First line of defense to counter calcium imbalance by CBP (excess calcium binds to CBP); Second line-Second line of defense to counter calcium imbalance by ER (taken up into ER by SERCA); Third line-Third line of defense to counter calcium imbalance by MT (taken up into MT by MCUs); Apoptosis-ER-stress and mitochondrial-induced apoptotic mechanism activated (Calcium-induced neurodegeneration). SNc, substantia nigra pars compacta; ER, endoplasmic reticulum; MT, mitochondria; Ca 2+ , calcium; Cytc, cytochrome c; GLC, glucose; O 2 , oxygen; PYR, pyruvate; ATP, adenosine triphosphate; CBP, calcium-binding proteins; SERCA, sarco/endoplasmic reticulum Ca 2+ -ATPase; MCU, mitochondrial calcium uniports.
intermediate levels. The reason behind the lower percentage of surviving SNc cells at higher values of δ DAA might be due to the excessive activation of autoreceptors on the SNc neurons. There are four types of D2 autoreceptors on the SNc neurons, where they regulate neuronal activity and control dopamine synthesis, release, and uptake (Ford, 2014). When these autoreceptors get activated, they result in reduced neuronal activity, dopamine synthesis, release, and uptake. In the present SNc model, we have considered all autoreceptors except one regulating dopamine uptake. The noisy responses for dopamine restoration therapy was due to activation of the various autoreceptors. At higher levels of dopamine, autoreceptors on the SNc neurons get overactivated, resulting in an overall decrease in dopamine release, which in turn reduced dopaminergic tone to STN neuronal population (Schapira and Olanow, 2003;Piccini and Pavese, 2006;Vaarmann et al., 2013). Decreased dopaminergic tone results in disinhibition of STN, which in turn leads to excitotoxic damage to SNc neurons (Rodriguez et al., 1998). Thus, dopamine restoration therapy proves to be neuroprotective to SNc neurons under energy deficiency. However, it can deteriorate the survivability of SNc neurons under high dosage (Lipski et al., 2011;Paoletti et al., 2019).
The relationship between the percentage of surviving SNc cells and δ LES is threshold-like, where there was a significant rise in the percentage of surviving SNc cells with a δ LES value of just 20%. The reason behind the lower percentage of surviving SNc cells at lower values of δ LES might be due to a mechanism known as "strong excitotoxicity" (Albin and Greenamyre, 1992). Due to this mechanism, overactivation of glutamatergic receptors leads to calcium accumulation in SNc neurons, which undergo degeneration due to excitotoxicity. It suggests that overexcitation from disinhibited STN is indeed a significant contributor to excitotoxic loss of SNc neurons ( Figure 5C). However, a δ LES value higher than 20%, the percentage of surviving SNc cells increased steadily with an increase in δ LES values. At this point, excitotoxic loss of SNc neurons is due to "weak excitotoxicity, " which is overcome by reducing excitatory drive from STN by subthalamotomy. Thus, subthalamotomy therapy proves to be neuroprotective to SNc neurons under energy deficiency, and in addition, it showed a substantial result with lower values of δ LES .
At lower values of δ ASF , the percentage of surviving SNc cells was below 50%. The reason behind the lower percentage of surviving SNc cells at lower values of δ ASF might be due to 'weak excitotoxicity' (Albin and Greenamyre, 1992). As δ ASF increased, the percentage of surviving SNc cells increased as a result of reduced excitation from disinhibited STN by DBS. We hypothesized that the neuroprotective effect of DBS therapy was a result of STN pathological oscillations being blocked from propagating to other nuclei; in other words increased axonal and synaptic failures of STN neurons due to DBS disrupts the information transfer through the stimulation site (Ledonne et al., 2012;Rosenbaum et al., 2014;Muddapu et al., 2019). Thus, DBS therapy appears to be neuroprotective to SNc neurons under energy deficiency, according to the proposed mechanism.
The relationship between the percentage of surviving SNc cells and δ CCB , δ ECP , and δ ASB is threshold-like, where there was a sudden rise in the percentage of surviving SNc cells at higher values (Figure 6). The reason behind the lower percentage of surviving SNc cells at lower values of δ CCB , δ ECP , and δ ASB might be due to "strong excitotoxicity" (Albin and Greenamyre, 1992). It suggests that overexcitation from disinhibited STN caused calcium accumulation in SNc neurons, which was overcome by reducing the influx of calcium into SNc neurons by blocking calcium channels at higher values of δ CCB (Figure 6A). Similarly, excess calcium accumulation in SNc neurons was overcome by enhancing CBP, which binds to free excess calcium and inactivated them at higher values of δ ECP (Figure 6B). It also suggests the neuroprotective effect of apoptotic signal blocker therapy at higher values of δ ASB is a result of blockage of proapoptotic players (such as caspase 9, caspase 3) in the SNc neuron which was activated due to excess calcium accumulation through ER-stress and mitochondrial-induced apoptotic mechanisms. Thus, calcium channel blocker, enhancement of CBP, and apoptotic signal blocker therapies have been proven to be neuroprotective to SNc neurons under energy deficiency. Moreover, apoptotic signal blocker therapy seems to halt the SNc cell loss.
In addition to the above-mentioned therapeutic approaches, we have also simulated other therapeutic interventions such as the mitochondrial calcium channel and permeability transition pore blockers, ATP supplements (phosphocreatine, creatine), and antioxidants which did not show neuroprotective effects in the present model. This could be due to the fact that direct ATP supplements and antioxidants might not be enough to show the neuroprotective effect in the present model.

Limitations and Future Directions
Though the proposed model captured the exciting results of excitotoxicity, it is not without limitations. In the proposed model, the ischemic condition was implemented by modulating glucose and oxygen levels, which can be extended by adding a blood vessel module (Cloutier et al., 2009) and varying cerebral blood flow to simulate the ischemia condition more realistically. In the proposed model, stress was monitored in SNc neurons alone, which can be extended to other neuronal types in the model by monitoring stress levels, where intracellular calcium accumulation can be a stress indicator (Bano and Ankarcrona, 2018). In order to do so, all neuronal types should be modeled as conductance-based models where calcium dynamics should be considered. The synaptic weights in the proposed model were not dynamic; we would like to incorporate learning principles such as spike time-dependent plasticity in the STN population, which can show the long-term effect of DBS treatment (Ebert et al., 2014;Iakymchuk et al., 2015). A more ambitious goal is to incorporate the present model into a large-scale model of basal ganglia (Muralidharan et al., 2018) in understanding the effect of therapeutics in terms of the behavioral response (Erdi et al., 2006). The future version of the model should be able to show the neuroprotective effect of other therapeutic interventions such as ATP supplements (phosphocreatine, creatine) and antioxidants.

Conclusions
In conclusion, we believe that the proposed model provided significant insight into understanding the mechanism behind excitotoxicity in SNc neurons under energy deficiency conditions. From simulation results, it is evident that energy deficiencies occurring at cellular and network levels could precipitate the excitotoxic loss of SNc neurons in PD. From neuroprotective studies, it was clear that glutamate inhibition and apoptotic signal blocker therapies were able to halt the progression of SNc cell loss, a definitive improvement over other therapeutic interventions that only slowed down the progression of SNc cell loss. Logically connected with our earlier studies (Muddapu et al., 2019), we are presently reinforcing the idea of regarding PD as a metabolic disorder where metabolic deficiencies occurring at different levels in the neural hierarchy-subcellular, cellular, and network levels-precipitate the pathological changes of PD at multiple levels. Furthermore, using such models, we hope to be able to design and develop better therapeutics which target the root cause rather than the multitude of effects.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
VM and VC conceived, developed the model, and prepared the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to acknowledge Shreya Jha and Pradeep Varathan for their contributions in developing the initial model. The pre-print version of the manuscript was released by Muddapu and Chakravarthy (2020a) bioRxiv for feedback from fellow neuroscientists.