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

Muddapu, Vignayanandam Ravindernath; Chakravarthy, V. Srinivasa

doi:10.3389/fninf.2020.00034

ORIGINAL RESEARCH article

Front. Neuroinform., 30 September 2020
Volume 14 - 2020 | https://doi.org/10.3389/fninf.2020.00034

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

Vignayanandam Ravindernath Muddapu

V. Srinivasa Chakravarthy^*

Laboratory for Computational Neuroscience, Department of Biotechnology, Bhupat and Jyoti Mehta School of Biosciences, Indian Institute of Technology Madras, Chennai, India

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, 2020a,b). 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 Muddapu and Chakravarthy (2020b), 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.

FIGURE 1

Figure 1. Model architecture of excitotoxicity in SNc. SNc, substantia nigra pars compacta; STN, subthalamic nucleus; GPe, globus pallidus externa; Ca²⁺, calcium; ATP, adenosine triphosphate; Cytc, cytochrome c; PYR, pyruvate; O₂, oxygen; GLC, glucose.

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

\begin{array}{l} C^{x} * \frac{d (v_{i j}^{x})}{d t} = 0.04 {(v_{i j}^{x})}^{2} + 5 v_{i j}^{x} + 140 - u_{i j}^{x} + I_{i j}^{x} + I_{i j}^{s y n} & (1) \end{array}

\begin{array}{l} \frac{d (u_{i j}^{x})}{d t} = a (b v_{i j}^{x} - u_{i j}^{x}) & (2) \end{array}

Resetting:

\begin{array}{l} i f v_{i j}^{x} \geq v_{p e a k}^{x} {\begin{matrix} v_{i j}^{x} \leftarrow c \\ u_{i j}^{x} \leftarrow u_{i j}^{x} + d \end{matrix}} & (3) \end{array}

where, $v_{i j}^{x}$ , and $u_{i j}^{x}$ are the membrane potential, and the membrane recovery variables of a neuronal type x at the location (i, j), respectively, $I_{i j}^{x}$ , and $I_{i j}^{s y n}$ 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_{p e a k}^{x}$ 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 Muddapu and Chakravarthy (2020b). The detailed biophysical model of SNc neuron consisted of 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,

\begin{array}{l} \frac{d (v_{i j}^{S N c})}{d t} = \frac{F * v_{c y t}}{C^{S N c} * A_{p m u}} * [J_{m, N a} + 2 * J_{m, C a} \\ + J_{m, K} + J_{i n p}] & (4) \end{array}

where, F is the Faraday's constant, C^SNc is the SNc membrane capacitance, v_cyt is the cytosolic volume, $A_{p m u}$ 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 Muddapu and Chakravarthy (2020b).

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). The intracellular calcium concentration dynamics ([Ca_i]) of the SNc neuron (Muddapu and Chakravarthy, 2020b) is given by,

\begin{array}{l} \frac{d ([C a_{i}])}{d t} = J_{m, C a} - J_{c a l b} - 4 * J_{c a m} - J_{s e r c a, e r} + J_{c h, e r} \\ + J_{l e a k, e r} - J_{m c u, m t} + J_{o u t, m t} & (5) \end{array}

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 non-specific 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:

\begin{array}{l} τ_{R e c e p} * \frac{d (h_{i j}^{x \to y})}{d t} = - h_{i j}^{x \to y} + S_{i j}^{x} (t) & (6) \end{array}

\begin{array}{l} I_{i j}^{x \to y} (t) = W_{x \to y} * h_{i j}^{x \to y} (t) * (E_{R e c e p} - v_{i j}^{y} (t)) & (7) \end{array}

The N-Methyl-D-aspartic Acid (NMDA) current was regulated by voltage-dependent magnesium channels which were modeled as,

\begin{array}{l} B_{i j} (v_{i j}) = \frac{1}{1 + (\frac{M g^{2 +}}{3.57} * e^{- 0.062 * v_{i j}^{y} (t)})} & (8) \end{array}

where, $h_{i j}^{x \to y}$ 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_{i j}^{x}$ 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_{i j}^{y}$ 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²⁺] 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.

TABLE 1

Table 1. Parameter values used in the proposed HEM.

Lateral Connections

The lateral connections in SNc, STN, and GPe, were modeled as Gaussian neighborhoods (Muddapu et al., 2019),

\begin{array}{l} w_{i j, p q}^{m \to m} = A_{m} * e^{\frac{- d_{i j, p q}^{2}}{R_{m}^{2}}} & (9) \end{array}

\begin{array}{l} d_{i j, p q}^{2} = {(i - p)}^{2} + {(j - q)}^{2} & (10) \end{array}

where, $w_{i j, p q}^{m \to m}$ 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,

