Simulating combined monoaminergic depletions in a PD animal model through a bio-constrained differential equations system

Introduction Historically, Parkinson's Disease (PD) research has focused on the dysfunction of dopamine-producing cells in the substantia nigra pars compacta, which is linked to motor regulation in the basal ganglia. Therapies have mainly aimed at restoring dopamine (DA) levels, showing effectiveness but variable outcomes and side effects. Recent evidence indicates that PD complexity implicates disruptions in DA, noradrenaline (NA), and serotonin (5-HT) systems, which may underlie the variations in therapy effects. Methods We present a system-level bio-constrained computational model that comprehensively investigates the dynamic interactions between these neurotransmitter systems. The model was designed to replicate experimental data demonstrating the impact of NA and 5-HT depletion in a PD animal model, providing insights into the causal relationships between basal ganglia regions and neuromodulator release areas. Results The model successfully replicates experimental data and generates predictions regarding changes in unexplored brain regions, suggesting avenues for further investigation. It highlights the potential efficacy of alternative treatments targeting the locus coeruleus and dorsal raphe nucleus, though these preliminary findings require further validation. Sensitivity analysis identifies critical model parameters, offering insights into key factors influencing brain area activity. A stability analysis underscores the robustness of our mathematical formulation, bolstering the model validity. Discussion Our holistic approach emphasizes that PD is a multifactorial disorder and opens promising avenues for early diagnostic tools that harness the intricate interactions among monoaminergic systems. Investigating NA and 5-HT systems alongside the DA system may yield more effective, subtype-specific therapies. The exploration of multisystem dysregulation in PD is poised to revolutionize our understanding and management of this complex neurodegenerative disorder.


Introduction
Common theoretical and empirical approaches studying Parkinson's disease (PD) focus on dysfunctions in dopamine (DA)-producing cells in the substantia nigra pars compacta.This area projects to the striatum, the principal input gate of the basal ganglia, subcortical nuclei critical to managing motor behavior (Tozzi et al., 2021;Ledonne et al., 2023).Thus, a consistent reduction of striatal DA levels causes malfunctioning of the basal ganglia circuits that, in turn, contribute to the emergence of different PD symptoms (Pare et al., 1990;Dovzhenok and Rubchinsky, 2012;Caligiore et al., 2019).The main PD motor symptoms include resting tremor, bradykinesia, rigidity, and freezing of gait (Jankovic and Kapadia, 2001;Obeso et al., 2010;Caligiore et al., 2019).Cognitive impairments might be evident at diagnosis, even though they significantly manifest in the later stage of the disease progression (Williams-Gray et al., 2007;Aarsland et al., 2009).Moreover, several recent studies suggest that psychiatric disorders, such as depression or anxiety, often develop several years before typical motor symptoms (Faivre et al., 2019).In particular, motivational system dysfunctions manifest early in PD (Pagonabarraga et al., 2015;Cohen et al., 2022;Favier et al., 2022).
Starting from this system-level perspective, we propose a bio-constrained computational model that, for the first time, explicitly investigates the neural mechanisms underlying interactions between DA, NA, and 5-HT in a PD animal model.Data on basal ganglia and monoamine areas physiology (Kang and Kitai, 1993;Szabo and Blier, 2001;Liu et al., 2002;Damodaran et al., 2014), including evidence on the effects of NA and 5-HT depletions in a PD animal model (Delaville et al., 2012), constrain key model parameters, hinting at a potential causal dynamical interaction between basal ganglia regions and the areas responsible for neuromodulators release.The model produces predictions on expected activity changes in other brain areas which are not investigated in the target experiments of Delaville et al. (2012).
It also highlights the potential efficacy of alternative treatments targeting LC and DRN, yet these are preliminary findings whose effectiveness needs further validation.A sensitivity analysis of the model parameters allows us to identify the critical features affecting the model activity.In this way, it could be possible to frame the most important features affecting the mechanisms underlying the activity of simulated brain areas.Finally, we made a stability analysis that confirms the soundness of the mathematical formulation used to design the model.This point could be critical to validate the effectiveness of the model (Fornari et al., 2020;Shi et al., 2020).It is a relevant aspect that is not readily apparent and seldom explored within computational neuroscience literature.While the soundness of the mathematical formulation is crucial, it becomes irrelevant if the model lacks robust theoretical foundations.Thus, both aspects, mathematical integrity and theoretical underpinnings, serve as indispensable pillars for developing a computational model.
The system-level computational model proposed in this article emphasizes that PD is a multifactorial disorder and opens promising avenues for early diagnosis and subtype-specific treatment tools, harnessing the intricate interactions among monoaminergic systems (Marras et al., 2020;Caligiore et al., 2021;Severson et al., 2021).The investigation into multisystem dysregulation in PD is poised to profoundly transform our comprehension and handling of this intricate neurodegenerative condition.

