Increasing Serotonin to Reduce Parkinsonian Tremor

While current dopamine-based drugs seem to be effective for most Parkinson's disease (PD) motor dysfunctions, they produce variable responsiveness for resting tremor. This lack of consistency could be explained by considering recent evidence suggesting that PD resting tremor can be divided into different partially overlapping phenotypes based on the dopamine response. These phenotypes may be associated with different pathophysiological mechanisms produced by a cortical-subcortical network involving even non-dopaminergic areas traditionally not directly related to PD. In this study, we propose a bio-constrained computational model to study the neural mechanisms underlying a possible type of PD tremor: the one mainly involving the serotoninergic system. The simulations run with the model demonstrate that a physiological serotonin increase can partially recover dopamine levels at the early stages of the disease before the manifestation of overt tremor. This result suggests that monitoring serotonin concentration changes could be critical for early diagnosis. The simulations also show the effectiveness of a new pharmacological treatment for tremor that acts on serotonin to recover dopamine levels. This latter result has been validated by reproducing existing data collected with human patients.


INTRODUCTION
Resting tremor is one of the most disabling features of Parkinson's disease (PD). People affected by resting tremor exhibit uncontrollable movements involving a body part, for example, the arm, when at rest. Tremor tends to decrease or stop when they deliberately move the affected part of the body (Deuschl et al., 2000;Kalia and Lang, 2015). Like other PD motor features, tremor is thought to result primarily from the death of dopamine-producing cells in the substantia nigra pars compacta (SNc), an area in the midbrain mainly targeting the striatum. This area is the main input gate of the basal ganglia, subcortical nuclei critical to managing motor behavior. Thus, a consistent reduction of striatal dopamine levels causes malfunctioning of the basal ganglia circuits that, in turn, may contribute to the emergence of tremor at rest (Pare et al., 1990;Wichmann and DeLong, 1999;Dovzhenok and Rubchinsky, 2012).
Based on this evidence, drug therapies for tremor often aim at recovering dopamine levels (Deuschl et al., 2000;Caligiore et al., 2019). However, while these approaches seem to produce amelioration for most PD motor dysfunctions, they generate variable responsiveness for resting tremor (Helmich et al., 2012;Wu and Hallett, 2013;Connolly and Lang, 2014). The lack of dopamine-based therapies consistency with tremor could be explained by considering that based on the dopamine response, PD resting tremor can be divided into various partially overlapping phenotypes that may be associated with different pathophysiological mechanisms (Zach et al., 2020). These mechanisms might be produced by a cortical-subcortical network involving areas traditionally not directly related to PD. Several data support this framework (Obeso et al., 2010;Caligiore et al., 2016Caligiore et al., , 2017bHelmich, 2018).
The membrane properties of the basal ganglia cells support pacemaking (Surmeier et al., 2005) but do not produce tremor oscillations in healthy basal ganglia circuits. In contrast in PD, tremor-related activity (i.e., neural activity in the tremor frequency band, correlated with the tremor movement) was observed in the globus pallidus (Hutchison et al., 1997) and subthalamic nucleus (STN) (Levy et al., 2000), in the thalamus (Thal) (Lenz et al., 1994), and in the cortex (Timmermann et al., 2003). Lesions in different parts of the basal ganglia-thalamocortical network [in the cortex (Deuschl et al., 2000), globus pallidus, and Thal (Mitchell and Ostrem, 2020), in the STN (Alvarez et al., 2005)] suppress tremor. The fact that breaking the loop at multiple sites leads to the same effect, tremor suppression, suggests that the loop itself, more than any of its parts, contributes to tremor generation.
The tremor network can also involve the dorsal raphe nucleus (DRN), which could modulate the release and concentration of dopamine through the serotonergic projections it sends to the striatum ( Bara-Jimenez et al., 2005;Di Matteo et al., 2008;Politis and Niccolini, 2015). It has been recently shown that aside from the dopaminergic dysfunction, in PD, there is a progressive loss of serotonergic terminals that has a slower pace and begins earlier than the dopaminergic one. This serotonergic dysfunction has been associated with the development of both non-motor and motor symptoms (De La Fuente-Fernndez et al., 2001;Lindgren et al., 2010;Politis and Niccolini, 2015) including tremor (Jankovic, 2018;Pasquini et al., 2018).
Overall, the vast range of conditions provoking or relieving tremor supports the involvement of a complex brain network, including both cortical and subcortical areas. Various dysfunctions within this network can produce different types of tremors. One mainly due to a direct dopaminergic system deficiency and for which dopamine-based actions work. Others involving the dopaminergic system and other brain areas, for which dopamine-based treatments are less effective (Pasquini et al., 2018;Zach et al., 2020). Starting from this system-level perspective, we propose a bio-constrained computational model to study the neural mechanisms underlying a possible type of PD tremor: the one mainly involving the serotonergic system. The model reproduces the interactions between basal ganglia direct (DP) and indirect pathways (IPs), primary motor cortex (M1) and Thal (Marreiros et al., 2013;Hintzen et al., 2018), and the critical role of the DRN, the brain region producing serotonin, on the basal ganglia circuits (Padovan-Neto et al., 2020). More in detail, the model reproduces the excitatory and inhibitory influences between these brain areas through an ordinary differential equations system (Reed et al., 2013). Moreover, it drives the movement of a simulated 2D anthropomorphic robotic arm, reproducing the main features of the human arm (Katayama and Kawato, 1993;Caligiore et al., 2014).
The simulations run with the model demonstrate that serotonin could affect dopamine levels to maintain homeostasis between the DP and IP signals (Di Giovanni et al., 2002;Guiard et al., 2008;Reed et al., 2013;Sgambato-Faure and Tremblay, 2018). When there is a dopaminergic dysfunction, the model shows that the serotonin-based compensation mechanism is not able by itself to bring the system back to physiological dopamine levels. This result agrees with recent studies (Roussakis et al., 2016;Jiménez-Sánchez et al., 2019). However, the model suggests that the compensation mechanism can nevertheless mitigate the progression of the disease in the very early stages, before the manifestation of overt tremor. The model also shows that serotonergic dysfunctions affect striatal dopamine release, indirectly contributing to the emergence of a serotonin-related tremor phenotype. Overall, these results could be critical to devising new early diagnosis strategies based on serotonin monitoring (Wilson et al., 2019).
The model allowed us to study the effect of a possible new drug therapy that builds on the serotonin-compensatory role to recover dopamine levels in simulated patients. This drug produces a boosting of the serotonin physiological compensatory effect that contributes to reduce tremor. This latter result has been validated by reproducing existing data collected with human patients, showing that a serotonin increase could produce a tremor decrease (Qamhawi et al., 2015).

