A Multiscale, Systems-Level, Neuropharmacological Model of Cortico-Basal Ganglia System for Arm Reaching Under Normal, Parkinsonian, and Levodopa Medication Conditions

In order to understand the link between substantia nigra pars compacta (SNc) cell loss and Parkinson's disease (PD) symptoms, we developed a multiscale computational model that can replicate the symptoms at the behavioural level by incorporating the key cellular and molecular mechanisms underlying PD pathology. There is a modelling tradition that links dopamine to reward and uses reinforcement learning (RL) concepts to model the basal ganglia. In our model, we replace the abstract representations of reward with the realistic variable of extracellular DA released by a network of SNc cells and incorporate it in the RL-based behavioural model, which simulates the arm reaching task. Our results successfully replicated the impact of SNc cell loss and levodopa (L-DOPA) medication on reaching performance. It also shows the side effects of medication, such as wearing off and peak dosage dyskinesias. The model demonstrates how differential dopaminergic axonal degeneration in basal ganglia results in various cardinal symptoms of PD. It was able to predict the optimum L-DOPA medication dosage for varying degrees of cell loss. The proposed model has a potential clinical application where drug dosage can be optimised as per patient characteristics.


INTRODUCTION
Parkinson's disease is the second most prominent neurodegenerative disease after Alzheimer's (Gonzalez-Rodriguez et al., 2020;Marino et al., 2020;Muddapu and Chakravarthy, 2021). The onset of the disease is characterised by shaky movements, the rigidity of joints, unregulated movements, and even loss of smell (Morley and Duda, 2010;Fullard et al., 2017;Armstrong and Okun, 2020;Balestrino and Schapira, 2020;Goldman and Guerra, 2020;Marino et al., 2020). The major cause of Parkinson's disease (PD) is the death of dopaminergic neurons in substantia nigra pars compacta (SNc) (Michel et al., 2016;Surmeier, 2018;Muddapu et al., 2020a). Dopamine (DA) deficiency due to SNc cell loss manifest as the cardinal PD symptoms that include tremor, rigidity, bradykinesia, and postural imbalance (Bereczki, 2010;Poewe et al., 2017;Balestrino and Schapira, 2020). Epidemiological data from the United States alone indicates that there has been an exponential growth of people suffering from PD over the last few decades (Dorsey et al., 2018;Marras et al., 2018). However, the exact cause of this cell death is still not known. Various lines of investigation, experimental and computational, are in progress and hopefully, we will be able to narrow down the roots of this disease (Pissadaki and Bolam, 2013;Pacelli et al., 2015;Fu et al., 2018;Giguère et al., 2019;Muddapu et al., 2019Muddapu et al., , 2020aAnilkumar et al., 2020;Gonzalez-Rodriguez et al., 2020;Muddapu and Chakravarthy, 2021). Understanding the cause and effect relationship between the underlying pathology and symptoms of any neurological disease has fundamental challenges since the roots of the disease are at the molecular and cellular level while the symptoms are seen at the behavioural level (Bakshi et al., 2019). Hence it is important to have a multi-scale model that spans molecular mechanisms to behavioural outputs. With this motivation in mind, we present a computational model that relates DA deficiency in PD to motor symptoms in ON and OFF conditions of medication. As an example of drug action, we simulate the effect of levodopa (L-DOPA) drug administration in our model.
The major ingredients of this computational modelling approach include a behavioural simulation of reaching task that replicates normal and parkinsonian movements, the basal ganglia cortico-motor circuitry that will drive the movements, a dopaminergic subsystem that modulates the control circuitry of the basal ganglia (BG), and pharmacological intervention of L-DOPA medication.

Behavioural Simulation and the Cortico-Basal Ganglia Circuitry
Reaching movements is considered as one of the signatures of planned coordinated movement. Insights into the modelling approaches to these coordinated movements of the arm using the control feedback approach were inspired from previous modelling studies (Fitts, 1954;Morasso, 1981;Knill and Pouget, 2004;Körding and Wolpert, 2004;Todorov, 2004;Shadmehr and Krakauer, 2008). These studies didn't represent the neural correlates underlying the functionality. Later, neural correlates behind coordinated movements of arm (Doya, 1999;Nakahara et al., 2001;Hikosaka et al., 2002) and a reinforcement learning-based two-link arm model with kinetic parameters were developed (Izawa et al., 2004). In this study, we are trying the bridge the molecular level representation of dopamine to the behavioural level representation of motor movements.

Basal Ganglia and Motor Learning
The interactions between the cortex and BG play a very important role in motor learning. It is through these interactions that the decision is made between two competing signals-one favouring the direction of movement and the other suppressing the movement. In order to facilitate this process, an action selection mechanism happens in BG subcortical structure, globus pallidus interna (GPi). The action selection mechanism of BG has been explored by various research groups (Gurney et al., 2001;Humphries et al., 2006;Bogacz and Gurney, 2007). Before the action selection takes place at GPi, the signal is forwarded to the GPi through distinct parallel projections from the striatum facilitated by DA receptor type 1 (D1) and DA receptor type 2 (D2) of medium spiny neurons (MSNs) that are modulated using the dopaminergic input from the SNc (Moustafa et al., 2016;Chakravarthy and Moustafa, 2018). In this study, we explored the role of DA in BG functioning and motor learning.

Dopamine and BG Pathways
Dopaminergic input from the SNc neurons modulates the DA receptors present in the striatal neurons, the input nuclei of the BG differentially. The striatum consists of the D1 and D2 expressing MSNs that project via two different pathways. D1-MSN neurons project along the direct pathway, D2-MSNs project along the indirect pathway. The direct pathway projects directly to the output nuclei, GPi and substantia nigra pars reticulata (SNr), whereas the indirect pathway projects to the output nucleus, GPi, via globus pallidus externa (GPe), and subthalamic nucleus (STN). DA release from SNc neurons maintains the balance between activation of direct and indirect pathways. In order to understand the effect of DA deficiency as in PD conditions or the mechanism of DA replenishment by administration of L-DOPA, we need to understand DA synthesis, uptake, and release, which was explored in this study (Chakravarthy and Moustafa, 2018;Muddapu and Chakravarthy, 2021).