Methods . Dynamical model
We introduce a novel model designed to explore the dynamic interactions of three prominent brain monoamines: dopamine (DA), noradrenaline (NA), and serotonin (5-HT), within both typical and pathological contexts.Figure 1 provides an overview of the neural circuitry we assessed.We mimic monoamine efflux by simulating the activation of critical brain regions responsible for initiating such release.Regarding 5-HT/DRN, it is reported that DRN sends projection to striatum (Vertes, 1991;Miyanishi et al., 2023), and we model these projection as excitatory on both direct and indirect pathway.DRN to GP is modeled as excitatory due to serotonin increasing the firing rate of GP neurons as reported in the literature (Chen et al., 2008;Rav-Acha et al., 2008).The projections from DRN to SNcVTA were reported by Gervais and Rouillard (2000) to be dense, and the authors underline the possibility of this connection to be mainly inhibitory; there is a possibility for this projection to be both excitatory and inhibitory but future studies are required to investigate this matter further.Here, we assume this connection to be inhibitory as previously reported by several studies (Dray et al., 1976(Dray et al., , 1978;;Kelland et al., 1990Kelland et al., , 1993)).The projection from DRN to LC is also inhibitory (Segal, 1979).StrD1/D2 has inhibitory effects on GP (Nicholson and Brotchie, 2002), while SNcVTA has inhibitory effects on strD2 and excitatory effects on strD1 (Reed et al., 2013); moreover, SNcVTA has excitatory effects on LC as suggested in the literature (Deutch et al., 1986;Lin et al., 2008).The connection from SNcVTA to ./fncom. .
DRN is modeled as inhibitory because it is reported that after DA depletion, there is an increase of spontaneous activity in DRN (Wang et al., 2009).The schema from Figure 1 implies the following system of equations:

StrD2
(3) where the abbreviated notation ẋ stands for dx dt and: • the time constants τ x are all positive and refer to a dampening term which brings back the activity of each area to its resting activation level in the absence of external stimulation (see Table 1); • the parameters α represent the linear components of the system, are all positive, and follow the notation: α from to ; α ext x are synthetic terms that implicitly account for the rest activation of each area and other external stimuli which are not part of the modeled circuit; β is also positive and follows the same notation β from to but accounts for non-linear effects; • the ratios of monoaminic projections from an area to its targets are assumed to be constant.

. Formalization
Let y be the status vector of the system of Equations 1-6; we also define s to be the size of y, hence the number of equations in the system.We therefore have The system is autonomous and can therefore be represented in the form: where each component of the function f :(R × R s ) → R s is defined by the corresponding equation in Equations 1-6, and where we assume the initial state y(t 0 ) = y 0 to be known.Starting from Equation 8, Equations 1-6 can be represented in matrix form: where a ij = 0 if the corresponding α j i is not defined and likewise b i = 0 if α ext i is not defined, and also c ij = β j i with c ij = 0 if the corresponding β j i coefficient is not defined; "•" indicates the element-wise vector product (or Hadamard).

. Modeling lesions
Each element within the status vector (Equation 7) corresponds to the average activation frequency of the corresponding brain region.This, in turn, serves as an indirect indicator of the production and projection of monoamines to the affected areas.We posit that a depletion in monoamines results from the death or temporary incapacitation of a portion of neurons within a given area, which is directly manifested as a decrease in the average activation frequency of that area.Consequently, when the SNcVTA is lesioned, we observe a reduction in DA levels; lesioning the DRN results in a decrease in 5-HT, and lesioning the LC leads to a reduction in NA levels.
Each equation of the model is composed by three conceptual blocks: a damping term, a constant stimulus, and a reaction to projections from other areas.The constant stimulus represents external and internal activation sources that are not directly accounted for in this model.Together with the damping term, the constant stimulus accounts for the resting behavior of the area: The area will stabilize to its rest activation frequency.In absence of reaction terms, each equation has an equilibrium point: The time constants τ are derived from the literature (see Table 1), and we assume them to be typical values for the specific kind of neuron found in an area; we therefore assume that they are not altered by the lesion.It is, however, reasonable to expect that the sensitivity of lesioned area to internal and external stimuli will change in such a way that the average activation frequency changes to the levels which have been experimentally measured.
We can now define multiple versions of the same model, which differ from the healthy model only for the constant term and reaction coefficients of the lesioned area.For example, suppose a healthy subject is modeled using model (Equation 9) by the coefficients held in A, C, b.Having received a dopaminergic lesion (hence, SNcVTA neurons are malfunctioning), the subject will now be modeled by the same equations of model (Equation 9) but this time with coefficients A LDA , C LDA , b LDA , which differ by A, C, b only by the values corresponding to the parameters of the equation for SNcVTA, namely, α DRN SNcVTA , α LC SNcVTA , β LC SNcVTA , α ext SNcVTA .Likewise, when the subject also receives a serotonergic lesion, there will be a third set of parameters A LDA+L5HT , A single subject is therefore represented by multiple versions of the parameters matrices A, C, b, each set corresponding to one particular state: healthy (also called SHAM), LDA, L5HT, LNE when only one of the lesions is applied, LDA + L5HT, LDA + LNE when lesions are combined, and so on.