\begin{array}{l} H_{\infty} = \frac{1}{1 + e^{(\frac{- (v_{i j}^{x} - θ_{g} - θ_{g}^{H})}{σ_{g}^{H}})}} & (11) \end{array}

\begin{array}{l} \frac{d s_{i j}^{x \to y}}{d t} = α * (1 - s_{i j}^{x \to y}) * H_{\infty} - β * s_{i j}^{x \to y} & (12) \end{array}

\begin{array}{l} I_{i j}^{x \to y} (t) = W_{x \to y} * s_{i j}^{x \to y} * (v_{i j}^{y} (t) - E_{G A B A}) & (13) \end{array}

where, $I_{i j}^{x \to y}$ 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_{i j}^{x}$ , and $v_{i j}^{y}$ 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_{i j}^{x \to y}$ 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, $θ_{g}^{H}$ , $σ_{g}^{H}$ are adopted from Rubin and Terman (2004) and given in Table 2.

TABLE 2

Table 2. Parameter values used for SNc lateral connections.

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,

\begin{array}{l} A^{S T N} = s_{max}^{S T N} * e^{(- c d_{s t n} * D A_{s} (t))} & (14) \end{array}

\begin{array}{l} A^{G P e} = s_{min}^{G P e} * e^{(c d_{g p e} * D A_{s} (t))} & (15) \end{array}

\begin{array}{l} A^{S N c} = s_{min}^{S N c} * e^{(c d_{s n c} * D A_{s} (t))} & (16) \end{array}

where, $s_{max}^{S T N}$ , $s_{min}^{G P e}$ , and $s_{min}^{S N c}$ 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.

TABLE 3

Table 3. Parameter values of dopamine effect on target neuronal populations.

The post-synaptic effects of DA in SNc, STN, and GPe were modeled similar to Muddapu et al. (2019),

\begin{array}{l} W_{x \to y} = (1 - c d 2 * D A_{s} (t)) * w_{x \to y} & (17) \end{array}

where, w_{x → y} is the synaptic weight (STN → GPe, GPe → STN, STN → STN, GPe → GPe, STN → SNc, SNc → SNc), cd2 is the parameter that affects the post-synaptic current, and DA_s(t) is the instantaneous dopamine level.

Total Synaptic Current Received by Each Neuronal Type

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.

\begin{array}{l} I_{i j}^{S N c s y n} = F_{S T N \to S N c} * (I_{i j}^{N M D A \to S N c} + I_{i j}^{A M P A \to S N c}) + I_{i j}^{G A B A l a t} & (18) \end{array}

where, $I_{i j}^{N M D A \to S N c}$ and $I_{i j}^{A M P A \to S N c}$ are the glutamatergic currents from the STN neurons corresponding to NMDA and AMPA receptors activation, respectively; $I_{i j}^{G A B A l a t}$ 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.

\begin{array}{l} I_{i j}^{G P e s y n} = I_{i j}^{N M D A \to G P e} + I_{i j}^{A M P A \to G P e} + I_{i j}^{G A B A l a t} & (19) \end{array}

where, $I_{i j}^{N M D A \to G P e}$ and $I_{i j}^{A M P A \to G P e}$ are the glutamatergic currents from STN neuron considering both NMDA and AMPA receptors activation, respectively; $I_{i j}^{G A B A l a t}$ is the lateral GABAergic current from other GPe neurons.

STN

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

\begin{array}{l} I_{i j}^{S T N s y n} = I_{i j}^{G A B A \to S T N} + I_{i j}^{N M D A l a t} + I_{i j}^{A M P A l a t} & (20) \end{array}

where, $I_{i j}^{G A B A \to S T N}$ is the GABAergic current from GPe neuron; $I_{i j}^{N M D A l a t}$ and $I_{i j}^{A M P A l a t}$ are the lateral glutamatergic currents from other STN neurons considering both NMDA and AMPA receptors activation, respectively.

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, calcium-binding 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,

\begin{array}{l} i f C a_{i j}^{E R} (t) > E R_{t h r e s}, t h e n v_{i j}^{S N c} (t) = 0 & (21) \end{array}

where, $C a_{i j}^{E R}$ is the calcium concentration in the ER, ER_thres is the calcium concentration threshold after which ER-stress induced apoptosis is initiated $(E R_{t h r e s} = 2.15 \times 1 0^{- 3} m M)$ , $v_{i j}^{S N c}$ is the membrane voltage of neuron at the lattice position (i, j).

The ER calcium concentration ([Ca_er]) dynamics is given by,

\begin{array}{l} \frac{d ([C a_{e r}])}{d t} = \frac{β_{e r}}{ρ_{e r}} * (J_{s e r c a, e r} - J_{c h, e r} - J_{l e a k, e r}) & (22) \end{array}

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,