Dopamine Deficiency and L-DOPA Medication
Dopamine (DA) deficiency due to SNc cell loss manifest as the cardinal PD symptoms that include tremor, rigidity, bradykinesia, and postural imbalance (Bereczki, 2010;Poewe et al., 2017;Balestrino and Schapira, 2020). L-DOPA is one of the first-line treatment methodologies for PD. The effect of L-DOPA medication on DA turnover processes in SNc terminal (Best et al., 2009;Reed et al., 2012) and its effects on neural systems of BG was modelled (Baston et al., 2016). However, the effect of L-DOPA medication at the behavioural level has not been explored. In this study, we explored the effect of DA deficiency and L-DOPA intervention on behavioural output.
In this paper, we present a multiscale model of the cortico-BG system to simulate arm reaching movements under normal, parkinsonian, and L-DOPA medication conditions. At the lowest level, the intracellular molecular pathways of SNc cells are modelled so as to capture dopamine synthesis, uptake, and release. At the next level, the BG circuitry is modelled using ratecoded neurons which are cast within the reinforcement learning framework with striatum acting as the neural correlate for critic and the direct and indirect pathways facilitating exploitation and exploration, respectively. At the highest level, arm reaching movements are modelled by a two-link arm model driven by a sensory-motor cortical loop.
This article is organised into multiple sections. Section Materials and methods describes the model architecture, equations, and methods. Here we discuss various functional loops that constitute the model and how they are interconnected. This section also covers the integration of pharmacological intervention. In section Results we showcase the results from the model starting with training the model, simulating the behaviour of a control subject, replicating the PD ON condition and some of the cardinal symptoms, assessing the performance in terms of reaching time, and verifying the effect of L-DOPA therapeutic intervention. The model results also gave an indicator of how to optimise the drug dosage. Section Discussion discusses the simulation results in detail and presents the potential future scope and based on that the conclusion derived is given in section Conclusions.

MATERIALS AND METHODS
The proposed multiscale cortico-basal ganglia (MCBG) model was able to simulate the arm reaching in normal and Parkinsonian conditions which include some of the cardinal symptoms of PD (Figure 1). In addition, the effect of L-DOPA medication on arm reaching in PD condition was simulated (Supplementary Figure 1).
The proposed model can be broadly described in three parts. (i) Outer loop-motor-sensory loop, (ii) Inner loopcortico-basal ganglia loop, and (iii) Central loop-nigrostriatal loop (Figure 2). The outer loop consists of the motor cortex (MC), motor neurons (MNs), arm, proprioceptive cortex (PC), and prefrontal cortex (PFC). The inner loop consists of MC, thalamus, and BG nuclei comprised of the striatum, GPi, GPe, and STN. The central loop consists of striatum and SNc, which plays an important role in simulating PD conditions, where nigrostriatal and nigrosubthalamic pathways are affected by SNc cell loss. For L-DOPA medication, a pharmacokinetic module was formulated where input will be L-DOPA dosage given to the PD patient and output will be the amount of DA released in the striatum during the medication. The subsequent sections describe the dynamics involved in each of these three loops.

Outer Loop-Sensory-Motor Loop
The functional pathway of the outer loop is shown in Figure 2A. The outer loop consists of a two-link arm model driven by MNs. MNs receive motor commands from MC. The end effector position of the arm is sensed by PC and it forwards the signal to MC, which receives signals from PFC and BG. MC issues the motor commands based on the integration of incoming signals.

Arm Model
The kinetic model of the two-joint arm simulates the movement of the arm in two-dimensional space (Izawa et al., 2004;Zadravec and Matjačić, 2013;Supplementary Figure 2). Each joint (shoulder and elbow) is controlled by an agonist (Ag) and antagonist (An) muscle pair where the shoulder joint is controlled by anterior deltoid (shoulder flexor, M 1 ) and posterior deltoid (shoulder extensor, M 2 ) and elbow joint is controlled by brachialis (elbow flexor, M 3 ) and triceps brachii (elbow extensor, M 4 ) (Jagodnik and van den Bogert, 2010). The activations to these muscle groups (φ MN ) are transformed into joint angles for both shoulder and the elbow as follows, where, θ JA S and θ JA E are the joint angles of shoulder and elbow with respect to the x-axis (Supplementary Figure 1) and shoulderlength l S , respectively, in two-dimensional space, φ MN S Ag is the muscle activation of shoulder agonist muscle, φ MN S An is the muscle activation of shoulder antagonist muscle, φ MN E Ag is the muscle activation of elbow agonist muscle, and φ MN E An is the muscle activation of the elbow antagonist muscle.
The coverage of the arm in two-dimensional space is controlled by these joint angles. The joint angles are used to calculate the muscle lengths for both shoulder and elbow as given below. X targ , the target position; X arm , the current arm position; φ MN , the motor neuron activations; M L , muscle lengths; I gaba , inhibitory GABAergic current; DA e , extracellular dopamine; G (t), the MC output; G (t + 1), the BG-derived activity of thalamus.

Proprioceptive Cortex
The proprioceptive cortex (PC) is modelled as a self-organising map (SOM) (Kohonen, 2001) of size N PC x N PC where the sensory map of the arm was generated. Using muscle length vector (M L (t)) from the arm model (Equation 9) as a feature vector, PC SOM was trained. The activation of a single node (i, j) in the PC SOM is given as, where, W PC,i is the weight of the connection between the muscle length vector and i th the neuron of the two-dimensional PC network, M L is the muscle length vector and σ PC is the width of the Gaussian response of PC SOM.

The Prefrontal Cortex
The prefrontal cortex (PFC) encodes the goal position where, in real-time, the goal information is formed using the visual sensory feedback, which is passed on to the frontal areas. In our current model, we fix the goal or target position and denote it by X targ . The motor command initially is driven by the PFC as the PFC specifies the goal to be reached. Similar to the PC, the PFC SOM is trained using the target position vector as a feature vector. The input features of the PFC are the spatial locations where the arm can possibly reach in the two-dimensional space. The target locations that produce the activation in the PFC network is given as.
where, W PFC,ij are the weight of the connection between the target position vector and (i, j) th the neuron of the twodimensional PFC network, X targ is the target position and σ PFC is the width of the Gaussian response of PFC SOM.

Motor Cortex
Motor cortex (MC) is modelled as a combination of SOM and continuous attractor neural network (CANN) (Trappenberg, 2005) of size N MC x N MC . This type of architecture of MC is used to account for two distinct characteristics of cortical areas viz., low dimensional representation of input space and dynamics based on the connectivity in these cortical regions. A dynamic model like CANN is employed to facilitate the integration of multiple afferent inputs received from the PC, the BG, and PFC. The output activity of the MC CANN (G MC ) is defined by, where, g MC is the internal state of MC CANN, N MC is the size of MC network, b MC is the constant term. The internal state of the MC CANN g MC is given by, Frontiers in Computational Neuroscience | www.frontiersin.org where, C MC is the weight kernel representing lateral connectivity in MC CANN, which determines the local excitation/global inhibition dynamics, G MC is the output activity of MC CANN, I MC is the total input coming into MC CANN from PC, PFC, and BG and ⊗ represents the convolutional operation.