. Stability conditions
The solution trajectories of a dynamical system are referred to as "stable" when small perturbations of the initial conditions lead to trajectories which have fundamentally the same behavior of the unperturbed ones and differ from the latter in a proportional way with respect to the the perturbation magnitude.An unstable system can otherwise have a high sensitivity to such perturbation, which originate wildly varying solution behavior.At the limit of instability, there are chaotic systems, for which very small perturbation of the initial conditions can lead to completely different dynamical behavior.Stability is a property of a particular solution (Kelley and Peterson, 2001;Lakshmikantham and Trigiante, 2002;Riley et al., 2006;Press, 2007;Butcher, 2016).
A system which has stable solution trajectories can have great predictive power since its low sensitivity to the initial conditions results in solutions which have comparable errors.Chaotic systems, on the other hand, are so sensitive to the errors on the initial conditions that it is effectively impossible to use them for obtaining useful predictions.For instance, weather is a chaotic system, and that is why it is so difficult to compute reliable long-term weather predictions.Fortunately enough, the aspects of brain chemistry that we are simulating in this study do not exhibit chaotic behavior; on the contrary, they exhibit self-regulation and great resilience to perturbations.It makes therefore sense to require them to be simulated using a system of equations having asymptotically stable equilibria.
Equation 9 has at least one equilibrium point ȳ such that: The stability properties of the equilibrium point in Equation 11 are analyzed in details in Appendix 1; the resulting stability conditions are employed as one of the components of the fitness measure of the model, as described in Section 2.7.In this way, it is ensured that all the solution trajectories simulated by the model are asymptotically stable.Consequently, the numerical trajectories exhibit the expected behavior even for long-time simulations.

Simulation setup . Free parameters and constants
As hinted in Sections 2.3 and 3.3, we require a model to be able to reproduce its target data in four different states at the same time: healthy (SHAM), dopaminergic lesion (LDA), noradrenergic lesion (LNE), and serotonergic lesion (L5HT).The combinations LDA + LNE and LDA + L5HT are instead constrained only to a target range, to be able to also serve as a prediction (and hence as a measure of the agreement of the model with experimental data).Below, we briefly outline the mathematical formalism used to assess the model's stability properties and to establish the simulation setup (Equations 12-18).
Let S i be the set of parameters that define the model representing test subject i. S i contains the following: The set S i therefore contains a total of 36 parameters, 30 of which must be optimized at the same time to fit the available data.Appropriate subsets of the parameters in S i are then used to build the corresponding matrices A, C, b to completely define (Equation 9).and hence compute its solution and properties.

. Model fitness definition
We will hereafter refer to S i as the complete model for subject i, since it is the set of parameters that completely define it.The variations S kind i , such as S SHAM i and S LDA i , will refer instead to the subset of parameters which are currently being applied to actually simulate the model.
We will denote one solution as Y is therefore the solution obtained by integrating the model in the interval [t 0 , T], with the starting vector y 0 , and using the SHAM subset of parameters.The matrix Y comprises s rows, and the i-th row corresponds to the i-th equation in the system of differential equations.
The number N of integration steps, as well as their size, is usually variable and chosen by the integration method case-by-case; hence, it can potentially be different for each subset of parameters.
Likewise, we will denote with T kind i the corresponding reference solutions that will be used to evaluate the fitness of the model: where J = [t 0 , ..., t n ] is a vector of times.The fitness of a model is finally obtained as a combination of many fitness figures which measure a wide range of properties of the simulated solution with respect to the reference ones.In more detail: Details about the fitness computations are reported in Appendix 2.

. Synthetic reference data
The average brain activation in healthy subjects, as well as the time constants, has been extracted from the literature and are reported in Table 1.Starting from these average values, we craft a synthetic activity profile for a population of subjects (240 simulated subjects).This is achieved by modeling a normal distribution centered around the average value, with a normalized maximum excursion of ±50%.
In the reference (or target) data that we aim to replicate using the computational model, a population of adult male rats is subdivided into six groups: All SHAM SHAM L5HT LNE LDA LDA L5HT LNE where: • SHAM: indicates that subjects are treated with saline;  The experimental data we have collected, summarized in Tables 1, 2, ultimately consist of normal distributions around their respective center values which are to be considered constant; we do not have explicit information about the dynamic behavior of the system or about the transition from a state to the other.Of course, each subject must go through a dynamic transition from a healthy state to a lesioned state, but both states must be asymptotically stable solutions for the model.Since there is no single study that lists all the required brain area activation values for a particular subject at the same time, we have no choice but to generate a synthetic population of virtual subjects with area activation values which lie within the distributions identified across the literature.The generated target value distributions for all cases are summarized in Table 3 and Figure 2.

Area
Value In particular, we generate a number of subjects S i and each of them is associated with a set of target values T i which are defined as described in Table 3, where x ← N(µ, σ ) indicates that x is a random number drawn from a normal distribution using the provided parameters.Given the synthetic nature of this data, it is reasonable to assume activities of each area to have the same distribution.We impose every value to lie within ±50% of the center value by imposing 4σ = 1 2 µ.Since data from literature can be interpreted as an average percentage change in activity for lesioned areas, we decided to treat the generated healthy value for an area as the reference level of an individual and use that as a base to generate the lesioned values that where needed.In this way, a subject that has a higher than average value for an area in healthy conditions will also have an higher than average level in the same area when lesioned, although the value will change by the required proportional amount.
Reference solutions for subject i are therefore composed of constant values: where d is one of SHAM, LDA, L5HT, LNE, LDA+L5HT, LDA+LNE; the appropriate value of each area for the respective lesion is chosen according to d when available, otherwise defaults to the SHAM (not labeled) value. .

