EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5562 ISSN 1307-5543 – ejpam.com Published by New York Business Global Stability Analysis of Parkinson’s Disease Model with Multiple Delay Differential Equations using Laplace Transform Method G. Veerabathiran1, G. Jagan Kumar1, Siriluk Donganont2,∗ 1 Department of Mathematics, Hindustan Institute of Technology and Science, Chennai, Tamil Nadu, India. 2 School of Science, University of Phayao, Phayao 56000, Thailand. Abstract. In this manuscript, we investigate the stability analysis of Parkinson’s disease model multiple delay differential equations (DDEs) utilizing the Laplace transform method. Delay differ- ential equations are often encountered in a wide range of scientific and engineering applications, such as signal processing, control systems, and population dynamics. These equations are formu- lated using delayed arguments. The analysis and solution of these equations are frequently made more difficult by the existence of delays. Here, we simplify the process of locating explicit solutions by converting the DDEs into algebraic equations using the Laplace transform technique. The sta- bility characteristics of the solutions to the DDE are crucial for understanding the progression of Parkinson’s disease and the effectiveness of treatment strategies. Stable and asymptotically stable solutions are associated with better control and management of the disease, while instability sug- gests a potential for rapid deterioration, requiring more intensive intervention. Understanding the impact of different delays τ1 and τ2 and coefficients on the stability can help in designing better therapeutic protocols, and potentially developing new treatments that target the specific dynamics of the disease as modeled by the DDE. According to our findings, the Laplace transform method offers a methodical and effective way to solve complicated delay differential equations, revealing important information and having potential uses in both theoretical and practical fields. Addi- tionally, a spatiotemporal model of dopamine concentration in Parkinson’s disease is developed using the Laplace transform method, demonstrating its potential to predict symptom fluctuations, treatment strategies, and improve disease understanding, ultimately enhancing patient manage- ment and quality of life. 2020 Mathematics Subject Classifications: 34K20, 44A10 Key Words and Phrases: Delay differential equation (DDE), Parkinson’s disease model, Laplace transform method, Stability ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5562 Email addresses: gvmathveera@gmail.com (G. Veerabathiran), ultimateg1990@gmail.com (G. Jagan kumar), siriluk.pa@up.ac.th (Siriluk Donganont) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 2 of 20 1. Introduction Delay Differential Equations (DDEs) are integral in modeling various dynamic systems across disciplines such as physics, engineering, biology, and economics. For solving DDEs, the Laplace transform method is an effective tool due to its ability to handle initial value problems and incorporate delay terms effectively. Erkan Cimen [1] discussed the solution of second order DDE using the Laplace transform method. Their work emphasizes the application of this method to derive solutions and validate them with examples, demon- strating the theoretical robustness of their approach. Michal Pospisil [2] expand on the Laplace transform by deriving closed - form solutions for systems of non homogeneous linear differential equations with multiple constant delays. They use unilateral Laplace transforms to unity recent findings, providing a comprehensive approach to solving these complex systems. Reem Alrebdi [3] focus on the Pantograph DDE, a fundamental model in delay differential equations. Daniela Marian [4] studies the Hyers - Ulam stability of DDEs using the Laplace trans- form. Alfred Daci [5] and Gilbert Kerr [6] introduce the Laplace transform in solving dif- ference and differential - difference equations, demonstrating its efficiency and speed. The fundamental concepts and properties of the Laplace transform, show casing its practical applications in solving complex equations. Michelle Sherman [7] propose a methodology for solving DDEs with Dirac delta function using the Laplace transform. Tamas Kalmar - Nagy [8] demonstrates the use of the method of steps combined with the inverse Laplace transform for stability analysis of DDE. H.N. Agiza [9] apply DDEs to model Parkinson’s disease. They transform these models using the Taylor series and validate stability conditions using Matlab, highlighting the biological relevance and practical applications of DDEs in medical research. Michelle Sherman [10] compare the performance of Maple and Matlab in solving linear DDEs using the Laplace transform method. They analyze computational time and accuracy, finding that Matlab is generally faster for linear non - neutral DDEs, while Maple performs better for more complex neutral DDEs. Andre A. Kelle [11] explores the application of DDEs in economic macro dynamics. By solving these equations using the Laplace transform and other numerical techniques, the study demonstrates the relevance of DDEs in modeling economic systems with constant or flexible lags [12–15]. Sohaly and Elfouly [16] explore the stability of PD models formulated as nonlinear delay differential equations. Their work highlights that aft er prolonged use of dopamine - enhancing drugs, positive feedback mechanisms may exacerbate patient tremors, leading to instability in the system. The oscillatory behavior of PD models with discrete and distributed delays was investigated by Chunhua Feng [17]. By converting distributed delays into equivalent discrete delays, the model can be linearized around equilibrium points. Feng demonstrates that the stability of the linearized system is a good indicator of the overall system’s stability. Importantly, this work shows that even small delays can destabilize the system, which is crucial for understanding the effects of delayed feedback in PD dynamics. Kayelvizhi and Pushpam [18] introduce a polynomial collocation method based on G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 3 of 20 successive integration techniques for solving nonlinear DDEs in PD models. They explore various classical orthogonal polynomials, such as Bernoulli, Chebyshev, and Hermite poly- nomials, to approximate solutions to the DDEs. Their numerical simulations show that the proposed method is highly effective and reliable compared to traditional step methods. The novelty of their approach lies in its simplicity and applicability to real - world prob- lems across various scientific and engineering domains. An interdisciplinary approach by combining control theory, computational neuroscience, and deep brain stimulation (DBS) for PD treatment examined by Schiff [19]. The study emphasizes the role of modern model - based control theory in improving the efficacy of nonlinear dynamic systems, particularly through the integration of computational models of neuronal networks. Dovzhenok and Rubchinsky [20] delve into the origin of tremors in PD by examining the basal ganglia - thalamo - cortical loop as a primary generator of tremor activity. Using a conductance - based model of subthalamo - pallidal circuits, they demonstrate how variations in dopamine - modulated connections within this loop lead to tremor - like burst firing. Their findings suggest that modulating these connections through dopaminergic therapy or disrupting the loop through surgical interventions can suppress tremor activity. This work provides mechanistic insights into potential therapeutic targets within the basal ganglia - thalamo - cortical loop for tremor suppression. The significance of beta band oscillations, which are associated with motor symptoms in PD, particularly when the patient is OFF medication examined by Duchet [21]. Their study focuses on the duration of beta bursts, analyzing how these bursts change between ON and OFF medication states. By investigating local field potentials from the subtha- lamic nucleus, they show that the non - linearity in beta oscillation dynamics correlates with motor impairments. Langston [22] provided a historical perspective on the discovery of MPTP, a compound that induces selective degeneration of the substantia nigra, mirroring PD’s effects. The discovery of MPTP revolutionized PD research, offering an animal model that mimics the disease’s progression. Lainscsek [23] presented an innovative approach to assessing the severity of PD motor symptoms, particularly finger - tapping movements, through the use of nonlinear delay differential equations. By fitting these equations to time series data from patients, the authors developed a six - dimensional numerical descriptor to rate motor symptoms algorithmically. The intersection of insulin resistance and Parkinson’s disease through a biochemical systems theory (BST) model was explored by Braatz and Coleman [24]. Their work high- lights the complex interactions between insulin signaling pathways and neurodegeneration, with implications for identifying effective treatment strategies. By using Matlab for model simulation, the authors provide insights into how delayed treatments impact disease pro- gression. This model underscores the importance of early diagnosis and offers a framework for testing potential treatments in a computational environment, further illustrating the application of mathematical models in PD research. Bocharov and Rihan [25] emphasize the importance of DDEs in biosciences, particu- larly for modeling phenomena in fields such as epidemiology, physiology, and neuroscience. Their review highlights how DDEs offer a richer mathematical framework compared to or- G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 4 of 20 dinary differential equations (ODEs), better capturing the temporal dynamics of biological processes. They discuss various numerical techniques for solving DDEs, which are essential for understanding complex biosystems. The application of these methods in PD research has contributed to a deeper understanding of the disease’s progression and potential in- terventions. A competition model of tumor growth incorporating the immune response and phase - specific drugs introduced by Villasana and Radunskaya [26]. DDEs model the phases ofthe cell cycle, and stability analysis reveals that the system can exhibit periodic solu- tions through Hopf bifurcations. This highlights the importance of considering delays in understanding tumor - immune interactions and treatment outcomes. Nelson and Perel- son [27] examined HIV - I infection models with intracellular delays. Their findings show that incorporating delays alters the kinetic parameters of the system, particularly the loss rate of infected T cells. The study provides a mathematical framework to understand the implications of imperfect drug efficacy and offers general stability results for nonlinear DDE infection models. Several studies apply DDEs to understand the dynamics of Parkinson’s disease. Lain- scsek [28] use nonlinear delay differential equations to rate finger - tapping movements in Parkinson’s patients, providing a more objective assessment of disease progression. Addi- tionally, Bocharov and Rihan [29] explore the use of DDEs in biosciences, including models of neural networks, which have implications for understanding diseases like Parkinson’s. Braatz and Coleman [30] propose a mathematical model for insulin resistance in Parkin- son’s disease, using Biochemical Systems Theory (BST) and Matlab for more flexible simulations. Their model emphasizes the importance of early diagnosis and treatment. Torelli [31] addresses the stability of numerical methods for DDEs, focusing on the backward Euler method. This work is critical in ensuring that numerical simulations of DDEs yield reliable results. Gopalsamy and Zhang [32] study the stability of impulsively perturbed DDEs, providing sufficient conditions for asymptotic stability and oscillatory behavior. This research extends the understanding of nonlinear dynamics in systems with both delays and impulses. Wolfrum [33] explore DDEs with large delays, introducing novel approaches to stability analysis. Their work on strong and weak instabilities offers insights into the complex dynamics of delay systems, including control systems and semiconductor lasers. Lin and Wang [34] and Li [35] investigate systems with multiple discrete delays. These studies contribute to understanding how multiple delays interact to influence the stability and bifurcation of solutions, with applications ranging from motor control to biological systems. Yan and Zhao [36] discuss the oscillation and stability of linear impulsive DDEs, finding that these properties are equivalent to those of corresponding non - impulsive DDEs. Their results can enhance the stability analysis of systems with sudden changes. Enright and Hayashi [37] develop a solver for neutral DDEs based on a continuous Runge - Kutta method with defect control. Their DDVERK algorithm ensures accurate error and step size control, offering a robust tool for solving DDEs numerically. Keane [38] review the use of DDEs in climate models, focusing on delayed feedback loops in global energy balance and the El Niño Southern Oscillation (ENSO) system. Their study G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 5 of 20 demonstrates that DDEs can capture complex climate dynamics with relatively simple models, offering insights into the predictability of climate phenomena. The reviewed literature underscores the versatility and efficacy of the Laplace trans- form method in solving a wide range of delay differential equations. From theoretical advancements and stability analysis to practical applications in biology, medicine, and economics, the Laplace transform continues to be a critical tool in the mathematical tool kit for addressing dynamic systems with delays. The comparative studies on computa- tional tools further highlight the importance of choosing appropriate soft ware for solving specific types of DDEs, ensuring accuracy and efficiency in obtaining solutions. In this paper, we investigate the analytical solution of the multi - DDE ψ′(t)− aψ(t)− bψ(t− τ1)− cψ(t− τ2) = g(t), t > 0, (1) ψ(t) = ϕ(t), t ≤ 0 (2) and using Laplace transform method to solve Parkinson’s disease model DDEs. 2. Parkinson’s disease model Neurodegenerative disorder characterized by motor symptoms such as rigidity, tremor, and brady kinesia (slowness of movement) is known as Parkinson’s disease (PD). One notable impairment in PD is the disruption of the temporal structure of hand movements. However, the specifics of these temporal distortions and the impact of dopamine replace- ment therapy on them remain largely unexplored. Claudia Lainscsek [39] investigated the use of nonlinear delay differential equations (DDEs) to analyze the repetitive finger tap- ping movements in patients with PD. The study aims to distinguish the spatiotemporal distortions in these movements among PD patients off and on dopamine replacement ther- apy and compares these findings to age-matched control subjects. This method is used to understand and model the temporal evolution of repetitive hand movements in patients with PD. To address this gap, the authors applied nonlinear time series analysis techniques, particularly focusing on DDEs, to examine the finger tapping movements of PD patients. This analysis aimed to uncover the underlying dynamical system’s spectral and topological properties, which are robust against noise. The methodology involved having subjects perform a finger tapping task, where they tapped their right index finger and thumb together rapidly for 10 seconds, repeated three times. The movements were captured using a 12 - camera phase space 3D motion mon- itoring system, providing high - frequency data at 120 Hz. This data was then analyzed using DDEs to model the temporal evolution of the movements. In the Parkinson’s disease model, state variables often represent quantities central to the disease’s progression and response to medication. Commonly used state variables might include: • Dopamine Level: This variable indicates dopamine concentration, which is crucial for motor control and is depleted in Parkinson’s disease. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 6 of 20 • Motor Symptoms: This variable represents the severity of symptoms, such as tremor, rigidity, or bradykinesia. • Medication Effectiveness: This measures how effective the current medication is at managing symptoms. It can vary over time due to dosing schedules and fluctuating responses, known as “on” an “off” periods. The delay differential equation of Parkinson’s disease is characterized by ψ̇ = c1ψ(t) + c2ψτ1 + c3ψτ2 , (3) where ψ(t) is the state variable that represents the system’s behavior at time t, c1, c2 and c3 are coefficients, and delays τ1 and τ2 represent the time lags in the system that affect how past states of the system influence its current state. The state variable ψ(t) represents the state of Parkinson’s disease and may be inter- preted as a measure of symptoms, the level of a key neurochemical (like dopamine), or a composite indicator of disease severity. By modeling ψ(t), we can track how symptoms evolve over time and potentially predict the efficacy of different treatment strategies. The coefficients c1, c2, and c3 describe the influence of the current state c1ψ(t) and the delayed states c2ψ(t− τ1) and c3ψ(t− τ2) on the rate of change of ψ(t) . These coef- ficients could reflect various physiological or pharmacological factors ψ(t) may represent the immediate rate of change due to the current state, capturing the direct, instantaneous effects on the system, while c2 and c3 represent the delayed responses. These delays could arise from physiological lags in response to stimuli, such as the body’s delayed reaction to medication or time - lagged feedback effects in the brain’s neurochemical systems. The delays τ1 and τ2 could correspond to different response times related to disease progression or the effects of medication. For example, the first delay τ1 might represent the lag in the body’s response to dopaminergic medication often observed as “on” periods when symptoms improve, while the second delay τ2 might represent a slower feedback effect or a delayed response to the decline in medication levels leading to “off” periods when symptoms worsen. The “on” and “off” states in Parkinson’s disease are critical for understanding the practical implications of this model. In the “on” state, medication is effective, and symp- toms are temporarily improved. This effect might be reflected in a specific combination of the coefficients or in an additional term in the model that represents an external effect on ψ(t) when medication is present. Over time, the medication’s effects can diminish, leading to worsening symptoms in the “off” state. This transition fr om “on” to “off” could be represented by a time - dependent change in the coefficients (for example, by decreasing c2, c3) to capture the diminishing efficacy of the medication. The temporal aspect of ψ(t) in Parkinson’s disease model captures how the state of the disease changes over time. This could include fluctuations in symptoms due to medication cycles (“on” and “off” states) or the long - term progression of the disease. The spatial aspect involves extending the model to account for how ψ(t) varies across different regions in space, such as distinct areas of the brain affected by Parkinson’s disease. This G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 7 of 20 spatial modeling could be achieved by defining ψ(t) as ψ(x, t), where ψ represents spatial coordinates. Using only one state variable implies that ψ(x, t) acts as a single indicator of disease progression and symptom severity across both time and space. For example, this single variable might represent the concentration of dopamine or a symptom severity index that changes over time and across different parts of the brain. Spatial interactions could be incorporated by allowing the state variable to depend on neighboring values (e.g., through diffusion terms if using a partial differential equation). Modeling with one state variable allows for a simpler and more interpretable represen- tation, but it also means sacrificing some complexity. By using one variable, we may lose detailed information about the multiple underlying physiological processes, but we gain a unified view of the disease’s impact. When a patient with Parkinson’s Disease is off medication, the delays τ1 and τ2 might represent the time it takes for symptoms to manifest aft er the last dose of medication wears off. These delays could also represent the time required for neurodegenerative processes to exert noticeable effects on motor functions. In this scenario, larger delays might indicate a slower progression of symptoms, as it takes longer for the disease to impact the patient significantly aft er the medication effect dissipates. Conversely, shorter delays could indicate a more rapid onset of symptoms once medication is stopped, leading to faster deterioration in motor control and other functions. When the patient is on medication, the delays τ1 and τ2 could reflect the time it takes for the medication to take effect and how long the benefits of the medication persist in the system. The delays might also capture the body’s response time to the medication. Shorter delays might suggest that the medication quickly stabilizes the patient’s symptoms, resulting in more immediate relief. However, if these delays are too short, it might also indicate a rapid onset of medication wearing off, requiring frequent dosing to maintain symptom control. Longer delays, on the other hand, might mean that the medication has a sustained effect, allowing for more extended periods of symptom relief but potentially leading to slower responses to dose adjustments. For healthy controls, the delays τ1 and τ2 might be negligible or significantly different from those in PD patients. These delays could reflect normal physiological time lags in motor control processes or other neurological functions that do not result in symptoms. The difference in delays between healthy controls and PD patients (on or off medication) can be used to distinguish between normal and pathological states. The delays τ1 and τ2 could potentially serve as bio-markers for the effectiveness of treatment or the stage of the disease. Researchers could use this information to develop new therapeutic strategies aimed to achieve better disease management. Insights into how these delays affect disease dynamics could lead to the development of more sophisticated models that account for individual differences in delay patterns, leading to better predictive models for disease progression and treatment outcomes. In this paper, we investigate the stability analysis of DDE (3) with the condition ψ(t) = ϕ(t), t ≤ 0. (4) G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 8 of 20 We employing the linear delay differential equation (3) to analyze the spatiotemporal pat- terns of repetitive hand movements in PD patients. DDEs are powerful mathematical tools that model the temporal evolution of a system by incorporating delays that corre- spond to past states of the system. This allows for capturing complex dynamics, such as those observed in PD. By focusing on the frequency components of hand movements and their interactions, we aim to distinguish between the movement patterns of PD patients on and off medication, and those of age - matched control subjects. The insights gained from stability analysis contribute to a deeper understanding of the dynamic nature of Parkinson’s disease. This can lead to new research avenues focused on identifying critical factors that influence the stability of the disease and exploring how these factors can be manipulated to improve patient outcomes. By understanding which conditions lead to stability or instability, researchers and clinicians can design more effective therapeutic in- terventions. The analysis involves extracting dominant frequencies and their non - linear couplings, which provide a comprehensive understanding of the underlying motor control disruptions in PD. This mathematical modeling approach not only helps in classifying the severity and treatment response in PD but also offers a potential pathway for developing objective diagnostic and monitoring tools for the disease. 3. Laplace transform method In this section, we provide the general solution of multi - DDE using Laplace transform method. Definition 1. Let ψ(t) be a function of t > 0. Then the Laplace transform of ψ(t) defined by L[ψ(t)] = ∫ ∞ 0 e−stψ(t)dt = Ψ(s), where L is Laplace transform operator, t is a time domain and s is a frequency domain. The general form of DDEs is ψ(t) = g(t, ψ(t), ψ(t− τ1), ψ(t− τ2) ..., ψ(t− τk). Definition 2. Let ψ(t− τ) be a positive real valued function. Then L[ψ(t− τ)] = ∫ ∞ 0 e−stψ(t− τ)dt = ϕ(s) + e−stΨ(s), where ϕ(s) = e−sτ ∫ 0 −τ e−srψ(r)dr. Theorem 1. [40] Let ψ(t), t ∈ (0, t0) be the real valued function. Suppose that ψ(t) is piecewise continuous of exponential order and e−αt|ψ(t)| < M , exists for some constants α,M > 0. Then the Laplace transform L of ψ(t) exists. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 9 of 20 Theorem 2. [15] Let G(t), h(t) ≥ 0 are continuous real valued functions on (0, ∞) . If u(t) ≤ h(t) +G(t) ∫ t 0 l(τ)u(τ)dτ . Then u(t) ≤ h(t) +G(t) ∫ t 0 h(τ)l(τ)e ∫ t τ l(ζ)G(ζ)dζdτ. Theorem 3. Let the DDE (1) satisfies Theorem 1. Then, L is the exact solution of (1) - (2). Proof. Integrating the equation (1) with in the limit (0, t), we have ψ(t)− ϕ(0)− a ∫ t 0 ψ(u)du− b ∫ t 0 ψ(u− τ1)du− c ∫ t 0 ψ(u− τ2)du = ∫ t 0 g(u)du. (5) Replacing integral variable by u− τi = ti, for i = 1, 2, ..., n (n delay terms)∫ t 0 ψ(u− τi)du = ∫ 0 −τi ϕ(ti)dti + ∫ t−τi 0 ψ(ti)dti. (6) These expression in equation (5), we can write, ψ(t) + l1(t)− a ∫ t 0 ψ(u)du− b ∫ τ−τ1 0 ψ(t1)dt1 − c ∫ t−τ2 0 ψ(t2)dt2 = ∫ t 0 g(u)du, (7) where l1(t) = −ϕ(0)− b ∫ 0 −τ1 ϕ(t1)dt1 − c ∫ 0 −τ2 ϕ(t2)dt2, |ψ(t)| ≤ |l1(t)|+ l2(t) + |a| ∫ t 0 |ψ(u)|du+ ∫ t 0 |g(u)|du, where l2(t) = |b| ∫ t 0 |ψ(t1)|dt1 + |c| ∫ t 0 |ψ(t2)|dt2. By Theorem 1, ψ(t) is piecewise continuous in the interval (0, t), then e−αu|g(u)| < M, (u > 0), for some constants α,M > 0. Hence, we have ∫ t 0 |g(u)|du ≤ M eαt α . Thus we can write |ψ(t)| ≤ |l1(t)|+ l2(t) + |a| ∫ t 0 |ψ(u)|du+ M eαt α . (8) Using the Gronwall’s inequality, we get |ψ(t)| ≤ C1 +D1 + M eαt α + C2 +D2 + |a|eC1t+D1t ( C1t+D1t+ M α2 ( eαt − 1 )) , G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 10 of 20 where C1, D1, C2, D2 are constants given as C1 = |ϕ(0)|+ |b| ∫ 0 −τ1 |ϕ(t1)|dt1, C2 = |b| ∫ t 0 |ψ(t1)|dt1, D1 = |c| ∫ 0 −τ2 |ϕ(t2)|dt2, D2 = |c| ∫ t 0 |ψ(t2)|dt2. Theorem 4. Let ϕ(t) be continuous on [−∞, 0] and Ψ(s) is the Laplace transformation of ψ(t) . Then, the exact solution of (1) - (2) is ψ(t) = L−1 [ G(s) + T (s) K(s) ] , where T (s) = ϕ(0) + be−sτ1 ∫ 0 −τ1 e−sτ1ψ(t1)dt1 + ce−sτ2 ∫ 0 −τ2 e−sτ2ψ(t2)dt2, K(s) = s− a− be−sτ1 − ce−sτ2 . Proof. From (1) - (2) by using Laplace transform method L[ψ′(t)]− aL[ψ(t)]− bL[ψ(t− τ1)]− cL[ψ(t− τ2)], (9) L[ψ′(t)] = sΨ(s)− ψ(0), L[ψ(t)] = Ψ(s), L[g(t)] = G(s), L[ψ(t− τi)] = ∫ ∞ 0 e−stψ(t− τi)dt. Replacing integral variable by t− τi = ti, for i = 1, 2, . . . , n( n delay terms) L[ψ(t− τi)] = e−sτi ∫ 0 −τi e−sτiψ(ti)dti + e−sτi ∫ ∞ 0 e−sτiψ(ti)dti, L[ψ(t− τi)] = ϕi(s) + e−sτiΨ(s). We can write equation (9), sΨ(s)− ϕ(0)− aΨ(s)− b[ϕ1(s) + e−sτ1Ψ(s)]− c[ϕ2(s) + e−sτ2Ψ(s)] = G(s), [s− a− be−sτ1 − ce−sτ2 ]Ψ(s) = G(s) + T (s), Ψ(s) = [ G(s) + T (s) K(s) ] . G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 11 of 20 Taking Laplace inverse on both sides ψ(t) = L−1 [ G(s) + T (s) K(s) ] . Which is the exact solution of (1) - (2). Theorem 5. Let ϕ(t) be continuous on [−∞, 0] and Ψ(s) is the Laplace transformation of ψ(t). Then, the exact solution of (1) - (2) is ψ(t) = L−1 [ T (s) K(s) ] , (g(t) = 0), where T (s) = 1 + 1 k − s [be−sτ1 + ce−sτ2 ]− 1 k − s [be −kτ1 + ce−kτ2 ], K(s) = s− a− be−sτ1 − ce −sτ2 . Proof. Consider (1) - (2). When t ≤ 0, we have ψ′(t) = aψ(t) + bψ(t) + cψ(t). Then ψ(t) = ekt, k = a+ b+ c. Taking log on both sides, we obtain ψ = ekt+k1 . From (1) - (2) by using Laplace transform method L[ψ′(t)]− aL[ψ(t)]− bL[ψ(t− τ1)]− cL[ψ(t− τ2)] = 0, (10) L[ψ′(t)] = sΨ(s)− ψ(0), L[ψ(t)] = Ψ(s), L[ψ(t− τi)] = ∫ ∞ 0 e−stψ(t− τi)dt. Replacing integral variable by t− τi = ti, for i = 1, 2, . . . , n (n delay terms) L[ψ(t− τi)] = e−sτi ∫ 0 −τi e−sτiψ(ti)dti + e−si ∫ ∞ 0 e−sτiψ(ti)dti, L[ψ(t− τ1)] = ϕ1(s) + e−sτ1Ψ(s), where ϕ1(s) = e−sτ1 ∫ 0 −τ1 e−sτ1ψ(t1)dt1, ϕ1(s) = 1 k − s [e−sτ1 − e−kτ1 ]. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 12 of 20 Similarly L[ψ(t− τ2)] = ϕ2(s) + e−sτ2Ψ(s), where ϕ2(s) = e−sτ2 ∫ 0 −τ2 e−sτ2ψ(t2)dt2, ϕ2(s) = 1 k − s [ e−sτ2 − e−kτ2 ] . We can write equation (10) as follows: sΨ(s)− 1− aΨ(s)− b[ϕ1(s) + e−sτ1Ψ(s)]− c[ϕ2(s) + e−sτ2Ψ(s)] = 0, [s− a− be−sτ1 − ce−sτ2 ]Ψ(s) = T (s), (11) Ψ(s) = T (s) K(s) . Taking Laplace inverse on both sides ψ(t) = L−1 [ T (s) K(s) ] . Which is the exact solution of (1) - (2). 4. Results and Discussion Consider the multi - DDE of Parkinson’s disease model, ψ′(t)− c1ψ(t)− c2ψ(t− τ1)− c3ψ(t− τ2) = 0. (12) Applying Laplace Transform and using theorem (3), we have Ψ(s) = ψ(0) s− c1 − c2e−sτ1 − c3e−sτ2 . (13) Applying inverse Laplace transform, we obtain the solution ψ(t) = L−1 [ ψ(0) s− c1 − c2e−sτ1 − c3e−sτ2 ] . Case 1: Let C1 = −0.6, C2 = 0.4, C3 = 0.3, ψ(0) = 0.5, τ1 = 0.1, τ2 = 0.2. Then the solution ψ(t) converges to 0 as t tends to infinity. Hence, the DDE is asymptotically stable. Case 2: Let C1 = −0.6, C2 = 0.4, C3 = −0.1, ψ(0) = 0.5, τ1 = 0.1, τ2 = 0.2. Then the solution ψ(t) converges to 0 as t tends to infinity. Hence, the DDE is asymptotically stable. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 13 of 20 Figure 1: Analytical solution of DDE: Asymptotically stable Figure 2: Analytical solution of DDE: Asymptotically stable. From figure 1 and figure 2, the asymptotic stability indicates that not only is the system stable, but it will also eventually return to a steady state after perturbations. Over time, the effects of small disturbances diminish, and the system stabilizes at a particular state. This might represent a more desirable condition in disease management, where the patient’s symptoms gradually stabilize over time, possibly due to effective long-term treatment. It suggests that, even if symptoms worsen temporarily, they will eventually return to a manageable state. Case 3: Let C1 = −0.1, C2 = 0.2, C3 = 0.1, ψ(0) = 0.5, τ1 = 0.1, τ2 = 0.2. Then the solution ψ(t) diverges as t tends to infinity. Hence, the DDE is unstable. From figure 3, if the solution is unstable, small disturbances can lead to significant G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 14 of 20 Figure 3: Analytical solution of DDE: Unstable. deviations in the systems behavior. This means that the disease dynamics are highly sensitive to changes, leading to unpredictable outcomes. In the context of Parkin- son’s disease, this could imply a worsening of symptoms, where the disease state could rapidly deteriorate due to small changes in the patient’s condition or response to medication. This might indicate a need for more careful monitoring and possibly adjustments in treatment. The solution of the DDE does not stable, because the solution depend on the initial guess. If the initial guess is zero, then the solution to the DDE is stable. This implies that the perturbations in the system (such as small deviations in physiological parameters or medication timing) will not lead to large deviations in the state of the system. The model suggests that the disease dynamics remain under control and predictable. This could correspond to a scenario where the symptoms of Parkinson’s disease are well - managed, either through medication or other therapeutic interventions, with the system returning to equilibrium aft er disturbances. Understanding the delays differ between the on - medication, off - medication, and control states can help clinicians tailor treatment plans. If a patient’s delays change significantly when off medication, this might suggest a need for a more continuous or sustained release of medication to prevent symptom onset. Tracking changes in these delays over time could provide insights into the progression of the disease. If delays are becoming shorter or more variable, this could indicate a worsening of the disease or a reduction in the effectiveness of the medication. The ability to quantify and compare these delays allows for more personalized treatment strategies. A patient with particularly short delays when off medication might benefit from a different therapeutic approach than one with longer delays. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 15 of 20 5. Developing the Spatiotemporal model of dopamine concentration in Parkinson’s disease Song and Qu [40] investigate the role of delayed global feedback in the genesis and sta- bility of spatiotemporal excitation patterns within biological excitable media. They focus on excitable media such as cardiac and neural tissues, where spatiotemporal dynamics are influenced by external stimuli, like pacing. Through computational models, they explore how delayed feedback mechanisms contribute to the stability of these patterns, revealing that delays can lead to complex behaviors such as oscillations and pattern destabilization. Long et al. [41] introduced a spatial-temporal delay differential equation model to predict traffic flow more accurately by incorporating delay effects due to vehicle interactions over time and space. Weyhenmeyer et al. [42] apply delay differential analysis to multimodal data to classify Parkinson’s disease. The study leverages DDEs to analyze time-series data from Parkinson’s patients, capturing the delayed responses in motor symptoms char- acteristic of the disease. This approach exemplifies the potential of DDEs in medical applications, especially for conditions where delayed physiological responses are prevalent. In this section, we will explore a spatiotemporal model of dopamine concentration in Parkinson’s disease. This model uses delay differential equations (DDEs) with multiple delays to capture the complex dynamics of dopamine regulation, including delayed feed- back, release, and reuptake. By incorporating these delays, the model provides a more comprehensive understanding of the underlying mechanisms and potential therapeutic tar- gets. The goal is to simulate and understand dopamine leivels in the brain, specifically in the context of Parkinson’s disease, a neurodegenerative disorder characterized by the loss of dopamine-producing neurons. In Parkinson’s disease, dopamine levels in different parts of the brain (particularly the substantia nigra and striatum) fluctuate due to the degeneration of dopamine - producing neurons. By letting ψ(x, t) represent dopamine concentration across various spatial loca- tions χ in the brain, this model can capture the spatial spread and temporal evolution of dopamine loss. Understanding how dopamine depletion varies over time and space can aid in predicting symptom onset and severity in different brain regions. It may also guide targeted therapies, such as localized deep brain stimulation (DBS) or focused drug delivery, to specific areas exhibiting the highest rates of dopamine loss. Consider the dopamine concentration model with two delay differential equations [40– 42] ψ′(t) = −0.5ψ(t) + 0.3ψ(t− 5) + 0.7ψ(t− 2) (14) with the initial condition ψ(0) = 1 (baseline dopamine level). Here, the state variable ψ(t) represents dopamine concentration in a specific region of the brain, such as the striatum, which is critically affected in Parkinson’s disease. The con- stants c1, c2, and c3 describe how current and delayed concentrations of dopamine influence its rate of change. τ1 is the delay associated with the brain’s slower dopamine production in response to neuron loss, while τ2 represents the medication response delay, indicating that the impact of medication on dopamine levels is delayed aft er administration. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 16 of 20 Assuming c1 = −0.5 represents a decay factor due to natural dopamine degradation, c2 = 0.3 captures the brain’s attempt to restore dopamine over time, albeit with a delayed effect (τ1) due to neuron degradation. Additionally, c3 = 0.7 represents the effect of dopamine - boosting medication, which also has a delay of τ2. Here, τ1 = 5 hours reflects the delay in the natural response, while τ2 = 2 hours indicates the medication response delay. To solve the model, we apply the Laplace transform to (14), transforming it into an algebraic equation that can be more easily manipulated. The Laplace transform of the function ψ(s, t) is defined as: Ψ(s) = ψ(0) s+ 0.5− 0.3e−sτ1 − 0.7e−sτ2 . Applying the Laplace transform allows us to express the solution in the s - domain, where we can solve for ψ(s, t) and subsequently perform the inverse Laplace transform to obtain ψ(t) = L−1 [ ψ(0) s+ 0.5− 0.3e−sτ1 − 0.7e−sτ2 ] . Once we have derived the solution in the s - domain, we can apply binomial theorem, such as the inverse Laplace transform, to obtain the time - domain solution ψ(t). Figure 4: Spatiotemporal model of dopamine concentration in Parkinson’s disease From Figure 4, increasing dopamine concentration can have several significant im- pacts, particularly in the context of Parkinson’s disease and other neurological conditions. Increased dopamine levels can alleviate motor symptoms such as tremors, rigidity, and bradykinesia (slowness of movement) in Parkinson’s disease. Higher dopamine concen- trations can also improve coordination and balance, which are oft en impaired in these patients. Furthermore, increased dopamine can enhance motivation and the ability to experience pleasure, thereby combating symptoms of depression and apathy that are com- mon in Parkinson’s disease. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 17 of 20 Additionally, elevated dopamine levels can improve the effectiveness of medications used to treat Parkinson’s disease, such as levodopa, which is converted into dopamine in the brain. However, excessively high dopamine concentrations can lead to side effects such as hyperactivity, impulsivity, and psychosis, making careful monitoring essential for patients receiving dopaminergic therapy. Overall, increasing dopamine concentration can significantly enhance the quality of life for individuals with Parkinson’s disease and related disorders. However, careful man- agement is crucial to balance the benefits against potential side effects. Strategies for increasing dopamine levels may include pharmacological interventions, lifestyle changes, and therapeutic approaches aimed at enhancing dopamine function within the brain. 6. Conclusion In this study, we explored the stability analysis of a Parkinson’s disease model governed by multiple delay differential equations (DDEs) using the Laplace transform method. Our approach effectively simplifies the complex nature of DDEs by converting them into al- gebraic equations, allowing for a more accessible path to explicit solutions. The stability characteristics of these solutions are critical in understanding the progression of Parkin- son’s disease and the effect of treatment strategies. Stable and asymptotically stable solu- tions indicate better disease management, while instability points to potential challenges that require more aggressive intervention. The analysis demonstrated that the Laplace transform method offers a systematic and powerful tool for solving complex DDEs, shed- ding light on the role of delays τ1 and τ2 and their impact on the system’s stability. This insight can be invaluable in therapeutic protocols, such as medication timing and dosing, to achieve better outcomes for patients with Parkinson’s disease. Our findings suggest that the method has broad applicability, not only in theoretical explorations but also in practical applications, where understanding the dynamic behaviour of delayed systems is essential. Finally, the impact of the spatiotemporal model of dopamine concentration in Parkinson’s disease using the Laplace transform method was successfully demonstrated, highlighting its potential to predict symptom fluctuations, treatment strategies, and en- hance the understanding of disease progression, ultimately improving patient management and quality of life. Author Contributions The authors equally conceived of the study, participated in its design and coordination, drafted the manuscript, participated in the sequence alignment, and read and approved the final manuscript. Funding This work was supported by the University of Phayao and Thailand Science Research and Innovation Fund (Fundamental Fund 2025, Grant No. 5020/2567). G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 18 of 20 Conflicts of Interest: The authors declare that they have no competing interests. References [1] Erkan Cimen and Sevket Uncu. On the solution of the delay differential equation via laplace transform. Communications in Mathematics and Applications, 11(3):379, 2020. [2] Michal Pospisil and Frantisek Jaros. On the representation of solutions of delayed differential equations via laplace transform. Electronic Journal of Qualitative Theory of Differential Equations, 2016(117):1–13, 2016. [3] Reem Alrebdi and Hind K Al-Jeaid. Accurate solution for the pantograph delay differential equation via laplace transform. Mathematics, 11(9):2031, 2023. [4] Daniela Marian. Laplace transform and semi-hyers–ulam–rassias stability of some delay differential equations. Mathematics, 9(24):3260, 2021. [5] Alfred Daci, Saimir Tola, and Roberto Datja. Laplace transform in difference equa- tions and differential–difference equations. Journal of Multidisciplinary Engineering Science Studies, 8(1):4245–4247, 2022. [6] Gilbert Kerr, Nehemiah Lopez, and Gilberto González-Parra. Analytical solutions of systems of linear delay differential equations by the laplace transform: Featuring limit cycles. Mathematical and Computational Applications, 29(1):11, 2024. [7] Michelle Sherman, Gilbert Kerr, and Gilberto González-Parra. Analytical solutions of linear delay-differential equations with dirac delta function inputs using the laplace transform. Computational and Applied Mathematics, 42(6):268, 2023. [8] Tamás Kalmár-Nagy. Stability analysis of delay-differential equations by the method of steps and inverse laplace transform. Differential Equations and Dynamical Systems, 17:185–200, 2009. [9] HN Agiza, MA Sohaly, and MA Elfouly. Small two-delay differential equations for parkinson’s disease models using taylor series transform. Indian Journal of Physics, 97(1):39–46, 2023. [10] Michelle Sherman, Gilbert Kerr, and Gilberto González-Parra. Comparison of sym- bolic computations for solving linear delay differential equations using the laplace transform method. Mathematical and Computational Applications, 27(5):81, 2022. [11] Andre Keller. Contribution of the delay differential equations to the complex economic macrodynamics. WSEAS Transactions on Systems, 9, 04 2010. [12] Samuel Bernard, Jacques Bélair, Michael C Mackey, et al. Sufficient conditions for stability of linear differential equations with distributed delay. Discrete and Contin- uous Dynamical Systems Series B, 1(2):233–256, 2001. [13] Sun Yi, A Galip Ulsoy, and Patrick W Nelson. Solution of systems of linear delay differential equations via laplace transformation. In Proceedings of the 45th IEEE Conference on Decision and Control, pages 2535–2540. IEEE, 2006. [14] D. S. Mitrinovic, J. E. Pecaric, and A. M. Fink. Gronwall Inequalities on Other Spaces: Discrete, Functional and Abstract, pages 433–467. Springer Netherlands, Dordrecht, 1991. G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 19 of 20 [15] Mohammed A. Sohaly and M A Elfouly. Stability analysis of two-delay differential equation for parkinson’s disease models with positive feedback. 2021. [16] C Kayelvizhi and AEK Pushpam. Application of polynomial collocation method based on successive integration technique for solving delay differential equation in parkinson’s disease. Indian Journal of Science and Technology, 17(2):112 – 119, 2024. [17] Chunhua Feng. Oscillatory behavior of the solutions for a parkinson’s disease model with discrete and distributed delays. Axioms, 13(2):75, 2024. [18] Steven J Schiff. Towards model-based control of parkinson’s disease. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sci- ences, 368(1918):2269–2308, 2010. [19] Andrey Dovzhenok and Leonid L Rubchinsky. On the origin of tremor in parkinson’s disease. 2012. [20] Benoit Duchet, Filippo Ghezzi, Gihan Weerasinghe, Gerd Tinkhauser, Andrea A Kuhn, Peter Brown, Christian Bick, and Rafal Bogacz. Average beta burst duration profiles provide a signature of dynamical changes between the on and off medication states in parkinson’s disease. PLoS computational biology, 17(7):e1009116, 2021. [21] J William Langston. The mptp story. Journal of Parkinson’s disease, 7(s1):S11–S19, 2017. [22] C Lainscsek, P Rowat, Luis Schettino, D Lee, D Song, C Letellier, and H Poizner. Finger tapping movements of parkinson’s disease patients automatically rated us- ing nonlinear delay differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(1), 2012. [23] Elise M Braatz and Randolph A Coleman. A mathematical model of insulin resistance in parkinson’s disease. Computational Biology and Chemistry, 56:84–97, 2015. [24] Gennadii A Bocharov and Fathalla A Rihan. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics, 125(1-2):183–199, 2000. [25] Minaya Villasana and Ami Radunskaya. A delay differential equation model for tumor growth. Journal of mathematical biology, 47:270–294, 2003. [26] Patrick W Nelson and Alan S Perelson. Mathematical analysis of delay differential equation models of hiv-1 infection. Mathematical biosciences, 179(1):73–94, 2002. [27] C Lainscsek, P Rowat, Luis Schettino, D Lee, D Song, C Letellier, and H Poizner. Finger tapping movements of parkinson’s disease patients automatically rated us- ing nonlinear delay differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(1), 2012. [28] Gennadii A Bocharov and Fathalla A Rihan. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics, 125(1-2):183–199, 2000. [29] Elise M Braatz and Randolph A Coleman. A mathematical model of insulin resistance in parkinson’s disease. Computational Biology and Chemistry, 56:84–97, 2015. [30] Lucio Torelli. Stability of numerical methods for delay differential equations. Journal of Computational and Applied Mathematics, 25(1):15–26, 1989. [31] K Gopalsamy and BG Zhang. On delay differential equations with impulses. Journal G. Veerabathiran, G. J. Kumar, S. Donganont / Eur. J. Pure Appl. Math, 18 (2) (2025), 5562 20 of 20 of Mathematical analysis and Applications, 139(1):110–122, 1989. [32] M Wolfrum, S Yanchuk, P Høvel, and E Schøll. Complex dynamics in delay- differential equations with large delay. The European Physical Journal Special Topics, 191:91–103, 2010. [33] Xihui Lin and Hao Wang. Stability analysis of delay differential equations with two discrete delays. Canadian applied mathematics quarterly, 20(4):519–533, 2012. [34] Xiangao Li, Shigui Ruan, and Junjie Wei. Stability and bifurcation in delay– differential equations with two delays. Journal of Mathematical Analysis and Ap- plications, 236(2):254–280, 1999. [35] Jurang Yan and Aimin Zhao. Oscillation and stability of linear impulsive delay differ- ential equations. Journal of Mathematical Analysis and Applications, 227(1):187–194, 1998. [36] Wayne H Enright and Hiroshi Hayashi. A delay differential equation solver based on a continuous runge–kutta method with defect control. Numerical Algorithms, 16:349–364, 1997. [37] Andrew Keane, Bernd Krauskopf, and Claire M Postlethwaite. Climate models with delay differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Sci- ence, 27(11), 2017. [38] Claudia Lainscsek, Luis Schettino, Peter Rowat, Elke van Erp, David Song, and Howard Poizner. Nonlinear dde analysis of repetitive hand movements in parkin- son’s disease. In Applications of Nonlinear Dynamics: Model and Design of Complex Systems, pages 421–425. Springer, 2009. [39] Joel L Schiff. The Laplace transform: theory and applications. Springer Science & Business Media, 2013. [40] Zhen Song and Zhilin Qu. Delayed global feedback in the genesis and stability of spatiotemporal excitation patterns in paced biological excitable media. PLOS Com- putational Biology, 16(10):e1007931, 2020. [41] Qingqing Long, Zheng Fang, Chen Fang, Chong Chen, Pengfei Wang, and Yuanchun Zhou. Unveiling delay effects in traffic forecasting: A perspective from spatial- temporal delay differential equations. In Proceedings of the ACM on Web Conference 2024, pages 1035–1044, 2024. [42] Jonathan Weyhenmeyer, Manuel E Hernandez, Claudia Lainscsek, Howard Poizner, and Terrence J Sejnowski. Multimodal classification of parkinson’s disease using delay differential analysis. In 2020 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pages 2868–2875. IEEE, 2020.