Lateral Connections in MC
The lateral connectivity in the MC CANN is characterised by short-range (local) excitation and long-range (global) inhibition whose dynamics are defined by the weight kernel W C MC is given by, where, i MC , j MC are the location of the nodes in MC, i h , j h corresponds to the central node, A C lat is the strength of the excitatory connections, K C is the global inhibition constant and σ C lat is the radius of the excitatory connections.

Total Inputs Into MC
The total input (I MC ) coming into MC CANN from PC (information about the current position of the arm), PFC (information about target position), and BG (error feedback signal) is given by, where, A PC , A PFC , A BG are the respective gains of PC, PFC, and BG, G PC , G PFC , G BG are the output activities of PC-derived SOM part of MC, PFC-derived activation part of MC, and BG-derived network activity of thalamus. PC activity is used to generate low-level feature maps in MC using the SOM algorithm. The activation of the (i, j) th node in the SOM part of the MC G PC,ij is given as, where, U C is the output activity of PC SOM network, W MC,i is the weight of the connection between the PC SOM network and i th neuron of the two-dimensional MC SOM network, and σ MC is the width of the Gaussian response of MC SOM. The input from PFC to MC (G PFC ) is the product of weight matrix (W PFC→MC ) and the output activity of PFC SOM is given by, where, U PFC is the output activity of the PFC SOM network, W PFC→MC is the weight matrix between PFC and MC.
In an earlier line of modelling studies, we have shown that the classical Go-NoGo interpretation of the functional anatomy of the BG must be expanded to Go-Explore-NoGo, in view of the putative role of the Indirect Pathway in exploration (Sridharan et al., 2006;Chakravarthy and Balasubramani, 2014). This series of models had resulted in the so-called Go-Explore-NoGo policy, which refers to a stochastic hill-climbing performed on the value function computed inside the BG (Magdoom et al., 2011). When the arm reaches the target (ǫ < 0.1), the connections between the PFC and MC are trained by, where, η PFC→MC is the learning rate between PFC and MC, G MC targ is the MC activation required for the arm to reach the target, and G MC PFC (G PFC ) is the MC activation due to PFC.

Motor Neurons
The output activity of MC CANN projects to the MN layer which consists of four MNs that drives four muscles of the arm whose activation is given by, where, A MN is the gain of MN, W MC→MN is the weight matrix between MC CANN and MN layer, and G MC is the output activity of MC CANN.
To close the sensory-motor loop, we perform a comparison with the initial activation to the MN layer that was used to produce desired activation φ MN D (t). The weights between the MN and MC are trained in a supervised manner by comparing the network-derived MN activation φ MN (t) to the desired activation φ MN D (t). This gives a loop that is consistent in mapping the external arm space to the neuronal space and vice. The connection between MC and MN is trained by, where, η MC→MN is the learning rate between MC and MN, φ MN D is the desired MN activation required for the arm to reach the target and φ MN is the network-derived MN activation due to MC, and G MC is the output activity of MC CANN. The training schema for the outer loop (sensory-motor loop) is described in Supplementary Information.

Inner Loop-Cortico-Basal Ganglia Loop
The functional pathway of the inner loop is shown in Figure 2B.
The inner loop consists of MC, BG, and thalamus. MC receives information from BG via the thalamus. MC sends information to BG based on the integration of incoming signals received from PFC (target goal position, X targ ), PC (current end-effector position of the arm, X arm ) and BG (via thalamus, error feedback signal, G BG ).

Basal Ganglia
Basal ganglia (BG) consists of the striatum, GPe, GPi, STN, and SNc. The output signal from BG provides the necessary control for the arm to reach the target by modulating the MC activity. The output of the MC is as given in Equation 12.

Value Computation and Stochastic Hill Climbing
The signal from the PC contains information about the current end-effector position of the arm (X arm ) whereas the signal from PFC contains the target goal position X targ . These two signals are combined in the BG to form a value function, V arm (t), that represents the error between the desired and the actual positions of the hand as, where, X targ is the target goal position, X arm is the current endeffector position of the arm, σ V is the spatial range over which the value function is sensitive for that particular target. The output of the BG performs a stochastic hill-climbing over the value function (Chakravarthy and Moustafa, 2018;Narayanamurthy et al., 2019) and drives the MC to facilitate the arm in reaching the target. The value difference (δ V ) which is computed by comparing the current and previous values is given as, where, V arm (t) is the current value and V arm (t − 1) is the previous value. Based on this value difference signal (δ V ), the striatum will send the inhibitory GABAergic current I gaba to the SNc neurons while the SNc neurons will in turn release dopamine into the extracellular space (DA e ), which is absorbed by the striatum. DA e is transformed into δ SNc V . δ SNc V modulates the selection of direct and indirect pathways in the BG. The dynamics between the striatum and the SNc are described in greater detail in the subsequent section, "The Central Loop."

Action Selection
Striatum: The resultant δ SNc V acts as a modulatory signal on D1R-MSNs and D2R-MSNs of the striatum, which processes the input signal, G MC (t), and send outputs y D1 & y D2 via direct and indirect pathways, respectively.
where, λ D1 and λ D2 represent the effect of dopamine (value difference) on the D1 and D2 MSNs, respectively, W CTX→D1 and W CTX→D2 represent connections between cortex and D1 MSNs and cortex and D2 MSNs, respectively, G MC is the output activity of MC, δ SNc V is the SNc-derived value difference, θ D1 and θ D2 are the thresholds of the D1 and D2 MSNs, respectively, a D1 and a D2 are the sigmoidal gains of the D1 and D2 MSNs, respectively. Since a D1 = −a D2 , the activation of direct and indirect pathways depends on the δ SNc V such that when δ SNc V is positive (negative) the direct (indirect) pathway is selected. Note that the λ D1 and λ D2 parameters in Equations (25-26) are dependent only on δ SNc V ("tonic dopamine") and not on its temporal derivative ("phasic dopamine").

STN-GPe subsystem:
In the indirect pathway, D2 MSNs of the striatum project to the GPe, where y D2 influences GPe neural dynamics, which in turn influences STN neural dynamics. STN-GPe forms a loop with inhibitory projections from GPe to STN and excitatory projections from STN to GPe. Such excitatoryinhibitory pairs of neuronal pools have been shown to exhibit limit cycle oscillations (Gillies et al., 2002) which was modelled as coupled Van der Pol oscillator (Kawahara, 1980). The dynamics of the STN-GPe system is defined as, where, x GPe and x TN are the internal states of GPe and STN neurons, respectively, y STN is the output of STN neuron, ε g and ε s are the strengths of lateral connections in GPe and STN networks, respectively, W glat and W slat are weight kernels representing lateral connectivity in GPe and STN networks, respectively, y D2 is the output of D2 MSN, τ GPe and τ STN are the time constants of GPe and STN, respectively, w sg is the connection strength from STN to GPe, w gs is the connection strength from GPe to STN, and λ STN is the parameter which controls the slope of the sigmoid in STN.
Lateral Connections in STN-GPe: The lateral connectivity in STN or GPe network is modelled as Gaussian neighbourhood (Muddapu et al., 2019) which is defined by the weight kernel W glat/slat as, where, d 2 i,j,k,l is the distance of neuron i, j from a centre neuron k, l , σ g/s lat is the spread of the lateral connections for GPe or STN network. The detailed analysis of the STN-GPe subsystem is described in section STN-GPe Dynamics of the Supplementary Material.
GPi: The signals arriving from D1 MSN y D1 and STN y STN via direct and indirect pathways, respectively, combines in GPi which is defined as, where, y D1 is the output of D1 MSN via direct pathway, y STN is the output of STN via indirect pathway, A D1 and A D2 are the gains of direct and indirect pathways, respectively.