Choosing the simulation time
The simulation time is arbitrarily set to 0.5s under the assumption that the basic behavior of each equation in the system will resemble (Equation 10); hence, the transition time between any state to a stable solution will be dominated by the slowest time constant (which is derived from the literature).In fact, since the solution to: assuming y(0) = 0, is we can compute the time it takes for the solution to grow past 99% of its limit value: (which in engineering contexts is broadly known as the "rule of the five taus").
In this specific case, the slowest τ ≤ 20ms, therefore we can assume the transient phase to be finished after 5 • 0.02s = 0.1s, and a time 5 times longer, 0.5s, should be adequate to see a long stable steady state, and we would expect the dynamic behavior to have stabilized already around 0.1s.

. Model parameters dimensionality analysis
Each component of the status vector y (Equation 9) directly represents the average activation frequency of a brain area and is therefore expressed in Hz.The derivative terms in each equation of Equations 1-6 are all derivatives with respect to time of a frequency; hence, they are all expressed in Hz/s (or 1/s 2 ).Consequently, the external stimulus parameters α ext must also be expressed in Hz/s, while the remaining α parameters must be 1/s, hence Hz.The second order term parameter β is instead a pure number, since Hz 2 =1/second 2 =Hz/s.Finally, all time constants τ are naturally expressed in seconds.

. Simulation and optimization
All simulations are obtained using a variable order, variable step scheme backward-differentiation formulas (BDF) (Byrne and Hindmarsh, 1975) integrator provided by the Python scipy library.The optimization of parameters is then performed using differential evolution (DE), also as provided by the Python scipy package, with a population of 240 simulated subjects, DE/best/1/exp strategy with C r = F = 0.95 initialized with a uniform Halton distribution.
An external optimization cycle which sets different random generator seeds has been used to retry the cases which did not find convergence.In so doing, all cases did eventually converge after a few attempts.The full codebase can be found at https://github.com/WohthaN/Simulating_noradrenaline_and_serotonin_depletions_in_parkinson.

Results
All models fit the corresponding target values as defined in Table 3 with a fitness f ≥ 1 − 10 −8 .Since the measure is dominated by the smallest fitness value being combined by definition, also the mean square difference of each component of the solution from its reference value is bonded by the same order of magnitude, which is a suitable precision for the purposes of this study.Figure 3 shows an overview of the simulated behavior of all six areas in all conditions.The model reproduces the target data with animals presented in Delaville et al. (2012).In all cases, the simulated values for the SHAM case overlap the target values with the imposed tolerance.In the lesion groups, all areas which do not have a target defined are model predictions.All instances of the model meet the stability conditions defined in Section 2.4 as imposed by the fitness measure.
Figures 4, 5 give more insight in the behavior of each modeled area, also offering a comparison between the distributions over the whole population (which unfortunately could never be measured in vivo) and the more realistic sampled population distributions.The left graphs show distributions over the whole population (synthetic result), while on the right populations are sliced to have each subject in only one group (as in Figure 3).This latter case reproduces real laboratory conditions where one subject can only be measured once and hence belongs to one group only.Values for SHAM, LDA, L5HT, and LNE are fitted exactly on the required value, while the respective combinations are instead predictions of the model, whose values are only rangeconstrained with the rules defined in F COMB (see Equation 55 in Appendix 2.5).

. Sensitivity analysis
We used the computational model to gain insights into the critical parameters that govern the activity of simulated brain regions.Specifically, we conducted a sensitivity analysis for each simulated brain area in relation to the various model parameters.Each parameter can be varied independently (within some acceptable range) to record its effects on the simulated brain areas.Similar observations can be replicated for all individuals of the available population, and the average excursion of each area can then be compared with the average excursion of the parameter to infer a sort of "parameter importance".Below, we briefly discuss the mathematical formalism used to conduct the sensitivity analysis (Equations 19, 20).
In particular, for each individual of the population and for each parameter, we: • Vary the parameter around its original value ±50% in 100 uniform steps: if v is the parameter value for individual S SHAM i , we produce the set V param,i = {( x 99 + 97 198 )v} ⊆ [0.5v, 1.5v], i = 1, ..., 100.We therefore have a V param,i set for each parameter of each individual.
• Simulate the model using each value in V param,i and save the final value for each brain area.If the simulation stops early (hence the simulation diverged or reached physically impossible states), the result is discarded and the parameter value removed from V param,i .For each parameter and each individual, we therefore obtain six sets, one for each brain area, which we call A area param,i .
The sets are then joined across the population: where param ∈ S SHAM is one of the free parameters as defined in Section 3.1, area is one of the six brain areas {GP, StrD1, StrD2, SNcVTA, DRN, LC}, and i points to the i−th subject in the population.A sensitivity index is then computed for each area by scaling both V and A by their respective median values and dividing the standard deviations: I param,area can naturally be seen as a sensitivity matrix, with one column per area and one row per free parameter.As a last step, I param,area is normalized with respect to its maximum value.Figure 6 shows the computed sensitivity matrix for the entire fitted population in the SHAM case; a value of 1 indicates the maximum measured sensibility, while a value of 0 would mean that a particular parameter has no effect on that area.It is evident from the matrix that all the areas are relatively sensitive to changes in noradrenalinergic balance (external activation of LC), and also, even if in a somewhat lesser extent, to changes in the serotonergic balance (external activation of DRN).It is therefore reasonable to expect the stimulation of LC and/or DRN to produce changes in the activation of all areas.

