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, 1978; Kelland et al., 1990, 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 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:

where the abbreviated notation ẋ stands for dxdt and:

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:(ℝ × ℝ^s) → ℝ^s is defined by the corresponding equation in Equations 1–6, and where we assume the initial state y(t₀) = y₀ to be known.

Starting from Equation 8, Equations 1–6, 8 can be represented in matrix form:

where a_ij = 0 if the corresponding αij is not defined and likewise b_i = 0 if αiext is not defined, and also cij=βij with c_ij = 0 if the corresponding βij coefficient is not defined; “°” indicates the element-wise vector product (or Hadamard).

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, αSNcVTADRN, αSNcVTALC, βSNcVTALC, αSNcVTAext. Likewise, when the subject also receives a serotonergic lesion, there will be a third set of parameters A_LDA+L5HT, C_LDA+L5HT, b_LDA+L5HT which again differ from A_LDA, C_LDA, b_LDA only by the parameters corresponding to the equation for DRN, and so on.

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.

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.

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.

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 Sikind, such as SiSHAM and SiLDA, 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₀, T], with the starting vector y₀, 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 Tikind the corresponding reference solutions that will be used to evaluate the fitness of the model:

where J = [t₀, ..., t_n] is a vector of times. The fitness of a model is finally obtained as a combination of many fitness figures f_i∈[0, 1], 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.

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:

where:

Table 2 summarizes the reference changes in GP activity.

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.

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σ=12μ. 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.

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.

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²). 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² =1/second² =Hz/s. Finally, all time constants τ are naturally expressed in seconds.

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.

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:

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.

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.

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:

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.

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