Thalamus
The combined inputs y GPi at GPi from direct y D1 and indirect y STN pathways are then passed on to the thalamus. To integrate and philtre the information from the GPi output, the thalamus was modelled as a CANN which is defined as, where, g thal is the internal state of thalamus CANN, N thal is the size of thalamus network, b thal is the constant term. The internal state of the thalamus CANN g thal is given by, where, W C thal is the weight kernel representing lateral connectivity in thalamus CANN, which determines the local excitation/global inhibition dynamics, G thal is the output activity of thalamus CANN, I BG is the total input coming into thalamus CANN from BG, y GPi is the output of GPi, G BG is the BG-derived network activity of the thalamus, and ⊗ represents the convolution operation.

Central Loop-Nigro-Striatal Loop
The functional pathway of the central loop is as represented in Figure 2C. The central loop consists of the striatum (the input nucleus of BG) and SNc. SNc projects to the striatum via dopaminergic axons (DA e ) and striatum in turn projects to SNc via inhibitory GABAergic axons I gaba . Based on the sensory feedback signal received from the PC (X arm ) and the target information from the PFC X targ , the striatum performs value computation (V arm ). Based on the values from current (V arm (t)) and previous (V arm (t − 1)) instants, the value difference (error, δ V ) is computed. Based on the value difference (δ V ), appropriate amount of GABAergic current I gaba is sent to SNc, which in turn releases dopamine into the striatum (DA e ) accordingly.

SNc Neuron (Soma)
The detailed single-compartmental biophysical model of the SNc neuron is adopted from Muddapu and Chakravarthy (2021). The model incorporates all the essential molecular level mechanisms such as ion channels, active pumps, ion exchangers, dopamine turnover processes, etc.
Based on the value difference signal (δ V ), the inhibitory GABAergic current I gaba , flows from the striatum to SNc. I gaba along with excitatory glutamatergic current I nmda/ampa contributes to the overall synaptic input current flux J syn to the SNc neurons. J syn regulates the membrane voltage of the SNc along with the sodium, calcium, and potassium fluxes as given by, where, F is the Faraday's constant, C snc is the SNc membrane capacitance, vol cyt is the cytosolic volume, AR pmu is the cytosolic area, J m,Na is the sodium membrane ion flux, J m,Ca is the calcium membrane ion flux, J m,K is the potassium membrane ion flux, J syn is the overall input current flux, δ V is the value difference, I gaba is the inhibitory GABAergic current flux, and I nmda/ampa is the excitatory glutamatergic (NMDA/AMPA) current flux. The membrane voltage of SNc (V SNc ) regulates the calcium membrane ionic flux which results in calcium oscillations inside SNc neuron. The calcium membrane ionic flux (J m,Ca ) is given by, (40) where, F is the Faraday's constant, z Ca is the valence of calcium ion, vol cyt is the cytosolic volume, I CaL is the L-type calcium channel current, I pmca is the ATP-dependent calcium pump current, and I NaCaX is the sodium-potassium exchanger current. The intracellular calcium oscillation or dynamics ([Ca i ]) is defined as, where, J m,Ca is the flux of calcium ion channels, J calb is the calcium buffering flux by calbindin, and J cam is the calcium buffering flux by calmodulin. A detailed description of the SNc neuron is provided in section Biophysical Model of SNc of the Supplementary Material.

SNc Terminal
The three-compartmental biochemical model of the SNc terminal is adopted from Muddapu and Chakravarthy (2021). The SNc terminal model incorporates all the necessary molecular-level mechanisms of the dopamine turnover process such as synthesis, packing, release, and reuptake.
Calcium-Dependent Dopamine Release: Dopamine (DA) synthesis and release by SNc terminal depend on calcium oscillations. The flux of dopamine release (J rel ) from the SNc terminal is given by, where, [Ca i ] is the intracellular calcium concentration in the SNc terminal, P rel is the release probability as a function of intracellular calcium concentration, n RRP is the average number of readily releasable vesicles, and ψ is the average release flux per vesicle within a single synapse.
Calcium-Dependent Dopamine Synthesis: The flux of calciumdependent dopamine synthesis is defined as, where, K synt is the calcium sensitivity, V synt is the maximal velocity for L-DOPA synthesis, and [Ca i ] is the intracellular calcium concentration. The flux of synthesised L-DOPA, J synt , whose velocity is the function of intracellular calcium concentration and L-DOPA synthesis is regulated by the substrate (TYR) itself, extracellular DA (via autoreceptors) and intracellular DA concentrations, is given by, where, V synt is the velocity of synthesising L-DOPA, [TYR] is the tyrosine concentration in terminal bouton, K TYR is the tyrosine concentration at which half-maximal velocity was attained, K i,cda is the inhibition constant on K TYR due to cytosolic DA concentration, where, J rel represents the flux of calcium-dependent DA release from the DA terminal, J DAT represents the unidirectional flux of DA translocated from the ECS into the intracellular compartment (cytosol) via DA plasma membrane transporter (DAT), J o eda represents the outward flux of DA degradation, which clears DA from ECS, and δ SNc V is the SNc-derived value difference. A detailed description of the SNc terminal is provided in section Biophysical Model of SNc of the Supplementary Material.
The cortical input to the striatum is modulated by the δ SNc V derived from the network of SNc neurons. When δ SNc V is high, the direct pathway will be selected, else the indirect pathway is selected.

Time Scales of Various Loops
The time scales of various loops in the model are as given in Table 1. The STN-GPe loop runs with a timestep (dt) of 0.02 ms. Once the STN-GPe loop runs for 2,500 iterations, one timestep of the cortico-BG loop is run. Simultaneously, the SNc-STR loop which provides the modulatory signal for the selection of the Go-NoGo pathway in the striatum is run. For each timestep of the cortico-BG loop, STN-GPe and SNc-STR loops run 2,500 and 2,000 iterations, respectively. The total simulation time for the arm reaching task (cortico-BG loop) is 5 s and if the arm doesn't reach the target in the stipulated timeframe of 5 s, the trial is considered non-reachable. At the spatial level, the details of different loops are given in Figure 2.