. Statistical analysis
Before performing the statistical analysis, all data were checked to verify normal distribution applying D'Agostino and Pearson's normality test, as provided by the scipy package.When appropriate, one-way ANOVA by post-hoc test using Tukey's honestly significant difference (HSD) (also as provided by the scipy package) was performed.Histograms are annotated according to the p-value of the HSD test as follows: * * * * ≤ 0.0001 < * * * ≤ 0.001 < * * ≤ 0.01 < * ≤ 0.05.

. Possible roads to a treatment
We used the model that accurately replicates existing data to forecast potential outcomes in scenarios where experimental measurements have not yet been conducted.Our primary emphasis was on exploring alternative treatments centered around the manipulation of monoamines.The robustness of the model data reproduction, while rigorously adhering to stability criteria, bolstered the strength of our predictions.
Comparing the brain area activation level distributions in the SHAM to the LDA groups (Figures 7, 8), it is evident that the dopaminergic depletion also inhibits DRN and hence provokes a statistically significant serotonergic depletion.This behavior is compatible with the serotonin measurements reported in Delaville et al. (2012).The administration of LDA, however, does not alter significatively the behavior of LC (and hence noradrenaline   production).In the context of dopaminergic depletion, particularly due to a lesion in the SNcVTA caused by LDA, we engage in simulations aimed at exploring strategies to alleviate activity dysfunctions by targeting alternative monoaminergic circuits.
According to the model schema in Figure 1, dopaminergic levels can potentially be altered in two ways: 1. Externally stimulate the LC to change its production of noradrenaline.2. Externally stimulate the DRN change its production of serotonin.
The stimulation could either be chemical, by providing the area of the precursors needed to generate monoamines, or electrical, to artificially alter the average firing rate of the neurons from that area (and hence producing and projecting more monoamines to the areas which receive projections from the stimulated one).According to the sensitivity matrix in Figure 6, although, it is reasonable to expect LC stimulation to be strongly influential on dopamine levels, but DRN stimulation should have a smaller effect on the activation of the SNcVTA and strong side effects instead, which would not be compatible with a successful treatment.

. . Treatment optimization
Whether the stimulation of LC, DRN, or both could potentially restore healthy levels of brain areas in depleted subjects can be verified through the optimization of a subset of the parameters of our model.In particular, we can try to optimize the external stimulation parameter of LC, DRN, or both in the LDA version of our simulated subjects.First of all, we need to extend a subject set of parameters S i , as previously defined in Section 3.1, with three new subsets of parameters, namely: • LDA + cLC with free parameters α ext LC and corresponding fitness measure The corresponding model matrices A, C, b are constructed by using as base the LDA (hence, dopamine depleted) set of parameters for a subject and leaves as the only free parameters the external stimulation of the areas being tested.The three populations need to have different fitness measures because they each stimulate a different area to simulate the treatment.The stimulated area must of course be ignored by the respective fitness measure.Let us examine in detail the measure F cLC .Similarly to the composed fitness measure described in Appendix 2.1 for the healthy model, this measure is defined as the composition of the following measures: • The mean square error of the area activation value, one measure per area, as defined for the SHAM case in Equation 38in Appendix 2.2, but excluding the area being stimulated (in this case, excluding LC).• A parameter constraint similar to the one defined in Equation 56in Appendix 2.6, but this time used to enforce the external  stimulus parameter to be equal or greater than the original once (hence forcing the optimization to choose a stimulation rather than an inhibition).In particular, the component is defined as in Equation 21: • An asymptotic stability constraint as defined in Equation 58in Appendix 2.7, where of course Ã is constructed using the current parameters subset S LDA+cLC .
The fitness measures F cDRN and F cCOMB are of course constructed in an analogous way.In the latter case, mean square errors for both stimulated areas are ignored in the measure.The optimization is finally performed independently on all subjects of the three groups using the same algorithm described in Section 3.6, including the outer optimization cycles.

. . Treatment e cacy
The optimizer could successfully restore healthy levels of the measured areas in the vast majority of subjects by stimulating LC or both LC and DRN, but it never succeeded by only stimulating DRN.In particular, there was convergence in 197 subjects (which we refer here as responders) and did not find a satisfactory solution in the remaining 43 subjects (which we refer here as non-responders).Figure 9 shows that in the combined treatment, which obtained very similar results to the stimulation of LC alone, the relative increment to the external stimulation parameter of DRN is in fact several orders of magnitude smaller than the one applied to the corresponding parameter for LC.We can therefore assume that while the combined stimulation may have resulted in a slightly better fitness from the purely numerical perspective, DRN stimulation is indeed not useful as a treatment also in combination to LC stimulation.
Figure 8 provides compelling evidence of the profound impact of statistically significant LC stimulation in reinstating the equilibrium of serotonin and dopamine levels, as reflected in the activation levels of DRN and SNcVTA, respectively, within the population of LDA subjects.
Figures 10-12 illustrate the changes of distributions in the parameter space lesion and the subsequent treatment.The right side of Figure 12 highlights the differences in parameter distributions between the subjects that have been successfully treated (in green), and the ones whose levels could not be successfully restored (in red).None of the parameters of the responder subjects are significantly different from the one of the non-responder ones.The only parameter that shows a small significance difference is the sensitivity of SNcVTA toward noradrenaline from LC, as is shown in Figure 11.However, the spread of the distribution of that parameter is very large, and the value of that specific parameter alone is not useful for predicting if a subject is a responder or not.An accurate statistical study of the parameter space would be necessary to determine whether a particular combination of parameters could be used for predicting the curability of a subject, but this lies outside the scope of this study.
Figure 12 also shows the parameter space comparison of responder (labeled "cured") and non-responder (labeled "not cured") individuals for all parameters.The treatment modifies the external stimulation to LC and DRN but is otherwise identical to the LDA case.This result agrees with Figure 9 demonstrating that there is a very small shift in the DRN stimulus distribution which is not appreciable in Figure 12.