Model Architecture
The model architecture reproduces the dynamical interaction between six components (Figure 1): M1, Thal, DRN, SNc, DP, and IPs. DP reproduces the overall activity of the striatum and the internal portion of the globus pallidus, whereas IP reproduces the activity of the circuit involving the external globus pallidus (GPe) and the STN (Smith et al., 1998;Gurney et al., 2001;Haber, 2003).

Model Equations
The following differential equations system simulates the activity of the model components and their dynamical interactions: (1) ). We obtained this system by critically manipulating the equations proposed by Reed et al. (2013). In particular, we introduced two innovations critical to investigate the role of serotonin in PD tremor: (i) a more detailed IP equation (Equation 8) to simulate the oscillatory behavior affecting the STN-GPe loop (Nambu et al., 2002;Stanford, 2003;Surmeier et al., 2005); (ii) a bio-constrained 2D dynamic arm, allowing FIGURE 2 | Representation of the dynamic arm used to test the model motor behavior. The task space is a horizontal plane in which the arm is free to move according to the inputs received from M1.
to record the effects of serotonin increase on overt tremor (Qamhawi et al., 2015). Although the reciprocal connections linking GPe and STN can produce and maintain low-frequency rhythms (≤5 Hz) (Plenz and Kital, 1999;Terman et al., 2002), they do not produce tremor oscillations in healthy basal ganglia circuits. In PD, the DA loss induces an excessive synchronization of the GPe-STN oscillatory activity, contributing to generate and propagate abnormal tremor-related oscillations over the striatalthalamo-cortical system (Nini et al., 1995;Guillery et al., 1998;Boraud et al., 2005;Gatev et al., 2006;Özkurt et al., 2011).
These equations allow studying the modulation between the various components of the model with simple linear functions, except for the oscillatory term (Equation 8), and for the term that represents the DA release in the striatum (Equation 6). Moreover, DA has an excitatory effect on the DP activity (Equation 7) and an inhibitory effect on IP (Equation 8). These two different types of dopaminergic modulations reflect what happens in the real basal ganglia (Cachope and Cheer, 2014;Fiore et al., 2015). Finally, each equation has a decay term, ensuring that the variable (firing rates or concentrations) goes to zero in the absence of input.