Simulating Parkinsonian Conditions
To simulate the Parkinsonian condition in the present model, the number of neurons in the SNc population (network) was reduced. In order to kill the SNc neuron, we clamped their membrane voltage (V SNc ) to resting membrane voltage (−80 mV). As the number of SNc neurons die the total amount of dopamine (DA e ) that is made available to the striatum decreases. This influences the selection of the indirect pathway in the BG system over the direct pathway resulting in pathological conditions. In the present model, two types of PD conditions were simulated: in the first type, SNc cell loss affects striatum alone (PD1) and in the second type, SNc cell loss affects both striatum and STN (PD2).
In normal conditions, the SNc-derived value difference δ SNc V will be similar to the actual value difference computed (δ V ).
In case of PD1, the SNc-derived value difference δ SNc V will be lesser than the actual value difference computed (δ V ). In the case of PD2, along with δ SNc V < δ V , δ SNc V impacts the STN lateral connections, thereby influencing the complexity of the STN-GPe subsystem. The STN-GPe subsystem is an integral component of the indirect pathway and is believed to play a major role in exploratory behaviour (Sridharan et al., 2006;Chakravarthy and Balasubramani, 2014).
In normal condition: δ SNc V = F (DA e ) DA e = SNc I gaba , P SNc ; P SNc = 100% In PD1 condition: δ SNc V = F (DA e ) DA e = SNc I gaba , P SNc ; P SNc < 100% In PD2 condition: where, δ SNc V is the SNc-derived value difference, δ V is the value difference computed, DA e is the extracellular dopamine, I gaba is the inhibitory GABAergic current from the striatum, P SNc is the percentage of SNc neurons, and ε s is the lateral connection strength in the STN network.

Levodopa Medication
When a drug is administered to a patient, the medication action is broadly classified into two major branches: pharmacokinetics (what the body does to the drug) and pharmacodynamics (what the drug does to the body) (Shanbhag and Shenoy, 2020).

Pharmacokinetics
Pharmacokinetics deals with the absorption, distribution, metabolism, and excretion of drugs. In the present study, we have adapted a two-compartment pharmacokinetic model of levodopa (L-DOPA) (Baston et al., 2016), which consists of central and peripheral compartments (Figure 3). Orally consumed L-DOPA is absorbed in the intestine and reaches the bloodstream. The bloodstream carries the drug all over the body. Proteins break down L-DOPA and around three-fourth of the drug is deactivated before it even reaches the brain. The central compartment where L-DOPA is administered and plasma L-DOPA concentration was measured which is defined as, compartment, k 21 is the rate constant from peripheral to central compartments, k 12 is the rate constant from central to peripheral compartments, and k 1e is the total clearance rate constant from the central compartment.
The interaction between plasma L-DOPA and other body fluids, which occurs in the peripheral compartment, is defined as, (51) where, V PC is the volume of the peripheral compartment, [LDOPA CC ] is the L-DOPA concentration in the central compartment, [LDOPA PC ] is the L-DOPA concentration in peripheral compartment, k 21 is the rate constant from peripheral to central compartments, and k 12 is the rate constant from central to peripheral compartments.

Pharmacodynamics
Pharmacodynamics deals with molecular, biochemical, and physiological effects of drugs, including drug mechanism of action, receptor binding (including receptor sensitivity), postsynaptic receptor effects, and chemical interactions. In the present study, we have adapted a three-compartment dopaminergic terminal model (Reed et al., 2012) which consists of extracellular, vesicular, and cytoplasmic compartments. When L-DOPA medication is administered, the flux of exogenous L-DOPA ([LDOPA CC ]) transported into the terminal through aromatic L-amino acid transporter (AAT) while competing with other aromatic amino acids [such as tyrosine (TYR) and tryptophan (TRP)] (Reed et al., 2012) is given by, where, K ldopa e is the extracellular L-DOPA concentration at which half-maximal velocity was attained, V aat is the maximal velocity

Timescales in the Model
Reaching movements, like several other behavioural events, involve dynamics at multiple timescales: the neuronal activity which is generally in milliseconds, and the actual movement which unfolds over the order of seconds. In the present model, the outer (sensory-motor) loop is assumed to run slightly slower than the inner (cortico-basal ganglia) and central (nigrostriatal) loops. As the dynamics of the STN-GPe loop in the indirect pathway needs some time to settle, we run this loop for 2, 500 iterations (dt = 0.02 ms), before sending the output to the MC (MC runs for 100 iterations with dt = 50 ms). Thus, a single update of the MC activity happens after every 50 ms during which the BG dynamics run. Similarly, since the dynamics of the SNc neuron needs some time to settle, we run the SNc neuron for 2, 000 iterations (dt = 0.025 ms), before sending the output to the BG. Thus, a single update of the MC activity happens after every 50 ms during which the SNc dynamics run. All the results presented are at the timescale of the MC. In the present model, the SNc neurons run in milliseconds timescale whereas the pharmacokinetic-pharmacodynamic model of L-DOPA medication runs in hourly timescale. In order to show the drug effect, we sample various points across the L-DOPA medication curve (Supplementary Figure 7.1) and simulated the MCBG model for the arm reaching task for each sampled point.

RESULTS
Here, we showcase the performance of the model by simulating the PD condition and read out their effects on behavioural outcomes (Figures 4, 5). Furthermore, demonstrates the effect of differential dopaminergic axonal loss manifest into some of the cardinal symptoms of PD (Figures 6, 7). Next, assessing the performance in terms of reaching time and verifying the effect of L-DOPA therapeutic intervention (Figures 8, 9). Finally, describing the model results which gave an indicator of how to optimise the drug dosage across the course of the disease progression (Figures 10-12).