\begin{array}{l} i f C a_{i j}^{M T} (t) > M T_{t h r e s}, t h e n v_{i j}^{S N c} (t) = 0 & (23) \end{array}

where, $C a_{i j}^{M T}$ 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_{i j}^{S N c}$ is the membrane voltage of neuron at the lattice position (i, j).

The MT calcium concentration ([Ca_mt]) dynamics is given by,

\begin{array}{l} \frac{d ([C a_{m t}])}{d t} = \frac{β_{m t}}{ρ_{m t}} * (J_{m c u, m t} - J_{o u t, m t}) & (24) \end{array}

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,

\begin{array}{l} W_{S T N \to S N c} (N_{s c}, t) = {\begin{cases} W_{S T N \to S N c}^{0}, N_{s c} (t) > T_{I} \\ W_{S T N \to S N c}^{0} * δ_{G I I}, N_{s c} (t) \leq T_{I} \end{cases} & (25) \end{array}

\begin{array}{l} T_{l} = P_{s n c} - (p c l * P_{s n c}) & (26) \end{array}

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_{S T N \to S N c}^{0}$ 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,

\begin{array}{l} D A (N_{s c}, t) = {\begin{cases} D A_{s} (t), N_{s c} (t) > T_{I} \\ D A_{s} (t) + δ_{D A A}, N_{s c} (t) \leq T_{I} \end{cases} & (27) \end{array}

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,

\begin{array}{l} i f N_{s c} (t) \leq T_{I}, t h e n v_{i j}^{S T N} (δ_{L E S}, t) = - 80 & (28) \end{array}

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_{D B S} = \frac{1}{T_{D B S}})$ , 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,

\begin{array}{l} P_{B W} (t) = {\begin{cases} A_{D B S}, t_{k} \leq t < t_{k} + \frac{δ_{D B S}}{2} \\ - A_{D B S}, t_{k} + \frac{δ_{D B S}}{2} \leq t < t_{k} + δ_{D B S} \\ 0, e l s e \end{cases} & (29) \end{array}

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,

\begin{array}{l} I_{i j}^{D B S - S T N} (t) = \sum_{β = 1}^{N_{c p}^{F C P}} {M_{β} (t) * P}_{B W} (t) * e^{\frac{- [{(i - i_{c})}^{2} + {(j - j_{c})}^{2}]}{σ_{D B S - S T N}^{2}}} & (30) \end{array}

where, $I_{i j}^{D B S - S T N} (t)$ is the DBS current received by STN neuron at position (i, j) considering the lattice position (i_c, j_c) as the electrode contact point at a time (t), $M_{β} (t)$ is the indicator function which controls the activation of stimulation site β, $N_{c p}^{F C P}$ is the number of activated stimulation contact points $(N_{c p}^{F C P} = 4)$ , P_BW(t) is the biphasic current waveform at time t, and σ_DBS−STN is used to control the spread of stimulus current in the STN network.

In the proposed excitotoxicity model, the DBS therapy was implemented by the following criterion,

\begin{array}{l} I_{i j}^{D B S - S T N} (N_{s c}, t) = {\begin{cases} 0, N_{s c} (t) > T_{I} \\ I_{i j}^{D B S - S T N} (t), N_{s c} (t) \leq T_{I} \end{cases} & (31) \end{array}

where, $I_{i j}^{D B S - S T N} (N_{s c}, 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,

\begin{array}{l} W_{S T N \to S N c} (S_{D B S}, t) = {\begin{cases} W_{S T N \to S N c}, S_{D B S} = O F F \\ W_{S T N \to S N c} * W_{A S F} (δ_{A S F}), S_{D B S} = O N \end{cases} & (32) \end{array}

where, $W_{S T N \to S N c} (S_{D B S}, t)$ is the change in synaptic weight of STN to SNc based $S_{D B S} = {O N, O F F}$ at a particular time (t), $S_{D B S}$ 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.

TABLE 4

Table 4. Parameter values of neuroprotective therapies.

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,

\begin{array}{l} I_{i j}^{C a L} (N_{s c}, t) = {\begin{cases} I_{i j}^{C a L} (t), N_{s c} (t) > T_{I} \\ I_{i j}^{C a L} (t) * δ_{C B}, N_{s c} (t) \leq T_{I} \end{cases} & (33) \end{array}

where, $I_{i j}^{C a L} (N_{s c}, 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 calcium-binding 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,

\begin{array}{l} [C B P_{i j}^{x}] (N_{s c}, t) = {\begin{cases} [C B P_{i j}^{x}] (t), N_{s c} (t) > T_{I} \\ [C B P_{i j}^{x}] (t) + δ_{E C P}, N_{s c} (t) \leq T_{I} \end{cases} & (34) \end{array}

where, $[C B P_{i j}^{x}] (N_{s c}, t)$ is the CBP concentration based on the number of SNc cells surviving at a particular time (t) (N_sc(t)), $[C B P_{i j}^{x}] (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,

\begin{array}{l} [A P O P_{i j}] (N_{s c}, t) = {\begin{cases} [A P O P_{i j}] (t), | N_{s c} (t) > T_{I} \\ [A P O P_{i j}] (t) * δ_{A B}, | N_{s c} (t) \leq T_{I} \end{cases} & (35) \end{array}

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.

FIGURE 2

Figure 2. Characteristic firing patterns of the SNc neuronal model. Oscillations of membrane voltage potential (A–C), intracellular calcium concentration (D–F), and extracellular dopamine concentration (G–I) in tonic spiking, tonic bursting, and phasic bursting. SNc, substantia nigra pars compacta; V, voltage; Ca²⁺, calcium; eDA, extracellular dopamine; mV, millivolt; mM, millimolar; sec, second.

FIGURE 3

Figure 3. Model response of SNc neurons under hypoglycemia and hypoxia conditions. (A) Average frequency of firing. (B) Burst index. (C) Intracellular calcium concentration. (D) Extracellular dopamine concentration. SNc, substantia nigra pars compacta; Ca²⁺, calcium; eDA, extracellular dopamine; DA, dopamine; Hz, hertz; mM, millimolar.

FIGURE 4

Figure 4. Response of the proposed hybrid model of excitotoxicity under energy deficiency. (A) Percentage of surviving SNc cells against the percentage of SNc cells under ED. (B) Percentage of surviving SNc cells for varying percentage of SNc cells under ED and synaptic weight from STN to SNc. STN, subthalamic nucleus; SNc, substantia nigra pars compacta; ED, energy deficiency.

FIGURE 5

Figure 5. Proposed model response for different neuroprotective therapies. Percentage of surviving SNc cells in glutamate inhibition (A), dopamine restoration (B), subthalamotomy (C), and deep brain stimulation (D) therapies. Red trace indicates the percentage loss of SNc cells at which intervention of a particular therapy was initiated. SNc, substantia nigra pars compacta; GII, glutamate inhibition; DAA, dopamine restoration; LES, STN lesioning; ASF, axonal and synaptic failures.

FIGURE 6

Figure 6. Proposed model response for different neuroprotective therapies. Percentage of surviving SNc cells in calcium channel blockers (A), enhancement of calcium-binding proteins (B), and apoptotic signal blockers (C) therapies. Red trace indicates the percentage loss of SNc cells at which intervention of a particular therapy was initiated. SNc, substantia nigra pars compacta; CCB, calcium channel blockers; CBP, calcium-binding proteins; ECP, enhancement of calcium-binding proteins; ASB, apoptotic signal blockers.

Characteristic Firing of SNc Neurons

SNc neurons exhibit two distinct firing patterns: low-frequency irregular tonic or background firing (3 − 8 Hz) and high-frequency 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⁻⁴ mM, and it rose to higher than 1 × 10⁻³ 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⁻⁶ mM, which was in the range of (34 − 48) × 10⁻⁶ 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⁻⁴ mM, and rose up to 1 × 10⁻³ 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⁻⁶ mM (Figure 2H), which was in the range of (7 − 30) × 10⁻⁶ 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⁻⁶ 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⁻⁴ mM, and rose to a higher value than 1 × 10⁻³ mM upon stimulation (Figure 2F). Dopamine released by the SNc neuron during phasic bursting peaked as high as ~ 150 × 10⁻⁶ mM (Figure 2I), which is in the range of (150 − 400) × 10⁻⁶ 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⁻³ mM (Figure 3C). Under energy-deficient conditions, average extracellular dopamine concentration released by the SNc neuron was very low (~ 4 × 10⁻⁶ mM) compared to normal conditions (~ 14 × 10⁻⁶ 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 (0 to 1 with a step size of 0.1). The percentage of surviving SNc cells, at low values of W_{stn → snc} (0.1 and 0.2) showed a linear decrease with an increase in the percentage of SNc cells under energy deficiency (Figure 4B). However, the pattern of the percentage of surviving SNc cells changed to sigmoidal at intermediate values of W_{stn → snc} (0.3 to 0.6) with an increase in the percentage of SNc cells under energy deficiency (Figure 4B). Beyond W_{stn → snc} values of 0.6, the percentage of surviving SNc cells approached zero, with an increase in the percentage of SNc cells under energy deficiency (Figure 4B).