Dynamic Arm
The model controls a simulated 2D dynamic arm (Flash and Hogan, 1985) implemented with the realistic biomechanical parameters proposed in Katayama and Kawato (1993) (as shown in this study for more details on the arm equations) (Figure 2). We numerically integrated the dynamical arm equations using a Runge-Kutta fourth-order method with an integration step set to 0.01. The arm shoulder is anchored at We used a proportional derivative controller (PD controller ) to capture some muscular visco-elastic properties of the human arm (Caligiore et al., 2014). This mathematical model computes the torques vector (T) necessary to move the arm joints, starting from the desired angles supplied by M1 (A des ), the current angular position (A cur ), and the current angular velocity (Ȧ cur ): where K p is the vector containing the values of the PD proportional constant for the shoulder and elbow joints (shoulder = 20.0 Nm/rad; elbow = 10 Nm/rad), whereas K d is the vector containing the values of the PD controller damping constant (shoulder = 1.5 Nm/rad; elbow = 1.0 Nm/rad). These values were chosen taking into account human arm visco-elastic properties (An et al., 1981).

Simulation Settings
The model was developed using the Python programming language. The code is available here https://github.com/ctnlab/ serotonin_PD_tremor_model.  to those found in experimental observations (Feldman et al., 1997;Segovia et al., 1997;Jones et al., 1998;Knobelman and Lucki, 2001;Mahon et al., 2006;Ohara et al., 2007;Shimamoto et al., 2013) (see Table 2). In particular, the parameters settings were done through an automatic optimization procedure, namely a "genetic algorithm" (Fortin et al., 2012;Caligiore et al., 2017a). The algorithm went through optimization cycles called "generations". The parameters search started from the values used in Reed et al. (2013). For each generation, a "population" of candidate solutions corresponding to different possible sets of model parameters was "evolved" to (1) minimize the mean absolute percentage error (APE) between the model steady-state values and those drawn from the experimental observations. A low value of such error indicated that the activity exhibited by the model was similar to the data collected from the real subjects;

Tuning the Model Parameters
(2) minimize the oscillation in the healthy model (meanPhys); (3) maximize the oscillation in the damaged model (meanPark). Each "generation" explores 30 candidate parameter populations, and hence models (during 300 generations of the genetic algorithm), to get a meanAPE = 0.003 (fitness = 0.9969), a value that guaranteed a good matching between the simulated and the real data; meanPhys = 0.0053 (fitness = 0.9947) to guarantee a low value of physiological oscillation; and (meanPark = 0.0051) to predict oscillation amplitude observed in PD subjects (fitness = 0.9949). The parameters shown in Table 1 are the best ones obtained with this procedure.