Testing Phase
The performance of the MCBG model was tested, the model results were compared to that of the CBG model (Muralidharan et al., 2018) and the experimental data (Majsak et al., 1998) for both control and PD conditions. In the MCBG model, PD conditions simulated were subdivided into two categories: in PD1, the SNc cell loss impacts only striatum whereas in PD2, the SNc cell loss impacts both striatum and STN. The MCBG and CBG models were tested and the performance was evaluated with respect to movement time, peak velocity, time-to-peak velocity, and average velocity along with the experimental results. In the control case, the MCBG model reaches the target in 0.46 ± 0.02 s compared to the CBG model and the experimental subject which reaches the target in 0.56 ± 0.1 s and 0.3432 ± 0.04 s, respectively ( Figure 4A, dark blue bar). The MCBG model obtained a peak velocity of 2.23 ± 0.05 ms −1 compared to the CBG model and experimental subject which obtained peak velocity of 1.88 ± 0.15 ms −1 and 2.15 ± 0.27 ms −1 , respectively, during the arm reaching toward the target in case of control ( Figure 4C, dark blue bar). The time taken to reach the peak velocity in the case of control was 0.21 ± 0.02 s for the MCBG model, 0.29 ± 0.09 s for the CBG model, and 0.19 ±0.02 s for the experimental subject ( Figure 4B, dark blue bar). Finally, the average velocities for MCBG and CBG models were found to be 1.49 ± 0.05 ms −1 and 1.26 ± 0.15 ms −1 , respectively, in the case of control (Figure 4D, dark blue bar).
In the case of PD, the experimental subject recorded an average movement time of 0.52 ± 0.63 s, respectively (Figure 4A, cyan bar), while the CBG model reaches the target in 1.17 ± 0.63 s (Figure 4A, cyan bar) whereas the MCBG model took 1.88 ± 1.42 s and 1.6 ± 1.35 s for PD1 ( Figure 4A, cyan bar) and PD2 (Figure 4A, yellow bar), respectively. The experimental subject recorded a peak velocity of 1.35 ± 0.18 ms −1 (Figure 4C, cyan bar) compared to the CBG model which obtained a peak velocity of 1.74 ± 0.13 ms −1 (Figure 4C, cyan bar) whereas the MCBG model obtained peak velocities of 1.18 ± 0.35 ms −1 (Figure 4C, cyan bar) and 0.98 ± 0.31 ms −1 (Figure 4C, yellow bar), respectively, during the arm trajectory toward the target. The time taken to reach the peak velocity in the PD case was 0.27 ± 0.03 s for the experimental subject ( Figure 4B, cyan bar), 0.35 ± 0.07 s for CBG model (Figure 4B, cyan bar) and 0.56 ± 0.28 s, and 0.79 ± 0.35 s in PD1 (Figure 4B, cyan bar) and PD2 (Figure 4B, yellow bar) cases, respectively, for MCBG model. Finally, the average velocity for the CBG model was found to be 0.77 ± 0.21 ms −1 in PD (Figure 4D, cyan bar), and the average velocities for the MCBG model were found to be 0.68 ± 0.27 ms −1 and 0.59 ± 0.23ms −1 in PD1 (Figure 4D, cyan bar) and PD2 (Figure 4D, yellow bar), respectively.

Simulating Parkinsonian Conditions
To simulate PD conditions in the model, SNc cells were killed and their effects on basal ganglia were considered in two aspects. In the first scenario, only the striatum is affected by SNc cell loss (PD1-cell loss affecting nigrostriatal pathway only) and in the second scenario, both striatum and STN are affected by SNc cell loss (PD2-cell loss affecting both nigrostriatal and nigrosubthalamic pathways).
FIGURE 4 | Comparison of performance of the proposed model (during the testing phase) with CBG model (Muralidharan et al., 2018) and experimental data adapted from (Majsak et al., 1998).

Effect of SNc Cell Loss on MCBG Behavioural Outcome
To assess the performance metrics with respect to dopaminergic cell loss affecting striatum and both striatum and STN, a comparison study was done with respect to the movement time, peak velocity, time required to peak velocity, and average velocity (Figure 5). In both cases (PD1 and PD2), the time required to reach the target ( Figure 5A) and time-to-reach the peak velocity ( Figure 5B) increases with an increase in SNc cell loss. In the PD1 case, the peak velocity increases with an increase in SNc cell loss when compared to the PD2 case where the peak velocity decreases with an increase in SNc cell loss ( Figure 5C). The reason behind this discrepancy in both cases will be explored in the next sections where one leads to tremor-like behaviour and the other leads to rigidity-like behaviour. In both cases, the average velocity across the trajectory decreases with an increase in SNc cell loss ( Figure 5D).

Differential Dopaminergic Axonal Degeneration Manifests Into Various PD Motor Symptoms
Both the PD scenarios (PD1 and PD2) simulated in the model can be attributed to differential degeneration of dopaminergic projections to various targets in the basal ganglia, and how degeneration manifests into various motor symptoms of PD. In the control case, the arm reaches the target in 0.55 s (Figure 6Ai) with the peak velocity of 1.91 ms −1 (Figure 6Aii). The population activity of STN exhibits desynchronous activity during the arm movement which is indicated in the STN spectrogram (Figure 6Aiii) and synchrony (average value = 0.03) (Figure 6Aiv) (synchrony measure is described in section Network analysis of the Supplementary Material). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 264 nM which was in the range of 150-400 nM (Schultz, 1998 ; Figure 6Av).
In 25% PD1, the arm reaches the target in 1.5 s (Figure 6Bi) with a reduced peak velocity of 0.71 ms −1 , exhibiting bradykinesia-like behaviour in the arm (Figure 6Bii). Population activity of STN exhibits greater synchrony compared to control case during the arm movement which is also indicated in STN spectrogram (Figure 6Biii) and synchrony with an average value of 0.11 (Figure 6Biv). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 148 nM which was lesser than in the control case (Figure 6Bv).
In 25% PD2, the arm reaches the target in 1.2 s (Figure 6Ci) with the peak velocity of 0.71 ms −1 , exhibiting bradykinesialike behaviour in the arm (Figure 6Cii). Population activity    of STN exhibits desynchronous activity, same as control case during the arm movement which is indicated in the STN spectrogram (Figure 6Ciii) and synchrony with an average value of > 0.01 (Figure 6Civ). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 154 nM which was lesser than the control case (Figure 6Cv).
In 50% PD1, the arm reaches the target in 2.7 s (Figure 6Di) with the peak velocity of 0.84 ms −1 showing tremor-like behaviour in the arm (Figure 6Dii). Population activity of STN exhibits low synchronous activity during the arm movement which indicates in STN spectrogram (Figure 6Diii) and synchrony with an average value of 0.17 (Figure 6Div). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 101 nM which was lesser than the control case (Figure 6Dv). In 50% PD2, the arm reaches the target in 4.7 s (Figure 6Ei) with the peak velocity of 0.84 ms −1 as a result of cogwheel-like behaviour in the arm (Figure 6Eii). The population activity of STN exhibits high synchronous activity during the arm movement which indicates in the STN spectrogram (Figure 6Eiii) and synchrony with an average value of 0.55 (Figure 6Eiv). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 90 nM which was lesser than the control case (Figure 6Ev).
In 75% PD1, the arm did not reach the target within 5 s (Figure 6Fi) with the peak velocity of 1.54 ms −1 displaying a tremor-like behaviour in the arm (Figure 6Fii). Population activity of STN exhibits low synchronous activity during the arm movement which indicates in STN spectrogram with increased power in 5 − 25 Hz region (Figure 6Fiii) and synchrony with an average value of 0.15 (Figure 6Fiv). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 51 nM which was lesser than the control case (Figure 6Fv). In 75% PD2, the arm did not reach the target within 5 s (Figure 6Gi) with zero peak velocity as a result of the rigidity-like state of the arm (Figure 6Gii). Population activity of STN exhibits high synchronous activity during the arm movement which is indicated in STN spectrogram with increased power in 15−50 Hz region (Figure 6Giii) and synchrony with an average value of > 0.99 (Figure 6Giv). Dopamine released by SNc neurons in the striatum during the arm reaching peaked at ∼ 13 nM, which was lesser than in the control case (Figure 6Gv).