Discussion
PD is a global health concern, impacting an estimated 10 million individuals worldwide (Balestrino and Schapira, 2020).Regrettably, a definitive treatment for PD remains elusive.However, there is optimism in the research community as numerous novel drugs are presently undergoing clinical trials.One of the pivotal facets of PD is its intricate association with various neural circuits, resulting in consequential alterations in the brain neurochemistry (Obeso et al., 2010;Caligiore et al., 2016Caligiore et al., , 2017;;Helmich et al., 2018).Our research endeavor confronts the challenge of exploring the interplay among these distinct circuits and how the monoamine system adapts during the progression of this condition.
To tackle this challenge, we employed a bio-constrained differential equation model of PD.This model enabled us to investigate how different brain regions respond to the individual or combined depletion of dopamine (DA), serotonin (5-HT), and noradrenaline (NA).Initially, we harnessed the model to replicate data obtained from Delaville et al. (2012).The simulation results remarkably mirrored the observed alterations in the firing  Relative increment applied by the optimizer to the external stimulation parameter of LC and DRN in the combined case, limited to the subjects that were successfully treated.The increments to the DRN stimulation are several orders of magnitude smaller than the ones applied to LC.The sole stimulation of DRN is not a viable treatment.The small changes applied to the DRN stimulation by the optimizer may therefore have contributed to a numerically better solution, which is, however, not substantially di erent from the one obtain by the sole stimulation of LC.The subjects which did not reach a fitness of (and hence are not to be considered successfully treated) have been excluded from this plot.The color scale gives an indication of the final fitness reached by the optimizer, blue is the lowest and yellow the highest.activity of the globus pallidus (GP) following 5-HT lesions in both control and dopamine-depleted rats.Notably, the model also yielded predictions of other deviations in firing activity within brain regions not examined in the original experimental setup involving rats.Ultimately, our model served as a valuable tool to simulate the effects of prospective therapies aimed at restoring the activity in the lesioned regions.To model the impairment of DA function associated with PD, we implemented a reduction in the activity of the SNcVTA (Berretta et al., 2005).As a result, we have a dysregulation of the baseline activity within the specific brain regions we examined.Notably, DRN exhibited a notable increase in firing activity (see Figures 3, 5) in line with the literature (Kaya et al., 2008;Zhang et al., 2018;Caligiore et al., 2021).The DRN intrinsic excitability (Prinz et al., 2013) projections back from the SNcVTA to the DRN (Kalen et al., 1988;Peyron et al., 1996;Kitahama et al., 2000) might support this effect.In addition, dopamine agonist increases excitation on DRN serotonergic neurons (Haj-Dahmane, 2001;Martın-Ruiz et al., 2001).However, there is also literature supporting that there is no change in the electrical activity of serotonergic-like neurons following L-Dopa injection (Miguelez et al., 2011(Miguelez et al., , 2016)), that lesion in VTA could lead to a reduction of serotonin in DRN, and that supports a DRN activity reduction after a dopaminergic lesion (Furlanetti et al., 2016).The different results showed in DRN firing activity showed in this literature might depend on several factors such as the type of lesion (unilateral or bilateral), the site of lesion, and time of recording after lesion.Moreover, we excluded from our model several brain regions that might affect DRN regulation, for example, the medial prefrontal cortex and subthalamic nucleus playing a critical role in the inhibition of DRN (Celada et al., 2001;Temel et al., 2007;Hartung et al., 2011).It could be possible that during the initial stages of the disease, the DRN activity increases, while it decreases with the disease progression.The inversed ushape followed by the DRN activity with the PD progression could

FIGURE
The only feature that di erentiates, albeit with low significance, individuals for which it was possible to find a treatment stimulating either LC or DRN is the sensitivity of SNcVTA to noradrenaline from LC, for both the linear (A) and non-linear (B) terms.** p ≤ . .