Simulation Trials
The model equations are numerically integrated with an integration step t = 0.1 (i.e., one simulation step corresponds to 0.1 s of real-time). The simulation runs for 300 s (3,000 simulation steps). At the beginning of the simulation, the shoulder and elbow angles of the arm are set to the starting position (α = 0 • ; β = 0 • ). Within the simulation time, we monitor the changes in the model components (i.e., activity of brain areas and neuromodulators concentrations) in five-time windows that we called "simulation trials": the first one lasts 75 s and it is used to simulate the behavior of the health system; the second one (75, 125) s is used to study the effects of the DA lesion (occurring at 75 s); the other three remaining trials (the first two lasting 50 s each, the last 75 s) are used to investigate the effects of increasing doses of the simulated serotonin-based treatment.

Experiment Repetitions
We collected data from 20 different simulated subjects, obtained by running the model with distinct seeds of a pseudo-random number generator, which in turn causes different values of a noise signal acting on the model. In particular, we added to the initial values of the variables of the equations (Table 2), a small random number drawn from a normal distribution having mean zero and SD set to the initial value times, a number increasingly chosen in the range (0.01, 0.02) (0.01 for the subject 1, 0.02 for the subject 20, with an increasing step set to 0.0005). At the same time, we added to the values of the parameters of the equations (Table 1), a small random number drawn from a normal distribution having mean zero and SD set to the value of the parameter times, a number increasingly chosen in the range (0.01, 0.02) (increasing step still set to 0.0005). In this way, we exploited the capacity given by the simulation approach to producing many subjects with low cost to have higher statistical reliability of the results.

RESULTS
This section presents the results obtained with the model and directed to (i) support the 5-HT role in maintaining homeostasis between DP and IP by affecting DA release (Reed et al., 2013; Sgambato-Faure and Tremblay, 2018); (ii) demonstrate that 5-HT dysfunctions can indirectly contribute to the emergence of tremor by producing a DA decrease; (iii) study the effect of a new pharmacological treatment for tremor that acts on 5-HT to recover DA levels. This latter result was validated by reproducing data collected with human patients showing that 5-HT enhancement could reduce tremor (Qamhawi et al., 2015).

Serotonin Compensates Dopamine Loss
To investigate how 5-HT affects DA release, we introduced an SNc malfunctioning. At the beginning of the simulation, the model has no dysfunctions. After 75 s, we introduced an SNc impairment by increasing the value of τ SNc (+25%) (Equation 5). This action produces an SNc activity decrease approximating the effects of a neural loss. At the same time, we increase the value of α IP (+11%) (Equation 8), simulating an efficacy decrease of the DA receptors that regulate the oscillatory component of IP.
The SNc spike frequency decreases, leading to a reduction of DA concentration (Figures 3A,B). This loss produces a chain of effects: Thal activity decreases ( Figure 3C) because DP activity goes down whereas IP activity goes up ( Figure 3D); M1 activity decreases ( Figure 3C) because it directly depends on Thal activity. Lower activation of M1 produces a lower inhibition of DRN that, as a consequence, increases its activity ( Figure 3B). Finally, higher activity of DRN supports an increase of 5-HT release that, in turn, partially recovers the DA level ( Figure 3A).