Quantifying Tremor-Like and Rigidity-Like Motor Symptoms
To quantify between tremor-like and rigidity-like motor symptoms of PD, root mean square (RMS) acceleration was computed across movement trajectory for various PD conditions where RMS acceleration can be used as an indicator of random non-deterministic movements (Figure 7). In the PD1 scenario, the RMS acceleration increases with an increase in SNc cell loss which indicates irregular changes in the velocity of arm movement ( Figure 7A). This irregular velocity profile in PD1 is a result of tremor-like motor behaviour. In the PD2 scenario, the RMS acceleration increases with an increase in SNc cell loss to 50%, and beyond 50% RMS acceleration decreases with an increase in SNc cell loss ( Figure 7B). The tremor-like motor behaviour is indicated by the RMS acceleration increases until 50% SNc cell loss and from there on, we can see a sudden decrease, which marks the onset of rigidity.

Effect of Levodopa Medication
In order to show the L-DOPA medication effect on the MCBG model, we simulated different scenarios where various L-DOPA  dosages were administrated across various PD conditions and movement time was monitored.

Comparison of MCBG Model With Experimental Results
The L-DOPA therapeutic effect was monitored by recording the performance in terms of the average movement time across the time course of the dosage for the next 10 h. The performance of the model was also recorded 2 h prior to the administration of the drug. The MCBG model results were compared with experimental studies where PD patients were evaluated based on UPDRS Part III score (Nomoto et al., 2018 ; Figure 8). The experimental PD subjects were categorised into two groups based on the UPDRS part III score (motor evaluation) where the group 1 PD subjects have a mean UPDRS III score of 28 (13-51) and the group 2 PD have a mean UPDRS III score of 30.3 (22-41) (Nomoto et al., 2018). An average L-DOPA dosage of 141 mg was given to both the experimental groups. The MCBG model was simulated with 62% SNc cell loss and 150 mg of L-DOPA administered at the second hour of the simulation. The PD1 MCBG model performance in terms of movement time ( Figure 8A, blue curve) matched with experimental group 2 result in terms of UPDRS III score ( Figure 8A, orange curve). Similarly, PD2 MCBG model performance in terms of movement time ( Figure 8B, blue curve) matched with experimental group 1 result in terms of UPDRS III score ( Figure 8B, orange curve).

Effect of L-DOPA Medication With Disease Progression
The effect of L-DOPA (150 mg) medication on the model performance was studied across different percentages (25, 37, 50, 62, and 75%) of SNc cell loss for both PD1 and PD2 scenarios. The L-DOPA medication was given at the second hour in the simulation. The simulated results show that as SNc cell loss increases, the model performance deteriorates, and also the therapeutic effect decreases as the disease progresses in both PD1 and PD2 scenarios (Figure 9). The maximum therapeutic effect of L-DOPA was seen for 50% and 62% SNc cell loss in both PD1 and PD2 scenarios (Figures 9E-H). In 75% SNc cell loss, the model performance was poor in case of PD1 when compared to PD2 (Figures 9I,J). The model performance was categorised into three regions based on the following criteria: If the arm reaches the target within 2 s, then that region was marked in green colour which indicates the normal movement. If the arm reaches the target between 2 and 4 s then that region was marked in yellow colour, indicating slow movement or bradykinesia. If the arm reaches the target beyond 4 s, then that region was marked in red colour which indicates very slow movement or akinesia. The simulated results show that as the SNc cell loss increases the movement time curve shift from green to the yellow region when medication was ON and the movement time curve shift from yellow to the red region when medication was OFF (Figure 9).

Effect of L-DOPA Dosage and SNc Cell Loss on Therapeutic Window
As discussed in the previous section, the model performance was categorised into three regions: green (normal movement), yellow (slow movement, bradykinesia), and red (very slow movement, akinesia). The therapeutic window is computed by taking the time difference between the points when the performance improved after taking medication and entered into the green shaded region until it started wearing off and crosses back to the yellow shaded region (where the effects of L-DOPA start wearing off).
In the case of 25%, SNc cell loss (PD1), as the L-DOPA dosage increases the therapeutic window (green region) decreases (Figure 10, first column). But at higher percentage loss of cells (37, 50, 62, and 75% SNc cell loss), as the L-DOPA dosage increases the therapeutic window (green region) increased (Figure 10). However, in the case of PD2 for all percentages of SNc cell loss, as the L-DOPA dosage increases the therapeutic window (green region) increased (Figure 11).

MCBG Model
The proposed model tries to present a biologically realistic model of the effect of L-DOPA on PD symptoms, specifically in terms of movement parameters. In our modelling approach, a large-scale cortico-basal ganglia model forms the backbone of our network. The two-link arm model that is interfaced to the MNs simulates the movement of the hand and the feedback related to the hand position and distance from the target is processed by the PC and passed on to MC. MC uses the corrective signals from the BG to initiate the next action. The BG dynamics are highly influenced by the dopaminergic input from the SNc and by incorporating a detailed biophysical model of the SNc into the network model, we were able to show the effect of loss of dopaminergic cells on the movement parameters. Going forward we aim to relate the pathological behaviour with respect to the dynamics at the molecular level happening inside the SNc.

Differential Projections and PD Symptoms
The proposed model was able to explain a wide range of pathological behaviours associated with PD by controlling the release of dopamine into the extracellular space and reducing the complexity of the STN-GPe network. We modelled differential projections from the SNc to the Striatum as well as from SNc to STN. By reducing the supply of dopamine through the SNc to Striatum projections, the slowness of movement or bradykinesia could be simulated, and in combination with modulating the complexity of the STN-GPe network through the SNc to STN projections, symptoms like tremor and rigidity were simulated. The complexity of the STN-GPe network was varied by controlling the dopaminergic projections of the SNc neurons toward the STN, thereby affecting the lateral connections within the STN subsystem. By progressively reducing the number of dopaminergic cells in SNc, we could replicate some of the cardinal symptoms of PD-bradykinesia, tremor, and rigidity.

