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₁) and posterior deltoid (shoulder extensor, M₂) and elbow joint is controlled by brachialis (elbow flexor, M₃) and triceps brachii (elbow extensor, M₄) (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, θSJA and θEJA are the joint angles of shoulder and elbow with respect to the x-axis (Supplementary Figure 1) and shoulder-length (l_S), respectively, in two-dimensional space, ϕSAgMN is the muscle activation of shoulder agonist muscle, ϕSAnMN is the muscle activation of shoulder antagonist muscle, ϕEAgMN is the muscle activation of elbow agonist muscle, and ϕEAnMN 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.

where, μAgS, μAnS, μAgE, and μAnE are the agonist and antagonist muscle lengths of shoulder and elbow, respectively, a_S is the distance between the shoulder join centre and M₁ or M₂ moment lever, b_S is the distance between the shoulder joint centre and M₁ or M₂ moment lever, a_E is the distance between elbow joint centre and M₃ or M₄ moment lever, and b_E is the distance between elbow joint centre and M₃ or M₄ moment lever.

Using these muscle lengths in the form of a four-dimensional vector (ML=[μAg SμAn SμAgE μAnE]), a sensory (proprioceptive) map of the arm was generated. The end effector position of the arm (Xarm=[x1arm x2arm]) in the two-dimensional space is calculated as,

where, θSJA and θEJA are the joint angles of shoulder and elbow with respect to the x-axis (Supplementary Figure 1) and shoulder-length (l_S), respectively, in two-dimensional space, l_S is the distance between the shoulder joint centre (S) and elbow joint centre (E), l_E is the distance between the elbow joint centre (E) and end effector (H), a_S is the distance between the shoulder joint centre and M₁ or M₂ moment lever, μAgS, μAnS, μAgE, and μAnE are the agonist (M₁ or M₃) and antagonist (M₂ or M₄) muscle lengths of shoulder and elbow, respectively.

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 (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 two-dimensional PFC network, X^targ is the target position and σ_PFC is the width of the Gaussian response of PFC SOM.

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,

where,  MCC 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.

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 ϕDMN(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 ϕDMN(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, ϕDMN 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.

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.

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, WthalC 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.

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 (δVSNc) will be similar to the actual value difference computed (δ_V). In case of PD1, the SNc-derived value difference (δVSNc) will be lesser than the actual value difference computed (δ_V). In the case of PD2, along with δVSNc<δV, δVSNc 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:

In PD1 condition:

In PD2 condition:

where, δVSNc 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.

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

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.

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

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⁻¹ (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⁻¹, 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⁻¹, exhibiting bradykinesia-like 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⁻¹ 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⁻¹ 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⁻¹ 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).

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.

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

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

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

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.

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.

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.

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

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