Serotonin Dysfunctions Produce Dopamine Loss
Starting from a healthy system, we produced a 5-HT impairment by gradually increasing the τ DRN parameter, approximating the effects of a DRN neural loss (Equation 3). Figure 4 shows that as we increase the degree of the simulated 5-HT impairment, the DA concentration in the striatum drops down. As a result, the oscillatory component of IP increase (Equation 8), feeding the emergence of an overt tremor (Figure 5). This outcome supports the hypothesis that 5-HT dysfunctions might indirectly contribute to the emergence of tremor.
The results shown in Figures 4, 5 were statistically validated through a repeated-measures ANOVA with the condition (HEALTH, DAMAGE 1, DAMAGE 2, and DAMAGE 3) as a within-subjects factor. The Bonferroni post hoc test was also conducted on significant interactions, and the Greenhouse-Geisser correction was used in the case of violations of sphericity. The ANOVA revealed that there is a statistically significant difference between DA concentrations in different conditions [F (1,016−19,299) = 61,683,499, Cohen's f = 0, p < 0.001, power = 1]. Moreover, the post hoc analysis showed that the health condition (i.e., with no DRN damage) is significantly different from the other three conditions, with a greater DA concentration (level of DA concentration: 2,723 vs. 2,435 vs. 2,313 vs. 2,204 nM, respectively). Furthermore, the DA concentration decreases with the damage entity, the greater the DRN damage, the lower the DA concentration. Similarly, the ANOVA showed that there is a statistically significant difference between the tremor rate in different conditions [F (1,011−19,204) = 6,967,198, Cohen's fr = 0, p < 0.001, power = 1]. The post hoc analysis confirmed that the health condition is significantly different from the other three conditions, the tremor, in this case, is reduced with respect to the DRN damage conditions. Furthermore, the greater the DRN damage, the greater the tremor.

