1. Introduction
Parkinson’s disease (PD) is a progressive neurodegenerative disorder characterized by the loss of dopaminergic neurons in the substantia nigra pars compacta and the accumulation of misfolded α-synuclein in Lewy bodies. Mounting evidence supports a prion-like mechanism [1] of α-synuclein propagation between neurons, contributing to regional neurotoxicity. The disease follows a stereotyped anatomical trajectory—initially affecting the peripheral and lower brainstem regions and eventually spreading to cortical areas—consistent with Braak staging [2]. Neuronal vulnerability is influenced by oxidative stress [3], mitochondrial dysfunction [4], and impaired proteostasis [5]. Although significant progress has been made in understanding these factors, many questions remain regarding the spatiotemporal evolution of neurodegeneration.
To address this, we designed a computational model to simulate α-synuclein spread and dopaminergic neuron degeneration. Our framework draws on cellular automata to represent neuron states, misfolded protein burden, and intracellular stress. The model incorporates experimentally informed parameters such as seed-dependent aggregation probability, oxidative stress thresholds, and Braak-stage connectivity. By varying key conditions such as mitochondrial resilience or seeding density, we aim to map out how early disruptions can shape long-term neurodegenerative trajectories. This approach supports hypothesis testing for disease initiation and provides a foundation for modeling therapeutic interventions.
2. Materials and Methods
This study employs a custom Python-based cellular automaton model to simulate the progressive degeneration [6] of dopaminergic neurons within a simplified 2D representation of the substantia nigra. Each grid cell represents a single neuron that can dynamically transition between four states: healthy, stressed, degenerating, or dead. The simulation aims to mimic the spatial and temporal characteristics of Parkinson’s disease pathology, particularly the prion-like spread of misfolded α-synuclein and the compounding effects of oxidative stress.
Each neuron interacts locally with its Moore neighborhood (the eight surrounding cells). At every time step, neurons evaluate their environment and update their state according to a set of biologically inspired rules:
Healthy (H) neurons can uptake misfolded α-synuclein from adjacent degenerating or dead neurons.
Once exposed, the neuron accumulates reactive oxygen species (ROS), modeled as a linear function of α-synuclein burden [7].
When accumulated ROS exceeds an assigned stress threshold, the neuron becomes stressed (S).
If the oxidative burden persists, the neuron progresses to a degenerating (D) state, followed by cell death (X).
Initial α-synuclein seeds are distributed in a 5% ventral-lateral cluster to simulate mid-stage Braak pathology. The misfolded protein propagates over time with a probability P_propagate, and ROS accumulation increases with each contact [8]. Key parameters used in the simulations are as follows:
Grid dimensions: 50 × 50 neurons
Initial α-synuclein seed density: 5% of grid
Propagation probability (P_propagate): 0.3
ROS increase per α-exposure event: 0.05 units
Oxidative stress threshold: Ranged from 0.3 to 0.7 in different runs
Time steps per trial: 500
Replicates per condition: 10 simulations to account for stochastic variation
Throughout each trial, key metrics such as total neuron survival, propagation radius, average time to degeneration, and the rate of spread were tracked. Interventions were modeled by adjusting the oxidative stress threshold mid-simulation to simulate antioxidant therapy, or by reducing P_propagate to mimic inhibition of α-synuclein uptake.
In our simulations, ROS accumulated additively over time and did not decay between time steps. The ROS level at time t + 1 was calculated using the formula:
ROSt+1 = ROSt + 0.05 * Nexposures
Baseline parameter values were chosen based on prior computational and experimental literature: a propagation probability of 0.3 [9] reflects moderate transmission efficiency reported in agent-based models [George & Brundin, 2019], a ROS increment of 0.05 approximates incremental oxidative burden [Cuddy et al., 2019], and a stress threshold of 0.5 balances sensitivity and resilience to ROS insult. The resulting neuronal state distribution and transitions across the grid are visualized in Figure 1.
Figure 1. This heatmap depicts a 50 × 50 neuron grid representing the substantia nigra. Colors correspond to neuronal states: healthy (0, blue), stressed (1, pink), degenerating (2, light red), and dead (3, dark red).
3. Results
The simulation revealed a nonlinear pattern of dopaminergic neuron degeneration driven by local α-synuclein propagation and cumulative oxidative stress. Across all baseline trials (with 5% initial seeding, propagation probability = 0.3, oxidative stress threshold = 0.5) [10], approximately 68.2% ± 4.9% of neurons were dead by step 500. The degeneration began in the ventral-lateral seeding zone and spread concentrically outward, consistent with Braak stage 3 - 5 topography.
3.1. Effect of Oxidative Stress Threshold
Varying the oxidative stress threshold significantly impacted total neuronal loss. At a low threshold of 0.3, neuron death was rapid and widespread, with over 85% of the grid degenerating within the first 300 steps, referencing Figure 2. In contrast, raising the threshold to 0.7 slowed progression considerably, with final degeneration limited to only 41.5% ± 3.1% of the neuron population, as seen in Table 1. These results suggest that neuronal resilience to ROS plays a critical role in shaping the trajectory of disease spread.
Figure 2. This figure shows the accumulation of intracellular oxidative stress in response to misfolded α-synuclein propagation. Lower thresholds lead to rapid ROS buildup, while higher thresholds offer neuroprotection by delaying cellular transition to degenerating states.
Table 1. Effect of oxidative stress threshold on neuron degeneration and time to cell death.
Threshold |
% Dead Neurons (step 500) |
Avg. Time to Death (steps) |
0.3 |
87.4% ± 2.7% |
164.3 ± 12.5 |
0.5 |
68.2% ± 4.9% |
233.6 ± 18.1 |
0.7 |
41.5% ± 3.1% |
305.2 ± 19.6 |
3.2. Impact of α-Synuclein Propagation Rate
Increasing P_propagate from 0.3 to 0.5 accelerated degeneration and widened the radius of spread. At the highest tested propagation rate (P = 0.5), complete grid involvement occurred by step 400 in 8 out of 10 replicates. Conversely, lowering P_propagate to 0.1 preserved spatial restriction, and neuron death was largely confined to the initial seeding region. These findings support the hypothesis that α-synuclein transmissibility significantly governs disease scale and speed.
3.3. Simulated Therapeutic Intervention
When antioxidant intervention was introduced at step 200 by raising the oxidative stress threshold from 0.5 to 0.7, the total number of dead neurons dropped from 68.2% to 52.4% ± 3.8%. Earlier intervention (at step 100) further improved outcomes, limiting degeneration to 39.1% ± 4.3%, seen in Figure 3. However, delayed intervention after step 300 showed minimal impact. These data highlight a critical therapeutic window during early to mid-propagation stages when ROS-buffering therapies may be most effective.
Figure 3. Line plot comparing neuron survival trajectories across baseline and two simulated antioxidant interventions. Early intervention (raising oxidative stress threshold at step 100) preserves viability significantly more than late intervention (raising at step 300).
3.4. Spatial Characteristics of Spread
In all baseline runs, the spread of degeneration maintained radial symmetry with minor irregularities due to random seed placement. Degeneration proceeded outward in concentric zones, matching known histopathological gradients of α-synuclein burden in PD. The final propagation radius under baseline conditions was approximately 21 ± 2 neurons, measured from the seed center.
4. Discussion
The results of this simulation suggest that dopaminergic neuron degeneration in Parkinson’s disease is not simply a linear function of α-synuclein load but a complex interaction between intracellular stress thresholds, local protein propagation, and spatial topology. The simulation reproduced several key features of PD progression described in clinical and pathological studies, including the characteristic ventral-to-dorsal and medial-to-lateral spread, as well as the nonlinear acceleration of neurodegeneration [11] once a critical burden of stress is reached.
A particularly significant observation was the sensitivity of disease spread to both oxidative stress thresholds and propagation probability. These findings support previous experimental studies showing that neurons with lower antioxidant capacity are more prone to degeneration in PD-affected regions. Furthermore, the results support the Braak hypothesis of prion-like α-synuclein transmission, in which localized pathology may expand into broader neurodegeneration.
The model also demonstrated the importance of timing in therapeutic intervention. Simulated antioxidant therapies introduced during early propagation [12] stages significantly reduced total neuronal loss [13], while delayed intervention had minimal effect. This aligns with human trial data showing that many neuroprotective treatments fail if initiated after motor symptom onset, suggesting the need for early biomarker-driven diagnosis.
While the model simplifies several biological features—including glial involvement, immune response, and intracellular degradation pathways—it provides a computationally efficient platform for visualizing disease dynamics and testing new hypotheses. Future versions could incorporate additional cell types and feedback loops to simulate immune activation [14], lysosomal degradation, and circuit-level effects [15].
Our model captures key spatial patterns consistent with postmortem histopathological data, including ventrolateral initiation and radial α-synuclein spread observed in Braak stage 3 - 5 regions [Del Tredici & Braak, 2009]. While the 50 × 50 grid serves as a useful proof-of-concept, it does not account for long-range brain connectome pathways, glial activation [16], or immune signaling. These biological features, along with circuit-level effects, represent important areas for future expansion.
We acknowledge that the current model is limited to local, neuron-to-neuron interactions. Including glial cells and simulating neuroimmune feedback may enhance accuracy. Sensitivity analyses with varied seeding densities, spatial layouts, and larger grids could improve scalability and realism. Each 500-step simulation ran in under two minutes on a 2.4 GHz processor, suggesting feasibility for real-time integration into larger-scale brain models.