L-DOPA Medication Effect
Once the PD condition and the associated symptoms were simulated, we integrated a pharmacokinetic-pharmacodynamic (PK-PD) model of L-DOPA medication (Baston et al., 2016;Véronneau-Veilleux et al., 2020), which showed improved results in reaching performance. L-DOPA medication is one of the first-line treatment methodologies for Parkinson's disease (Suzuki et al., 2020). Our model incorporates the medication effect by interfacing the SNc with the PK-PD model of L-DOPA drug administration. Depending on the dosage of the drug administered, L-DOPA is absorbed into the blood. After interacting with other bodily fluids, a portion of the L-DOPA crosses the BBB and gets absorbed by the dopaminergic terminals. Our results show that consumption of L-DOPA improves the PD symptoms to a great extent. Using our model, we could also see that the extent of improvement on the PD condition depend on the dosage (Figure 12). A higher level of serum L-DOPA results in dyskinesias and a low-level result in wearing off. Hence, an optimum dosage of medication has to be selected. In order to optimise the drug dosage, we performed our tests with various dosages of L-DOPA medication. We could see that as the percentage of SNc cell loss increases, a higher dosage of L-DOPA was required to sustain the medication effect. With the increase in the percentage of SNc cell loss, the therapeutic effect keeps decreasing. Hence our study focused on the variation of therapeutic effect with respect to the varying percentage SNc cell loss and L-DOPA dosage. The results observed are promising enough to suggest optimal tuning strategies of drug dosage for PD patients (Figure 12). The performance characteristics with respect to the variation in cell loss and the dosage help us to tune the optimum dosage in terms of the quantity and the frequency of dosage.

Side-Effects of L-DOPA Medication
From the simulation results, we can explain the L-DOPA wearing-off mechanism to a great extent. Our hypothesis is that the natural progression of the disease characterised by the increase in loss of SNc cells is one of the mechanisms that contribute to L-DOPA wearing off. There could be other factors as well that can accelerate this wearing-off phenomenon. Another hypothesis is that the loss of dopaminergic terminals will lead to synchronised activity in STN which in turn causes overexcitation of SNc neurons resulting in a phenomenon called excitotoxicity in SNc (Muddapu et al., 2019;Muddapu and Chakravarthy, 2020). Thus, fewer dopaminergic terminals and higher L-DOPA dosage result in an accelerated loss of the dopaminergic terminals leading to a faster wearing-off. There might be other contributing factors as well that may advance the shortening of the therapeutic window. There is the potential scope of carrying out a detailed study on the various causes of the L-DOPA wearing off and we believe our model serves as a good platform to conduct such comprehensive research. As shown in our results indicated in Figure 12, we also observe a decrease in performance and reduction in the size of the therapeutic window with an increase in LDOPA dosage beyond a certain value. For example, we can observe the plots of Figure 12A, for 37% cell loss, and Figure 12B, for both 37 and 50% cell loss that if LDOPA dosage is increased beyond 250 mg, the therapeutic window reduces. This reduction in performance can be attributed to Dyskinesias. Hence this model simulation also helps us to optimise the drug dosage with respect to the severity of the disease and the dosage of medication.
The comparison of previous computational models of BG were shown in the Table 2. We can see that the proposed model covers almost all aspects including simulating the PD condition and Behavioural outcomes. It is able to simulate all the cardinal symptoms of TD except the postural imbalance. It is more biologically plausible as a detailed model of SNc for DA release is used and it can simulate the medication effect. Compared to the other models mentioned the current model proposed is having better coverage.

Future Scope
We could reliably replicate some of the cardinal symptoms of PD using our MCBG model. Along with simulating the PD ON/OFF mechanisms, our model could also successfully demonstrate the medication effect of L-DOPA. With the L-DOPA PK-PD model integration with the MCBG model, we could also explain the side effects of L-DOPA medication such as dyskinesias and wearing off. Hence this model has the scope to be developed into a test bench for PD. The current diagnostics and treatment methodologies for PD are based on direct observation and therefore suffer from subjectivity (Nair et al., 2022) and we believe that our model can be developed further to provide a more quantitative approach to diagnose PD symptoms and optimise therapeutic interventions.

Understanding Causes of Wearing off Mechanism and Dyskinesias
We hypothesise that the natural progression of the disease and the excitotoxicity could be potential factors that result in L-DOPA wearing off. An increase in cytosolic DA will lead to excitotoxicity as unregulated cytosolic DA leads to neurodegeneration (Chen et al., 2008). In this line, the pharmacological model can be extended by incorporating the administration of other drugs that block the vesicular transporter (Pregeljc et al., 2020). In addition to dopamine-induced excitotoxicity, L-DOPA-induced toxicity can also cause neurodegeneration (Fahn, 2005;Lipski et al., 2011;Witt and Fahn, 2016;Muddapu et al., 2020b). However, there could be other contributing factors too and this model can serve as a starting step to explore research in a similar direction. As highlighted in the discussion section, a more detailed study of the L-DOPA wearing-off mechanism can be carried out to understand the mechanism and devise alternate or improved medication regimes. Another line of extension is to explore the phenomenon of different types of dyskinesias such as peak dosage and diphasic dyskinesias (Kim Y. E. et al., 2019).

Incorporating DBS to Address Dyskinesias
We also want to extend the model to show the effect of deep brain stimulation (DBS) on motor deficiencies in PD condition and explore the comorbidity effects of both L-DOPA and DBS on PD motor symptoms (Muthuraman et al., 2018;Muddapu et al., 2019;Muddapu and Chakravarthy, 2020;Mueller et al., 2020). One of the limitations of our model is that our model does not consider the influence of the hyperdirect pathway, which involves direct cortical connections to the STN (Nambu et al., 2002;Cai et al., 2019). Also, the model does not take into consideration the influence of cholinergic interneurons in the striatum (Crossley et al., 2016;Kim T. et al., 2019). These can be considered as further enhancements to the current model. Currently, our model is focusing on the motor deficiencies in the PD pathology. It would be interesting to model PD non-motor symptoms (Goldman and Postuma, 2014;Goldman and Guerra, 2020).

CONCLUSION
A comprehensive test bench for demonstrating the effect of drug action on symptoms can be a powerful tool in the therapeutic toolkit of neurodegenerative diseases such as Parkinson's disease. Our model is the first step toward this bigger goal. In the current study, we were able to successfully simulate the relationship between drug dosage, cell loss, and PD ON and OFF conditions. We could also demonstrate some of the cardinal symptoms of PD. We also integrated a PK-PD model of L-DOPA medication, which enabled us to simulate the medication effects of the L-DOPA. We also simulated various combinations of L-DOPA medication and percentage of SNc cell loss which enabled us to understand the general trends in drug effects.
These modelling results have the potential to optimise the medication in terms of the amount of dosage and the frequency of dosage.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Further, inquires can be directed to the corresponding author. The MATLAB code of the proposed MCBG model (http:// modeldb.yale.edu/266907) is available on the ModelDB server (McDougal et al., 2017) and an access code will be provided on request.