A New Treatment for Tremor Based on Serotonin Increase
In this section, we show that it is possible to build on the neural mechanisms underlying the physiological influence of 5-HT on DA release (sections 3.1, 3.2), to test a new treatment for tremor based on a 5-HT increase. Figure 6 shows the effects of this treatment at the level of neural activity, whereas Figure 7 shows the effects on tremor.
Starting from a healthy model, after 75 s we introduced an SNc impairment by increasing the value of τ SNc (+70%) (Equation 5) and by increasing the value of α IP (+333%) (Equation 8). Since these impairments are stronger than those produced  in section 3.1, they produce a stronger DA loss (compare Figures 6A,B vs. Figures 3A,B). In this case, the physiological compensatory mechanisms involving 5-HT in between 75 and 125 s only allow to partially recover the DA loss, and it is not strong enough to avoid tremor (case "SNc DAMAGE" on Figure 7).
Then we boost the physiological influence of 5-HT on DA release by gradually reducing the value of τ 5−HT (Equation 4) between 125 and 300 s. This parameter could be a proxy of the selective serotonin inhibitors (SSRIs) effects on serotonin transporters, which carry 5-HT in the DRN pre-synaptic neurons (Serretti et al., 2006;Dogan et al., 2008;Murphy and Lesch, 2008). Thus, a gradual decrease of τ 5−HT produces an increase of 5-HT release simulating the effects of a treatment based on a gradual increase of SSRIs. More in detail, we simulated three subsequent increases of SSRIs which produce the following chain of effects: Thal activity gradually increases ( Figure 6C) because DP activity gradually goes up whereas the IP activity goes down (Figure 6D); as a consequence, M1 activity increases ( Figure 6C) because it directly depends on Thal activity; a greater activation of M1 produces greater inhibition of DRN that, as a consequence, decreases its activity ( Figure 6B); a lower activity of DRN supports a lower inhibition of SNc (Equation 5) that, in turn, recovers the DA level reaching similar levels observed before the lesion (Figure 6A). The reduction of IP oscillatory activity together with the increase of M1 activity produce a gradual decrease of tremor (Figure 7). FIGURE 7 | Oscillation amplitude of the bio-mimetic arm before and after the DA impairment (τ SNc + 70%, α IP + 333%) and following SSRIs treatment (τ 5−HT − 20%, τ 5−HT − 30%, τ 5−HT − 40%).
Data were analyzed with repeated-measures ANOVA with the condition (HEALTH, SNc DAMAGE, TRMT 1, TRMT 2, and TRMT 3) as a within-subjects factor. The Bonferroni post hoc test was also conducted on significant interactions, and the Greenhouse-Geisser correction was used in the case of violations of sphericity. The ANOVA revealed that the main effect of SSRIs treatment [F (1,192−22,654) = 3,475,671, Cohen's f = 0, p < 0.001, power = 1] is significant. The post hoc test showed a significant increase of tremor oscillation after the damage (0.005 vs. 0.032 m), furthermore, the treatment improves the tremor, in particular, increasing the dose decreases the tremor (0.024 vs. 0.019 vs. 0.016 m).

DISCUSSION
In this study, we proposed a system-level computational model to investigate the role of the serotonergic system in PD tremors. The results shown in Figure 3 suggest that at the very early stage of the disease, a physiological increase in 5-HT concentration, which indirectly affects DA release, could partially recover a low reduction of DA. Even though this compensatory mechanism is not able by itself to bring DA back to physiological levels, it contributes to maintaining DA at values that do not produce overt tremor in the patient ( Figure 3D). Figure 4 further supports the critical role 5-HT has in the modulation of DA concentration by showing that 5-HT impairments produce DA loss and the emergence of tremor ( Figure 5). These results agree with empirical data showing that high levels of 5-HT could promote DA release from SNc dopaminergic projections (Blandina et al., 1989;Bonhomme et al., 1995;De Deurwaerdère et al., 2002). Interestingly, the tremor entity (amplitude of oscillations) produced by the 5-HT dysfunction is lower than the tremor entity obtained by the impairment of the dopaminergic system ( Figure 5 vs. Figure 7). The 5-HT dysfunction, indeed, only acts on SNc DA release. By contrast, the tremor entity shown in Figure 7 is due to both the SNc DA release dysfunction (τ SNc ) and the efficacy reduction of D2 DA receptors that regulate the oscillatory component of IP (α IP ). For this reason, the 5-HT-based treatment only partially contributes to recover the IP health activity ( Figure 6D). Hence, the model suggests that 5-HT impairments alone contribute to the emergence of weak overt tremor symptoms, like the one that one could show at the very beginning of the disease progression. To have large tremor oscillations, it is also critical to consider the IP D2 DA receptors responsiveness dysfunctions. This result agrees with data we obtained using a spiking neural network model to study the differences between tremor and patients with akinetic PD  and agrees with experiments suggesting the involvement of both 5-HT and DA dysfunctions in the progression of resting tremor (Pasquini et al., 2018).
The results summarized in Figures 6, 7 demonstrate the effectiveness of a possible treatment based on the gradual increase of SSRIs on 5-HT transporters carrying 5-HT on DRN (Serretti et al., 2006;Dogan et al., 2008;Murphy and Lesch, 2008). Figure 6 shows that the physiological increase of 5-HT concentration operating before the simulated treatment (75, 125) s produces a compensatory mechanism that increases the DRN activity to partially recover DA levels. By contrast, the treatment (125, 300) s leads to a DRN reduction, producing a higher DA level recovery and a tangible reduction of overt tremor (Figure 7).
Interestingly, the simulations show that the SSRIs treatment triggers different network dynamics than those involved in the 5-HT physiological compensatory mechanisms operating before the treatment. The physiological compensation could be biologically explained by a negative feedback mechanism mainly due to increased activation of the 5-HT1A autoreceptors present in the axonic terminals of 5-HT DRN neurons (Hajós et al., 1995). This mechanism produces an increase in DRN activity ( Figures 3B, 6B up to 125 s). By contrast, the simulations run with the model to test the effect of the SSRIs treatment suggest that the negative feedback between the concentration of 5-HT and the DRN activity has to follow a different dynamical trajectory, producing a decrease in DRN activity to get a high DA level recovery. The reduction of the 5-HT degradation constant (τ 5−HT ) causes an increase in the 5-HT extracellular levels, which leads to an increased DA release by SNc. The effect is an excitatory modulation of DP and an inhibitory modulation of IP. This dynamic produces an increase in the firing rate of Thal and, in turn, of M1 with consequent reduction of the DRN activity.
Parkinson's disease tremor treatments based on dopamine medication (e.g., Levodopa administration) lead to a widespread and non-specific increase of DA concentration in the SNc (Connolly and Lang, 2014;Lee and Yankee, 2021). By contrast, the treatment with SSRIs proposed in this study would allow acting in specific ways on the area of interest. In this case, we increase the intrinsic concentration of 5-HT. By directing the drug delivery toward specific cellular markers (Huang et al., 2019), it is possible to inhibit only the SSRIs of our interest, allowing a more focused and not widespread increase in neurotransmitter concentration (Valentino and Commons, 2005;Wang et al., 2005;Rao and Nanda, 2009;Dankoski et al., 2016).
Despite these encouraging results, the current version of the model mainly represents a tool to test the effects of neurotransmission modulation rather than test new treatments based on 5-HT manipulation. Recent influential data demonstrate that 5-HT loss could precede PD motor symptoms and could be critical to study the PD pathophysiology and to develop novel therapeutic strategies (Wilson et al., 2019). However, evidence obtained from clinical research found several hurdles that, to date, have impaired the clinical development of serotonergic drugs (Huot et al., 2017;Muñoz et al., 2020). For example, the studies conducted in patients with PD and animal models mainly focussed only on 5-HT1A and 5-HT2A receptors, but the 5-HT system consists of many receptors (at least 14). It is not yet clear, what is the role of these other receptors. With 5-HT1A stimulation, it seems not easy to separate any benefit from interaction with L-DOPA anti-Parkinsonian action. In addition, the use of 5-HT1A receptor partial agonists also seems marred by a potentially detrimental effect on L-DOPA anti-Parkinsonian action, suggesting that such drugs may have a narrow therapeutic window (Bezard et al., 2013).
Hence, modulation of the 5-HT system is a promising way to address several manifestations of PD, but further investigation is required to clarify mechanisms of neurotransmitter interactions and to determine optimal compounds and doses for effective therapies producing the maximal benefit with minimal adverse events (Huot et al., 2017;Muñoz et al., 2020). New data from clinical research will support the design of an updated version of the model to make more consistent the simulations of the potential role of 5-HT for early diagnosis and therapy.

CONCLUSIONS
This article proposes a computational model to study how serotonin could affect dopamine release within the basal ganglia and how serotonin dysfunctions indirectly contribute to the emergence of tremor. The computer simulations run with the model suggest that serotonin changes monitoring could help in early PD diagnosis. They also demonstrate the effectiveness of possible new pharmacological treatments for tremor acting on serotonin to recover dopamine levels. However, evidence obtained from clinical research found several hurdles that, to date, have impaired the clinical development of serotonergic drugs (Huot et al., 2017;Muñoz et al., 2020). Therefore, more clinical evidence is needed to confirm or disprove the results obtained with the model.
The work focuses on the role of serotonin in Parkinsonian tremor, for this reason, studies how serotonin affects dopamine release through the meso-striatal circuits, including the SNc. Other studies on depression, mild cognitive impairment, and Alzheimer's disease, instead, have investigated the interaction between serotonin and dopamine release through circuits involving the meso-cortico-limbic network, including the ventral tegmental area (Martorana and Koch, 2014;Silvetti et al., 2019;Wang et al., 2019;Caligiore et al., 2020). Future studies on the relationship between these two different ways the serotonin affects dopamine could be critical to understanding the systemlevel neural mechanisms underlying the comorbidities between PD, depression, mild cognitive impairment, and Alzheimer's disease (Marsh, 2013;Caligiore et al., 2016Caligiore et al., , 2017bCummings et al., 2019).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://github.com/ctnlab/ serotonin_PD_tremor_model.

AUTHOR CONTRIBUTIONS
DC: conceptualization, formal analysis, investigation, methodology, resource, validation, software, visualization, writing -original draft, writing -review and editing, funding acquisition, project administration and supervision. FM: data curation, formal analysis, investigation, writing -review and editing, visualization, validation and software. SB: formal analysis, investigation, methodology, resource and software. AC: data curation, formal analysis, investigation, methodology, resource, writing -review and editing, visualization, validation and software. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Flora Giocondo for helping in statistical analysis.