FIGURE
Distribution comparison of model parameters.In (A), a distribution comparison of model parameters in LDA subjects relative to treated subjects is presented.In (B), the parameter distribution of treated subjects is depicted, distinguishing between those successfully cured (green) and those who did not attain the desired fitness (red).Notably, among the non-responsive subjects, one exhibited reduced sensitivity of SNcVTA to noradrenaline.be similar to that followed by LC activity during PD and AD progression (Caligiore et al., 2020(Caligiore et al., , 2022;;O'Callaghan et al., 2021;Ye et al., 2023).Further research could confirm or refute this hypothesis.
The simulations run with the model show a reduction in striatal D1 activity and an increase in striatal D2 activity (see Figures 3,  4).This result agrees with the literature about hypokinetic Parkinsonian syndrome due to the activity dysregulation of the two populations of medium spiny neurons (MSNs).Dopamine D1 receptor-expressing MSNs (direct) become hypoactive, whereas dopamine D2 receptor-expressing MSNs (indirect) become hyperactive.The global effect is the excitation of the direct pathway and the inhibition of the indirect pathway, both conditions necessary for correctly managing motor output (Albin et al., 1989;Mallet et al., 2006;Parker et al., 2018;Caligiore et al., 2019).
The simulations also show a reduction of LC activity after dopaminergic depletion (Figures 3, 4).The DA excitatory effect on LC reproduced in our model could explain this result.However, the literature supports both decrease and increase in LC, depending on the PD animal models used for experiments (Ranjbar-Slamloo and Fazlali, 2020).Thus, future in vivo studies should confirm or refute the model result.The serotonin lesions allow reproduction of the GP firing activity decrease reported in the target experiment (Delaville et al., 2012).Moreover, the lesioned model shows an increase in SNcVTA activity (Figures 3, 5).The 5-HT inhibitory effect on dopaminergic neurons reproduced in our model could explain this result.However, future studies are required to confirm these data.Indeed, 5-HT could modulate SNpc and VTA DA neurons oppositely (Gervais and Rouillard, 2000).

. Exploring potential monoamine-based therapies
Levodopa is the most common medication used in PD.However, this drug has a wide range of adverse effects, most notably motor fluctuations and dyskinesias (for a rev see Lang and Obeso, 2004).The discovery of alternative treatment that not only targets dopaminergic system but also noradrenergic or serotonergic system is a big challenge (Politis and Niccolini, 2015;Wilson et al., 2019;Caligiore et al., 2021Caligiore et al., , 2022)).Figures 9-12 collectively indicate that the restoration of activity in the LC is sufficient to restore the regions under examination, while the DRN does not yield the same effect.This observation is consistent with current theories that underscore the central role of the LC in the progression of PD, particularly in relation to the early-stage emergence of non-motor symptoms (Bjerkén et al., 2019;Butkovich et al., 2020).Moreover, there is literature supporting the notion that the overexpression of specific transcription factors directly within the LC can effectively restore noradrenergic function and facilitate the recovery of the dopaminergic system (Cui et al., 2021).Notably, in the context of human studies, mounting evidence underscores the fact that LC degeneration can manifest at an earlier stage and with greater severity compared to the SNcVTA (German et al., 1992;Del Tredici et al., 2002;Rommelfanger and Weinshenker, 2007;Vermeiren and Deyn, 2017).Taken together, these findings strongly suggest that the LC may play a pivotal role in the emergence and progression of PD and present an intriguing avenue for therapeutic interventions targeting both the dopaminergic and noradrenergic systems.
The feasibility of this discovery is not straightforward.One plausible explanation for the lack of productivity of novel drugs targeting monoaminergic function could be the complexity of the brain neurotransmitter systems.While these systems play crucial roles in regulating mood and cognition, their functions are highly interconnected and often involve feedback loops and compensatory mechanisms.It is possible that the drugs developed to target these systems may inadvertently disrupt the delicate balance of neurotransmitter activity, leading to unintended side effects or diminishing efficacy over time.In addition, the existing drugs may primarily target specific receptor subtypes within the noradrenergic systems, leaving other potential therapeutic targets unexplored.This limited scope could hinder the development of more effective treatments (Tan et al., 2022;Pardo-Moreno et al., 2023).Finally, while the model suggests promising alternative treatments targeting the LC and DRN, it is crucial to consider the potential side effects associated with these approaches.Given the critical roles these regions play in regulating mood, arousal, and autonomic functions, interventions could lead to increased anxiety, sleep disturbances, and cardiovascular issues.These potential side effects highlight the need for comprehensive experimental validation and clinical studies to assess the safety and efficacy of these alternative treatments.Addressing these concerns will be essential for translating our model predictions into practical therapeutic strategies.

Conclusion
This article presents a novel system-level computational model that, for the first time, explores the roles of three important monoamines (dopamine, serotonin, and noradrenaline) in a PD animal model.The model assumes a direct relationship between reduced activity in specific brain areas (SNcVTA, DRN, and LC) and decreased neurotransmitter levels (DA, 5-HT, and NE) in target brain regions.In addition, it simulates VTA and SNc together and excludes some potentially relevant areas, such as the prefrontal and motor cortices, thalamus, and substantia nigra pars reticulata.Despite these simplifications, the model successfully reproduces some data on combined monoamine depletion, generates predictions regarding changes in other unexplored brain regions, suggesting avenues for further investigation, and highlights the potential efficacy of alternative treatments targeting the locus coeruleus and dorsal raphe nucleus.This suggests that the model captures crucial aspects of the interactions between the considered areas.Furthermore, it may also imply that the aspects not yet modeled might only play secondary roles in the system behaviors.
We do not subscribe to this belief.Our plan is to expand the model to incorporate more areas involved in PD (Dirkx et al., 2017;Helmich et al., 2018), to more accurately delineate macro-areas like SNc/VTA (Ledonne et al., 2023), and to include additional interaction paths.The upgraded versions of the model may be wellsuited for delving into more mechanisms, such as those elucidating the effects of dopamine receptor activation (Mailman et al., 2021;Lewis et al., 2023) or exploring monoamine interactions within various brain regions (Oh et al., 2022).At this stage, the model represents an initial effort to employ a system-level computational approach to address the intricacies of PD as a systemic disease, with a focus on the interplay of monoamines.It serves as an initial attempt to demonstrate the feasibility of this approach, employing a rigorous mathematical framework specifically, the study of stability, a seldom-used approach in computational neuroscience.The results obtained through the model are preliminary, and their effectiveness requires further validation.
In this respect, a series of experimental approaches could validate the model predictions.In vivo, neurochemical assessments using microdialysis in PD animal models could measure extracellular levels of DA, NE, and 5-HT in the basal ganglia.Pharmacological manipulations could selectively alter NE and 5-HT levels, observing the effects on PD symptoms and neurochemical dynamics to validate the model predictions about these systems.Neuroimaging studies using functional magnetic resonance imaging (fMRI) and positron emission tomography (PET) can assess functional connectivity and receptor density changes in brain regions highlighted by the model.This analysis could help validate the predicted alterations in neuromodulator levels and activity patterns.In addition, analyzing functional connectivity using high-density electroencephalography (HD-EEG) could validate the model on a macroscale by measuring changes in band-specific cortico-cortical and subcortico-cortical connectivity before and after manipulations of the monoaminergic system (Conti et al., 2023).Moreover, electrophysiological recordings of neuronal activity in relevant brain regions before and after such manipulations could validate the model's predictions regarding neural circuitry alterations.
The model proposed here could enhance our comprehension of interactions between brain regions in both normal and pathological conditions, potentially aiding in the restoration of damaged brain regions to reestablish balance.The model could be adapted to study other neurodegenerative and neuropsychiatric disorders involving similar monoaminergic systems.For instance, there is increasing evidence of the critical role of monoamines in AD (Caligiore et al., 2022) and multiple sclerosis (Carandini et al., 2021).Adapting the model to these disorders would involve modifying the parameters to reflect the specific pathophysiological mechanisms and neurochemical dynamics unique to each condition.Future research could address these points.

FIGURE
FIGUREConceptual model schema.The average activation frequencies of six brain areas are modeled (rounded rectangles); some interactions are modulated by monoamines (circles).Arrows represent positive (excitatory) e ects, while circles represent negative (inhibitory) e ects.Noradrenaline has a non-linear (both excitatory and inhibitory) e ect on SNcVTA which is indicated by a bar.Each area has a corresponding stimulus (ovals) which represents self-activation as well all any other stimulus the area might receive from the rest of the brain which is not modeled.Lastly, hexagons serve as indicators for identifying the areas influenced by the administration of specific drugs.

FIGURE
FIGUREGenerated target value distribution for all cases: in black, the SHAM case.The LDA dopaminergic lesion target (red) di ers from SHAM only for SNcVTA and LC.The serotonergic L HT lesion target (yellow) di ers from SHAM only in for GP and DRN.Finally, the noradrenergic LNE lesion target (blue) di ers from SHAM only for LC.Values for each individual are generated according to Section .: Each area follows a normal distribution around a center value with a maximum spread of ± % ( σ = .µ). Lesioned activation values of an area, when defined, are scaled according to Table.

FIGURE
FIGUREDistribution of simulated equilibrium points for the GP (A), StrD (C), and StrD (E) regions across the entire synthetic population.On the right, the simulated equilibrium points for GP (B), StrD (D), and StrD (F) are depicted within separate bins, each containing samples.

FIGURE
FIGUREDistribution of simulated equilibrium points for the SNcVTA (A), DRN (C), and LC (E) regions across the entire synthetic population.On the right, the same simulated equilibrium points for SNcVTA (B), DRN (D), and LC (F) are displayed within separate bins, each containing samples.

FIGURE
FIGUREComparative sensitivity of each area to individual parameters among healthy subjects.A higher value denotes a more pronounced impact (in absolute magnitude) of a parameter on the activation frequency of an area.Values are normalized to unity, highlighting the relative magnitude of e ects obtained from each parameter.

FIGUREE
FIGUREE ects of the dopaminergic (Parkinsonian) lesion to SNcVTA and LC induced by LDA.SNcVTA activation, and consequently the production of dopamine, is drastically lowered compared to the reference (gray) levels.LC activity is also lowered to % of its SHAM value.GP values remain unaltered as constrained by the fitness function; SNcVTA and LC are also subject to a softer constraint (see Appendix .). ).

FIGURE
FIGURELesion and treated values for GP (A), StrD (B), StrD (C), SNcVTA (D), DRN (E), and LC (F).A statistically significant boost of LC average activity (and hence of noradrenaline levels) can restore the activity (and hence monoamine production levels) of all the areas that were significatively impacted by LDA to SHAM levels.** p ≤ .and **** p ≤ ..

FIGURE
FIGURE

FIGURE
FIGUREComparison of parameter distribution among model subjects: SHAM subjects represented in black, while LDA subjects under two conditions, either LDA (A) or treated (B), depicted in orange.As outlined in Section ., only four parameters, impacting the SNcVTA equation, undergo alteration in the LDA scenario compared to SHAM.In addition, the treatment induces changes in external stimulation to LC and DRN.
TABLE Average activity target of the simulated healthy subjects, and time constants used in simulated healthy subjects.The arbitrary tolerance of ±0.3 has been added to account for the variability observed in the data found in the literature.

Table 2
TABLE Reference data for area interaction and lesion e ects.
summarizes the reference changes in GP activity.