A Data-Driven Biophysical Computational Model of Parkinson's Disease based on Marmoset Monkeys

07/27/2021 ∙ by Caetano M. Ranieri, et al. ∙ Universidade de São Paulo University of Campinas Heriot-Watt University 0

In this work we propose a new biophysical computational model of brain regions relevant to Parkinson's Disease based on local field potential data collected from the brain of marmoset monkeys. Parkinson's disease is a neurodegenerative disorder, linked to the death of dopaminergic neurons at the substantia nigra pars compacta, which affects the normal dynamics of the basal ganglia-thalamus-cortex neuronal circuit of the brain. Although there are multiple mechanisms underlying the disease, a complete description of those mechanisms and molecular pathogenesis are still missing, and there is still no cure. To address this gap, computational models that resemble neurobiological aspects found in animal models have been proposed. In our model, we performed a data-driven approach in which a set of biologically constrained parameters is optimised using differential evolution. Evolved models successfully resembled single-neuron mean firing rates and spectral signatures of local field potentials from healthy and parkinsonian marmoset brain data. As far as we are concerned, this is the first computational model of Parkinson's Disease based on simultaneous electrophysiological recordings from seven brain regions of Marmoset monkeys. Results show that the proposed model could facilitate the investigation of the mechanisms of PD and support the development of techniques that can indicate new therapies. It could also be applied to other computational neuroscience problems in which biological data could be used to fit multi-scale models of brain circuits.



There are no comments yet.


page 12

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Parkinson’s disease (PD) affects more than 3% of people over 65 years old, with figures set to double in the next 15 years [Poewe2017]. It is a neurodegenerative disease, whose symptoms include cognitive and motor deficits. In late stages, it can possibly also lead to depression and dementia [Shulman2011]. There is still no cure, and current therapies are only able to provide symptomatic relief.

PD is characterised by a dopaminergic neuronal loss within the substantia nigra pars compacta (SNc), which leads to a dysfunction of the basal ganglia-thalamus-cortex (BG-T-C) circuit. The BG-T-C circuit is a neuronal network with parallel loops that are involved in motor control, cognition, and processing of rewards and emotions [Obeso2009TheObservations, Redgrave2010]. There are also links between the degeneration of dopamine neurons within those brain regions and changes on electrophysiological behaviour [galvan2008pathophysiology].

A commonly used model to explain how PD affects the neural connections within this circuit, also known as the motor loop, is the so-called classic model, illustrated in Figure 1

a. It consists of projections from primary motor (M1) and somatosensory cortical areas to BG input structures, specifically the putamen (PUT) and the subthalamic nucleus (STN). In PUT, the cortical projections establish excitatory glutamatergic synapses with medium spiny neurons (MSNs).

The MSNs establish two distinct pathways to the BG output nuclei (globus pallidus pars interna – GPi and substantia nigra pars reticulata – SNr). The MSNs from the direct pathway (dMSN) directly project to the GPi/SNr, while the MSNs from the indirect pathway (iMSN) project to the globus pallidus pars externa (GPe), which in turn send projections to the GPi/SNr directly or indirectly via the STN (for reviews, see Obeso et al. [Obeso2009TheObservations], Lanciego et al. [Lanciego2012], and McGregor and Nelson [McGregor2019CircuitDisease]).

The cortical projection to the STN establish a third pathway, often called the hyperdirect pathway [Nambu2002]. Activation of the direct pathway facilitates movement by inhibiting the activity of GPi/SNr, thus reducing the inhibition of the ventral anterior nucleus (VA) and the ventral lateral nucleus (VL) and increasing the excitatory thalamic input to the motor cortex. Activation of the indirect and hyperdirect pathways, on the other hand, inhibit movement by increasing the inhibitory activity of the GPi/SNr over the VA/VL, hence decreasing the excitatory thalamic input to the motor cortex.

The activity of the motor loop is modulated by dopaminergic projections from SNc to PUT. The main effect of dopamine (DA) release in PUT is movement facilitation, since DA increases the excitability of the dMSNs and decreases the excitability of the iMSNs.

In PD, the depletion of striatal DA leads to an enhanced activation of the indirect pathway and a decreased activation of the direct pathway, resulting in the characteristic motor symptoms of this neural disorder [Surmeier2014]. In addition to changes in firing rates, the functional imbalance within the motor loop in PD also disrupts the firing patterns within each nucleus and amongst the structures of the BG-T-C circuit, increasing neuronal synchronisation, neuronal bursting, and enhancing the oscillatory activity at the beta frequency band [Galvan2015].

Figure 1: Models of the basal ganglia-thalamus-cortex (BG-T-C) circuit, central to the underlying mechanisms of Parkinson’s Disease (PD), including excitatory (blue) and inhibitory (red) connections between the regions involved. (a) Classical model of BG-T-C circuit. The motor loop in the mammalian brain is formed by the motor cortex (M1), the thalamus (TH) - composed of structures such as the ventral anterior nucleus (VA), the ventral lateral nucleus (VL), and the ventral posterolateral nucleus (VPL) -, and the basal ganglia (BG), the latter composed of a subset of structures: the striatum, which itself includes the putamen (PUT) and the caudate nucleus, the globus pallidus, divided into pars interna (GPi) and pars externa (GPe), the subthalamic nucleus (STN), and the substantia nigra, divided into pars compacta (SNc) and pars reticulata (SNr). PD is caused by the loss of dopaminergic neurons in the SNc, which weakens the connections represented by dashed lines and leads to malfunctioning of both direct and indirect pathways. (b) BG-T-C network used in this work, based on [Kumaravelu2016ADisease]. The cortex is represented by regular spiking (CtxRS) excitatory neurons and fast spiking (CtxFSI) inhibitory interneurons. The direct and indirect pathways in the striatum were modelled separately, representing the medium spiny neurons (MSNs) modulation by D1 and D2 dopamine receptors, respectively.

Brain regions linked to PD present complex interactions, with mutual excitatory and inhibitory feedback loops, which limit a comprehensive understanding of the physiopathology of the disease. Studies aiming at investigating the mechanisms underlying PD often use animal models. In classic animal models of PD, symptoms are elicited by delivering neurotoxins that damage the SNc dopaminergic neurons, such as 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) and 6-hydroxidopamine (6-OHDA), or chemicals that transiently inhibit dopamine production, such as alpha-methyl-p-tyrosine (AMPT) [Koprich2017AnimalDevelopment]. Also, antipsychotics like haloperidol have side effects that may promote dystonia and parkinsonism [KHARKWAL2016].

Bilateral 6-OHDA lesions in the marmoset medial forebrain bundle induce several PD motor symptoms, including impairments in fine motor skills, limb rigidity, bradykinesia, hypokinesia, and gait impairments. Alpha-methyl-p-tyrosine (AMPT) administration to 6-OHDA lesioned marmosets can transiently increase the severity of these symptoms. However, like MPTP macaques, these animals do not exhibit resting tremor. Santana et al. [Santana2015] provides an extensive characterisation of these symptoms, that were quantified through manual scoring (adapted version of the unified PD rating scale for marmosets), automated assessments of spontaneous motor activity in their home cages (using actimeters), and automated motion tracking while the animals explored two experimental apparatuses.

To date, no animal model of PD fully reproduces human features of the disease. In addition, due to experimental limitations, animal data often include only a limited set of PD-related brain regions, with subjects engaged in different behavioural settings. In this context, computational models, with biologically informed constraints that can be selectively altered, are a promising, complementary approach to advance our knowledge about PD beyond that obtained from anatomical and physiological studies [Humphries2018, navarro2020dynamical]. Some PD-related anomalies observed in animal models, and efforts to reproduce those in computational models, are presented by Rubin et al. [rubin2012basal].

Computational models are established tools to facilitate understanding of neural disorders [Schroll2013, Sanger2018, Pena2020] and, in the context of PD, accommodate several levels of description and range from focusing on disease mechanisms to understanding anomalous neuronal synchronisation [Humphries2018].

For instance, Pavlides et al. [Pavlides2015] conducted a detailed study to help unveiling the mechanisms underlying beta-band oscillations in PD and compared computational model predictions with experimental data. Muddapu et al. [Muddapu2019] studied loss of dopaminergic cells in the SNc due to neural dynamics between SNc and STN, shedding light on the relevance of ongoing neural activity and neural loss. Gurney et al. [Gurney2004] described mounting evidence relating the BG-T-C network and action selection mechanisms; actually, computational models showed a close relationship between action selection and BG-T-C oscillatory activity [Humphries2006, Holgado2010, Merrison2017].

Moren et al. [moren2019dynamics] proposed a model of the spiking neurons within the BG-T-C circuit, in order to observe the asynchronous firing rates around the 15 Hz beta-range oscillations, as well as on lower frequency bands. Terman et al. [Terman2002] developed a conductance-based computational network model which shed light on the mechanisms underlying the neural dynamics of STN and GPe, a model which was further developed by Rubin et al. [Rubin2004] to investigate the effects of deep brain stimulation (DBS) to eliminate anomalous synchronisation within the BG-T-C network in PD condition. In fact, one of the key areas in which computational models serve as an invaluable tool for developing novel therapies is that related to predicting the effects of DBS [Humphries2012, Lu2020].

Finally, based on a collection of previously published studies, Kumaravelu et al. [Kumaravelu2016ADisease] developed a computational model of the BG-T-C network tuned for the 6-OHDA rat model of PD (Fig. 1b). Compared to other computational models [Humphries2018], it was the first to specifically consider 6-OHDA and a single species.

Most computational models related to BG-T-C dynamics rely on rodent data [Humphries2006, Lindahl2016, Koelman2019], with only a handful focusing on primate data [Lienard2014, Topalidou2018]. The research by Shouno et al. [shouno2017computational], for instance, provided a spiking neuron model of the recurrent STN-GPe circuit for studying dysfunctions in oscillations within the 8-15 Hz frequency band for PD primate models.

All mammals have a similar set of BG structures that are similarly connected with thalamic and cortical structures. Nevertheless, recent studies suggest subtle differences between species [Lienard2014, Hardman2002, VanAlbada2009, Beul2017], also in the neuropathophysiology of PD [Koprich2017AnimalDevelopment, Dawson2018], with primates (including marmosets) being more similar to humans than rodents. For example, there are differences in the distribution of dopaminergic neurons in substantia nigra of rats and primates, and the subthalamic nucleus and internal globus pallidus of rats have less neurons containing parvalbumin than primates [Hardman2002]. Thus, a primate computational model of PD is of paramount relevance.

In this work, we developed a new computational model of PD based on published data from the BG-T-C brain circuit of marmoset monkeys  [Cyranoski2009]. We built upon the neuronal computational model of rat models of PD [Kumaravelu2016ADisease], and adjusted its parameters to match the electrophysiology data from 6-OHDA+AMPT marmoset model of PD [Santana2014SpinalDisease, Santana2015].

It is important to highlight that, in our work, we are using the LFP signal data to tune and validate our model, not spikes or other biosignals, thus the whole optimisation framework relates to LFP-based metrics. We are aware that there are several simplifications in our computational model, nevertheless results were shown that LFP power spectral densities at frequencies of interest, firing frequency dynamics, and spike coherence resemble those from healthy and PD marmosets.

The main contributions of this paper are: (i) the first computational model of PD validated on simultaneous, multi-site electrophysiological recordings (e.g., LFP recordings) from a marmoset monkey model of the disease, and (ii) an optimisation framework that can easily include novel biophysical parameters as soon as they become available.

This paper is organised as follows. In Section 2, the building blocks of the computational model are depicted, as well as the free parameters that were optimised, the algorithm to update those parameters, the experimental setup, and the evaluation protocol. In Section 3

, the results are presented regarding the optimisation process, the parameters learnt by the machine learning algorithms, and the metrics observed on the simulations of the computational models provided, considering spectral densities from simulated LFP, dynamics of the firing rates from simulated neurons, and coherence analyses. In Section 

4, a discussion is presented in order to contextualise our results and compare them with the expectations from the data from animal models, and knowledge from the literature. In Section 5, a conclusion is presented with a brief summary of what was presented.

2 Methods

To provide a computational model of the BG-T-C circuit for PD-related features in primates, we began by re-writing the code by Kumaravelu et al. [Kumaravelu2016ADisease], originally implemented in Matlab. We have ported the original code to the Python programming language, with the NetPyNE framework and the libraries from the NEURON simulator [Hines2009, Dura2019]. Then, we performed a series of adaptations and employed a data-driven approach to calibrate a set of parameters, in order to derive a model that resembles local field potentials from marmoset data [Santana2014SpinalDisease, Santana2015]. More specifically, we employed an optimisation technique called differential evolution (DE), an algorithm based on evolutionary computation [ashlock2006evolutionary]. This approach consists of optimising a predefined set of parameters (i.e., genotype) by gradually adapting them through successive steps (i.e., generations), providing variability and selection of the best solutions (i.e., individuals) through mechanisms analogous to biological evolution.

In the model by Kumaravelu et al. [Kumaravelu2016ADisease], no noise was introduced in the simulations. Neuronal connectivity and membrane initial conditions can be stochastic, and neuronal models include synaptic transmission delay. The dataset employed in our work was suitable to calibrate such a model, since it was collected from marmoset monkeys that were not engaged in any particular task, that is, they were moving freely, without any time-marked events such as sensory or artificial stimuli.

After having calibrated our marmoset model, different analyses were performed in order to enhance and validate it. The dataset used as ground-truth for adjusting the parameters of the computational model is not publicly available due to legal restriction, but it is available from the corresponding author on reasonable request. The next subsections will provide a detailed description of the methods employed. The code to reproduce the results from this paper, including the machine learning framework and the analyses of the results, is publicly available at https://github.com/cmranieri/MarmosetModel.

2.1 Computational Model

The computational model was based on Kumaravelu et al. [Kumaravelu2016ADisease]. Their model was build to reproduce the neurophysiological behaviour from rats based on data from healthy and 6-OHDA-lesioned individuals. As an initial step, we did an alternative implementation for their model within the NetPyNE framework, and we validated this implementation by comparing its outputs with those reported in [Kumaravelu2016ADisease].

Briefly, eight brain structures were modelled and connected based on a simplified version of the classic model (Figure 1b). In particular, the direct and indirect pathway in the striatum were modelled separately representing the MSN modulation by D1 and D2 dopamine receptors, respectively [McGregor2019CircuitDisease]. The cortex is represented by regular spiking (RS) excitatory neurons and fast spiking (FSI) inhibitory interneurons. Neurons from all but cortical regions were modelled using a biophysically based Hodgkin–Huxley [hodgkin1952currents] single-compartment model, whereas cortical neurons were constructed based on the computationally efficient Izhikevich’s model [Izhikevich2003]. The reasoning for different neuronal models lies on the fact that PD effects are captured by altering specific conductances in selected structures (see below), thus a conductance-based model is more suitable at these locations. Finally, a bias current was added in the TH, GPe, and GPi, accounting for the inputs not explicitly modelled. Remarkably, even though no oscillatory inputs are present in the model, synaptic delays and network interactions by means of recurrent connections promote sustained firing rate oscillations. For a detailed description of connectivity schemes and other implementation details, the reader is referred to Kumaravelu et al. [Kumaravelu2016ADisease].

The computational model described above can shift from the simulations of the healthy to the PD conditions by altering three conductances [Kumaravelu2016ADisease]: decreasing the maximal M-type potassium conductance in direct and indirect MSN neurons (MSN firing disfunction) from 2.6 to 1.5 ; decreasing the maximal corticostriatal synaptic conductance (reduced sensitivity of direct MSN to cortical inputs) from 0.07 to 0.026 ; and increasing the maximal GPe axonal collaterals synaptic conductance from 0.125 to 0.5 (increase of GPe neuronal firing). This is implemented in the model with a control flag.

One major addition to the model developed here is the simulation of local field potentials (LFP). These measurements are related to the extracellular activity produced by action potentials of the neurons within a brain region [gold2006origin]. A discussion on the behaviour of LFP signals within the basal ganglia and its consequences to humans, especially regarding conditions such as PD, was presented by Brown and Williams [brown2005basal]. The NetPyNE function for LFP calculation is based on the work of Parasuram et al. [Parasuram2016]. The LFP amplitude in each simulated electrode is obtained by summing the extracellular potential contributed by each neuronal segment, calculated using the line source approximation and assuming an Ohmic medium with conductivity 0.3 mS/mm. Thus, the electrical activity of neurons from each brain region contributes to the peaks and valleys recorded at each electrode (subject to extracellular medium attenuation).

In our work, first, each simulated brain region is assigned to a spatial 3D coordinate that matches that used in the stereotaxic surgery where electrodes were placed in the real marmoset monkeys [Paxinos2012, Santana2014SpinalDisease]. Then, a simulated electrode is placed at the centre of each region. In our model, each neuron is represented as a single cylindrical compartment with a membrane area of 100

. For each electrode, NetPyNE estimates the simulated LFP by summing the extracellular potential contributed by each neuronal segment (based on the transmembrane current generated from the single cylindrical source neuron), calculated using the “line source approximation” method and assuming an Ohmic extracellular medium with conductivity


2.2 Dataset and preprocessing procedures

The dataset we used in the present work is based on a previous study by Santana et al. [Santana2014SpinalDisease]. Our dataset includes data from three adult males and one adult female common marmosets (i.e., Callithrix jacchus). Data from two males were part of the aforementioned study; data from one male and one female are novel and followed exactly the same experimental procedures. A short summary is presented in the next subsection, followed by the preprocessing steps.

2.2.1 Dataset

The animals, weighing 300–550 g, were housed in a vivarium with natural light cycle (12/12 hr) and outdoor temperature. All animal procedures followed approved ethics committee protocols (CEUA-AASDAP 08/2011, 11/2011 and 03/2015) strictly in accordance with the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals. PD symptoms were elicited in all three male animals with injections of 6-OHDA toxin in the medial forebrain bundle under deep anaesthesia [Santana2015, Santana2014SpinalDisease]. Prior to neural recordings, animals that received 6-OHDA were subjected to acute pharmacological inhibition of dopamine synthesis (subcutaneous injections of AMPT mg/kg) to further exacerbate PD motor symptoms, mimicking a more severe stage of the disease. Although 6-OHDA lesions impact on both behavioural and electrophysiological features in all animals [Santana2015, Santana2014SpinalDisease], there are individual differences at earlier stages of dopaminergic depletion that could hinder our model development considering the relatively low number of subjects.

Both healthy and PD animals were implanted each with two custom-made microelectrode arrays composed of 32 microwires (one array in each hemisphere). The wires were 50 in diameter and were organised in bundles aimed to reach distinct areas of the BG-T-C system. Before the surgery, the animals were sedated with ketamine (10-20 mg/kg i.m.) and atropine (0.05 mg/kg i.m.), followed by deep anesthesia with isoflurane 1-5% in oxygen at 1-1.5 L/min. The arrays were then implanted using a stereotaxic manipulator to position electrodes at the targeted BG-T-C coordinates, which were determined using Stephan et al. [Stephan1980] and Paxinos et al. [Paxinos2012] stereotaxic atlas. The microelectrode array and the implant procedures were thoroughly described in Budoff et al. [Budoff2019].

Once the animals recovered from the surgery, recording sessions were performed in fully awaken animals behaving freely. LFPs were sampled at 1,000 Hz and recorded using a 64 multi-channel recording system (Plexon). The position of the recording microelectrodes were verified postmortem through either tyrosine hydroxylase (TH) staining or Nissl staining. Similarly, the extent of dopaminergic lesions were verified through the quantification of striatal fiber density and dopaminergic midbrain cells in TH-stained sections. Further experimental details are described in Santana et al. [Santana2014SpinalDisease].

2.2.2 Preprocessing

For our study, in total, 14 and 16 recording sessions were taken for the healthy and PD conditions, respectively, considering the brain hemispheres independent from each other. For the PD condition, we recorded from M1, PUT, GPe, GPi, ventral lateral (VL) and ventral posterolateral (VPL) thalamic nuclei, and STN, whereas for the healthy animal regions M1, PUT, GPe, and GPi were recorded. The raw data was organised so that, for each recording session, a data structure with was provided, where is the number of electrodes recorded and is the number of samples of the recording session (variable but typically lasting for several minutes).

In Figure 2

, are illustratrated the preprocessing steps adopted after data acquisition. For each channel, the pipeline began with a zero-lag low-pass filter (cutoff frequency of 250 Hz) and a high-pass filter (cutoff frequency of 0.50 Hz), to eliminate frequencies that are outside the LFP scope and may relate to electrical or mechanical interference. Then, we minimised power grid interference (hum) with a notch filter centred at 60 Hz and its harmonics (120 Hz and 180 Hz). Each resulting signal was then scaled according to a z-score normalisation, to account for the possible differences in signal amplitude due to different electrode impedance.

Figure 2: Data acquisition and preprocessing steps implemented in out method. Depending on the monkey condition, different regions of the brain were recorded. The input data was composed of a whole recording session, with variable lengths and numbers of channels (i.e., electrodes) per region. After preprocessing, the data was transformed into 2-seconds-long segments with seven channels, each related to one of the regions analysed.

In the next step, we computed the cross-correlation matrix according to Equation 1, where is the covariance matrix of the filtered and z-scored signals from electrodes and , which are located exclusively within a brain region. Channels within each region with mean correlation coefficient below the threshold of 0.70 were discarded. This procedure was employed because electrodes in each recorded region are placed very close to each other (see electrode and surgical procedures above), thus we expect LFP signals to be highly correlated (if they are not, it may relate to a noisy electrode signal) [Buzsaki2012TheSpikes.].


All remaining LFP channels within a brain region were averaged, which provided one data matrix for each recording session with dimensions , where is the number of brain regions recorded. These average LFP values were computed based solely on channels within each region. Next, we segmented each time-series in 2-second segments, which was the same length of the computational model simulations (Section 2.3

will bring the details). Considering the data sampling rate (1,000 Hz) and frequencies of interest (up to 50 Hz), 2-second segments provide enough data for our analyses. Prohibitively noisy segments were discarded using two criteria: first, segments with abnormal amplitudes, detected using an upper threshold of 0.20 for the absolute value of the mean of the signal over time; second, segments with limited (abnormal) oscillatory patterns, detected using a lower threshold of 0.10 for the amplitude standard deviation and a minimum threshold of 10 amplitude peaks. For each recording session, our preprocessed dataset had a final shape of

, where is the resulting number of segments.

In the dataset adopted for this work, whether animals were still or moving could have a profound effect on brain oscillatory activity and synchronisation metrics, because all animals were behaving freely and were not engaged in any particular behavioural task during the recording sessions. In fact, especially in motor and pre-motor regions, modulations in neural oscillatory dynamics linked to motor activity are well characterised (see Armstrong et al. [Armstrong2018] for a review), and recent studies show that even breathing can modulate neural oscillations [Tort2020]. However, we understand that action initiation, movement, or breathing have low influence on averaged LFP amplitude values computed, given that the 2-second window segments were randomly selected without time alignment to any specific movement or action.

2.3 Evolutionary Algorithm

Evolutionary algorithms are optimisation techniques in which a set of parameters, called genotypes, are gradually combined and changed according to mechanisms analogous to those of biological evolution, in order to maximise a fitness function dependent of those parameters [ashlock2006evolutionary]. Differential evolution (DE) [RAINERSTORN1997] was employed to fit the computational model parameters so that it matches the LFP beta-band power spectrum observed in the marmoset data.

The overall structure of the model was preserved from Kumaravelu et al. [Kumaravelu2016ADisease], while a set of conductances, background currents and synaptic modulations, as well as the numbers of neurons in each region of the BG-T-C circuit, were calibrated through the evolutionary algorithm. The connectivity, the delays, the synaptic mechanisms, the remaining conductances, and all other parameters were kept as in the original model (see Section 3 from the Supplementary Material).

More specifically, fourteen parameters compose the set of parameters to be optimised (i.e., the genotype). Parameter () relates to cerebellar input currents to the thalamus, which are linked to sensorimotor inputs [Manto2013]. Parameters () and () relate to currents at GPe and GPi, respectively, from all sources that were not explicitly modelled. The next two parameters, (nS/cm2) and (nS/cm2), refer to the maximum slow potassium conductance yielding afterhyperpolarization (AHP) at the STN and the calcium-activated potassium conductance at GPe and GPi, respectively. The sixth parameter, (nS/cm2), modifies the synaptic conductance from cortex (CTX) to striatum (STR). Finally, parameters seven to 14 map to the number of neurons in each modelled region. All of the aforementioned parameters were chosen because they have a direct influence on the firing rates of neurons within each region, which in turn affect the LFP [Parasuram2016]. Also, comparing marmoset with rodent literature, there is very limited quantitative work on the anatomical and neurophysiological parameters of the BG-T-C neuronal network.

ID Parameter Range Description
1 Background currents at TH ()
2 Background currents at GPe ()
3 Background currents at GPi ()
4 –dependent AHP conductance
at STN ()
5 –dependent AHP conductance
at GPe and GPi ()
6 Synaptic modulation from cortex to
striatum ()
7 Number of GPe neurons
8 Number of GPi neurons
9 Number of TH neurons
10 Number of StrD1 neurons
11 Number of StrD2 neurons
12 Number of CTX_RS neurons
13 Number of CTX_FSI neurons
14 Number of STN neurons
Table 1: Free parameters of the computational model, optimised by DE to fit the marmoset data.

In the DE, each individual from the population was a model that consisted of an adaptation of the model of Kumaravelu et al. [Kumaravelu2016ADisease], in which the parameters of Table 1 were set to the values defined by genotype . Each model was simulated for milliseconds, and the spike trains from each neuron and LFPs from each virtual electrode were recorded. The LFP recordings were applied to calculate the fitness function as follows.

Given a categorical set containing brain regions, the mean power spectral density (PSD) of the LFP from the electrode placed in region is denoted by and defined in Equation 2, where is the frequency interval of interest and is the periodogram computed with the Welch’s method [Rao2018SpectralSignals].


According to the literature on the electrophysiology of PD [Poewe2017, Tinkhauser2017BetaMedication], a noticeable abnormality is observed typically at the centre of the beta frequency band of LFP recordings from the basal ganglia of PD individuals. This frequency band corresponds approximately to the interval [13,30] Hz, although this range varies within human patients and animal model species. For the formulation of the fitness function, let a coefficient be the summation of the beta-band mean PSD plus the mean PSD of adjacent bands, composing the interval [8,50] Hz, normalised by the mean PSD of all frequencies up to Hz, as stated in Equation 3. This broader interval was defined to account for possible wider spectrum modulations in adjacent bands.


The fitness function is defined in Equation 4, where is the average value of Equation 3 calculated from the preprocessed data of all marmosets of PD condition, and is calculated considering the simulated LFP of a computational model . Notice that the healthy marmoset condition lacks readings from TH and STN regions (i.e., no electrodes were implanted in these regions). In addition, the dataset includes three PD model animals. For this reason, DE optimised parameters for mimicking the PD condition. Fitness values vary from 0, if simulated and marmoset data LFP in all brain regions substantially differ, to , if they match.


Eight brain regions are simulated, thus . PSD target values for the simulated regions StrD1 and StrD2 are drawn from marmoset LFP PSD values for PUT. Simulated TH is tuned based on the average PSD from marmoset VL and VPL, and simulated CtxRS and CTxFSI are tuned based on marmoset M1. Simulated GPe, GPi, and STN LFP PSDs are matched to the respective marmoset LFP PSDs.

The DE initial population was set to 200 individuals, whose initial parameters were drawn from a random uniform distribution in the interval

. Parameters were normalised to the ranges listed in Table 1

(i.e., the actual values set in the computational model) only at simulation time. In each DE generation, a set of 20 individuals were selected through tournaments of size two. Pairs of those selected individuals were randomly chosen, in order to generate two offspring by applying uniform crossover. This led to a child population of size 20. The mutation rate was set to 10% and followed a normal distribution

. The DE implements generational replacement with elitism, with only one elite individual of the parent population being kept, resulting in a population size of 21 individuals. Each model , where , was evolved for generations. We have performed 150 evolutionary runs, so that the highest fitness individual of each run was selected to compose the set of evolved genotypes.

2.4 Clustering

Upon completion of parameter optimisation by the DE, we investigated whether high fitness individuals had different genotypes. The rationale is that different parameter sets, even if biologically plausible, could lead to incompatible healthy and PD network dynamics [Bahuguna2017]. Considering that the fitness function was computed based on LFP values of the PD condition only, and that the healthy condition was obtained by changing the same parameters listed by Kumaravelu et al. [Kumaravelu2016ADisease]

, there was no guarantee that the genotypes evolved would lead necessarily to models that resemble the healthy and PD conditions of the animal models. For this reason, we performed a clustering analysis 

[Verma2012AMining] to the set of evolved genotypes, which we could then evaluate separately based on their spectral densities. This validation step is based on the fact that PD individuals present a peak at the beta band (13-30 Hz) when compared to healthy individuals [Tinkhauser2017BetaMedication].

Let be a set of clusters, with , where is the number of clusters, , and is the total number of genotypes within cluster . Considering to be the sample silhouette [rousseeuw1987silhouettes] of genotype with respect to , consider for all

, it is, each cluster is ordered from highest to lowest silhouette. In exploratory experiments (not shown), we investigated different clustering paradigms, namely K-means, density-based spatial clustering of applications with noise (DBSCAN), and agglomerative clustering. Based on these experiments, we opted for the K-means algorithm with two centroids (i.e.,

), because this configuration led to the highest mean silhouette score. Hence, the K-means algorithm was fed with all the individuals with the highest fitness per evolutionary run (i.e., set ), and the euclidean distances for the algorithm were computed on the 14 normalised parameters of the genotype.

2.5 Computational model spike and LFP analysis

The different clusters of genotypes were compared with respect to their parameter values, spike firing rates and LFP power spectra. For each cluster, the 50 highest fitness genotypes were chosen for the following analyses. Spectral analysis was performed by simulating , for milliseconds, in both healthy and PD conditions. Thus, for each condition, simulated LFP recordings were analysed per cluster for each condition. Since we simulated the same individuals (i.e., sets of parameters), with the same seeds for generation of random numbers, in each of the conditions (healthy and PD), the samples across these conditions were considered to be dependent. The PSDs were computed and evaluated with respect to the mean of the density spectrum per cluster, and the average power at the beta band.

For PSD analyses on the LFP of either the animal and computational models, to highlight the presence of a peak in the beta band in the PD condition, a ratio was defined as in Equation 5, where and are the mean spectral power across the PD and healthy models, respectively, for frequency . A lower threshold value was defined because, for denominators too close to zero, the ratio may lead to high values that actually has little meaning for interpretation. For the analyses with the animal models, was defined as the median power across the mean spectrum of the healthy condition. For the computational models, it was set to the percentile of the healthy spectrum.


Regarding spike dynamics, the models within each cluster were simulated for milliseconds with time step size milliseconds, always with the same seed for random number generation, and the firing frequency of all neurons was calculated in time bins, each corresponding to milliseconds.

2.6 Computational model validation

Considering that different currents, conductances, and numbers of neurons may influence the firing rate in each simulated brain region, which in turn modulates the LFP power spectra, one may conclude that even if there are different clusters, their neural dynamics are comparable because both clusters are formed by high fitness individuals. However, even if our computational model was optimised to replicate the LFP power spectra from marmoset animal models of PD, it should also mimic the power spectra from healthy marmosets (by changing selected conductances, see Section 2.1). In other words, if the computational model accurately captures the physiological phenomena responsible for the different beta-band centred LFP power spectra from PD marmoset monkeys, it should also replicate the healthy spectra (a scenario in which it was not evolved).

Therefore, we first confirmed that our marmoset animal model of PD presented frequency spectra in accordance with previous works, following Section 2.2. Then, we investigated whether the computational model would also capture this phenomena. For that, for each genotype cluster found (Section 2.4), we compared the LFP power spectra from the evolved PD computational model with that from the healthy model. This was performed by modifying a predefined set of conductances in the simulation (Section 2.1). To highlight the differences, we first analysed the ratios between the mean PSD of the PD and healthy simulated individuals from each cluster.

During evolution, fitness is given by LFP PSD in the vicinity of the beta band calculated in the whole sequence, hence it is possible that the same spectra relate to different LFP rhythms over shorter time scales. Thus, different neuronal spiking dynamics may lead to similar LFP dynamics over time. Moreover, spikes from single neurons are noisy and vary considerably over time and over repeated simulations. With large recordings, joint neuronal averages over time may hinder comprehension of neural population dynamics. Finally, one of the advantages of computational models such as the one used here is the direct access of each neuron state at any given time, but it is not trivial to interpret the dynamics of large populations of neurons over time. To clarify these issues, we studied low-dimensional neuronal trajectories for both healthy and PD computational model conditions [Cunningham2014].

To compute the neuronal trajectories, we first calculated the firing frequencies for all neurons from each simulated model in a particular cluster and condition (i.e., healthy and PD), based on the mean firing rates (MFR) taken from bins of size ms. Since the number of neurons within each region varies from 10 to 30 (see Table 1

), and there are eight regions considered for the computational model, this procedure generates time series with high dimensionality, ranging from 80 to 240, which would be difficult to visualise and analyse. To reduce the dimensionality, we employed principal component analysis (PCA)

[wold1987principal], that is, we analysed neural trajectories by projecting high-dimensional neural population activity in a 3D space using PCA of the spike MFR time series.

However, what if, instead of clearly occupying different regions in the state space, neuronal responses from the same conditions result in similar paths in the reduced dimensional space? To address this hypothesis and to compare PCA trajectories, we used Dynamic Time Warping (DTW) with Euclidean distance [Muller2007]. DTW finds the optimum non-linear alignment between two time series, hence it can estimate whether neuronal trajectories share a similar path, regardless of initial conditions. In the analysis performed, we employed the fastdtw Python package, which implements the method proposed by Salvador and Chan [salvador2007toward]. Each pair of three-dimensional time series, computed from the MFR signals and dimensionally reduced with PCA, was fed to the algorithm, which provided, as output, a scalar proportional to the dissimilarity between the two time series being compared.

More specifically, we compared the similarity of all possible pairs of neural trajectories considering all individuals within the clusters (healthy and PD dynamics). We compared all pairs of trajectories generated by individuals within the same condition (healthy or PD), which gave a measurement of how different the healthy or PD individuals are compared to each other (i.e., within-group comparison), and we compared pairs of trajectories between healthy and PD conditions (i.e., between-groups comparison).

Finally, one of the hallmarks of PD is the anomalous widespread synchronisation in the BG-T-C network. To validate our model in that aspect, we calculated the magnitude-squared coherence between nuclei and intranucleus. Based on a similar analysis performed in healthy and PD marmosets reported in Santana et al. [Santana2014SpinalDisease], we expect a widespread increase in this metric. The magnitude-squared coherence was calculated from the spike trains of neurons of each nucleus using Welch’s method with Hanning windowing without overlap and with spectral resolution of 1 Hz. The average was taken as recommended by Bendat and Piersol [Bendat2010]: the squared value of the average of the cross spectra divided by the product of the mean values of the auto spectra of each nucleus.

The value of the magnitude-squared coherence between brain regions and , defined as , was computed as in Equation 6, where is the number of neurons in region , and is the number of neurons in region , and is the cross spectrum between the spike trains from the -th neuron from region and the -th neuron from region .


Then, we considered the peak of the coherence in the 7-30 Hz band to highlight PD-related effects [Santana2014SpinalDisease]. The significance level for coherence was defined as  [Rosenberg1989], with and , because the windowing was done with 100 segments and we adopted as 95% the significance level. As the computational models have eight nuclei, an

matrix was constructed, representing the coherence between each pairs of nuclei. The median of this matrix was considered as the global coupling metric between nuclei in each simulation, because it is less sensitive to outliers than the mean.

3 Results

Based on two-seconds-long segments, computed according to the data preprocessing steps described in Section 2.2, (see Figure 3a for a sample), the PSDs of LFPs from healthy and PD marmosets were computed (see Figure 3b for the average spectrum). In all regions of the PD subjects, an increased PSD magnitude from 5 Hz to 25 Hz was observed, which is in accordance with the reported electrophysiological signatures of PD [Tinkhauser2017BetaMedication].

Figure 3: Animal data from marmoset monkeys, collected through electrodes implanted to each region of the BG-T-C circuit in a previous study [Santana2014SpinalDisease], and made available for our research. (a) Example of a two-second time-window of the preprocessed LFP of a PD-induced (i.e., 6-OHDA lesioned) marmoset. For a clearer visualisation, signals were bandpass filtered to the [8,50] Hz interval, only for this panel. (b) Top two panels show the mean power density spectra (PSD) over all segments for the healthy (blue) and PD (red) marmosets (data for each individual marmoset is included as supplementary material). For thalamic regions (i.e., VL and VPL) and STN, only 6-OHDA lesioned hemispheres are represented, since these regions were not recorded in the healthy marmoset. PSDs were normalised by the maximum PSD value for each time-window. The bottom panel shows the ratio (R) between PD and healthy PSD for each frequency (see Equation 5). To improve visualisation, is set to the median of the healthy spectrum. a.u.: arbitrary units.

From the estimated LFP power spectra from PD marmosets, the target LFP power spectra values for the computational marmoset model were computed as in Equation 3. The results, presented in Table 2, were fed to the DE fitness function (Equation 4).

0.44 0.44 0.38 0.46 0.42 0.39 0.39 0.37
Table 2: Target values obtained from the marmoset LFP data from the animal model, PD condition.

3.1 Evolutionary algorithm successfully found high fitness genotypes

After running the DE times, the resulting set of high fitness individuals (i.e., the highest fitness individual in the population at the end of each of the 60 generations at each evolutionary run) was analysed. The fitness values of all individuals were recorded at all generations of each evolution.

Figure 4 reports the best and mean individual fitness across generations, and the distribution of those values at the end of the evolutionary runs. Concretely, the best individual in a given generation is the set of parameters that led to the highest fitness value according to Equation 4. The mean individual fitness across generations refers to the average fitness of all individuals achieved at each generation.

Regarding the best individual fitness curve, results show that, at every evolutionary run, the initial population contained at least one individual with fitness value close to 6, and that value improved by approximately 1 at the end of evolution (the maximum fitness value possible is , see Equation 4). Considering the whole population, the initial average fitness was low (approximately 4.5), reaching a plateau close to 5.75 as evolution progressed. The mean fitness across individuals and the best individuals fitness have marginal improvement after generation 40, thus the DE was stopped at generations.

Figure 4: Fitness values (Equation 4) per generation of the evolutionary algorithm (box-plots regarding the runs at each generation). The genotypes (i.e., parameter sets for the free parameters elicited in Table 1) were meant to maximise , which, by definition, would be upper bounded at

. The upper panel refers to the highest fitness individuals at each evolutionary run, and the lower panel, to the mean fitness values of all individuals. (a) Box-plots of the best (upper panel) and mean (lower panel) fitness values at each generation. Outliers were represented by black diamonds. (b) Probability distribution of the best (upper panel) and mean (lower panel) fitness.

For all , we looked at the distribution of parameter values for clusters and , represented in Figure 5

. Both clusters present similar distributions for most of the parameters, either with small variance (e.g., the numbers of neurons at the cortex populations) or more uniform distributions with high variance (e.g., the number of neurons at the striatum). Other parameters, such as

and , had a clear mean peak and reduced variance in the distribution for , but a large variance for .

Figure 5: Violin plots showing the distribution of each free parameter (see Table 1) across the best individuals found at each run of the evolutionary algorithm employed for optimising this set of parameters (i.e., genotype). Although scales vary across parameters (see Table 1), all parameters were linearly scaled (i.e., normalised) to the interval at evolution time. For example, for parameters 7-14 (i.e., the numbers of neurons), a value of zero corresponds to the lower bound of the parameter interval, that is, 10 neurons. a.u.: arbitrary units.

3.2 High fitness genotypes form two clusters

A set of 150 high fitness individuals was generated by repeatedly running the evolutionary algorithm with different seeds. It is possible that high fitness individuals do not have a unique parameter distribution, and diverse parameter settings could lead to high fitness values. To investigate this issue, we performed a clustering analysis based on the evolved individuals.

Following the methods from Section 2.4, the K-means algorithm was employed to determine clusters. Figure 6 provides a radar plot representation of genotypes learnt for each cluster, and the correspondence between the mean value of each parameter and those of the rat computational model by Kumaravelu et al. [Kumaravelu2016ADisease].

Figure 6a shows 4 representative genotypes , chosen based on the highest silhouettes with respect to each cluster. For comparison, the parameters from the rat model [Kumaravelu2016ADisease] are superposed with the mean values between all individuals from both clusters in Figure 6b. This representation highlights substantial differences between clusters. For instance, the is at its maximum value in , while it shows a much lower value for . On the other hand, the number of neurons at the GPe is higher in than in .

Figure 6: Radar representations of the genotypes (i.e., sets of parameters, see Table 1) from individuals at each cluster obtained by applying the K-means algorithm, applying these parameters as features for the clustering technique. Although scales vary across parameters (see Table 1), all parameters were linearly scaled (i.e., normalised) to the interval at evolution time. For example, for parameters 7-14 (i.e., the numbers of neurons), a value of zero corresponds to the lower bound of the parameter interval, that is, 10 neurons. The first row represents cluster and the second row cluster . (a) Four individuals with the highest silhouettes with respect to each cluster. Data at the left refers to the fitness computed as in Equation 4. (b) Comparison between the parameters of the rat model and the mean values from each cluster. As in Figure 5, parameter values were scaled to the ranges shown in Table 1, except parameter 4 () of the rat model, whose original value is 1.0 .

3.3 Healthy and PD spectral signatures from computational model resembles those from marmoset monkeys

Regarding the spectral analyses of simulated sessions of the computational model, we employed the same procedure for normalisation as we did for the spectra of the animal model (see Figure 3), that is, we normalised each data segment by the maximum value. The sample signals of Figure 7a, shown as an example, were bandpass-filtered to the same range as in Figure 3a to the interval [8-50] Hz. The mean spectral power and the ratio are shown for the healthy and PD conditions for each cluster in Figure 7b (see Equation 5).

In , results show higher magnitudes of most frequencies up to 50 Hz for PD models, a fact that is less visible for . The mean PSD ratio from genotypes is close to 1 regardless of frequency range and brain region, whereas genotypes show prominent peaks in beta frequencies. A detailed analysis of box-plots (Figure 7c) confirm the significant differences in the beta band for cluster only. Considering the animal spectra (Figure 3b), in which we observe a significant difference in the beta band of the LFP, results displayed in Figure 7c confirm that spectral signatures from genotypes in resembles those from marmoset monkeys. Notice that the LFP mean PSDs from the computational model (Figure 7b) has a different shape compared to that from the animal LFP (Figure 3b), but the spectral signature is similar in both healthy and PD conditions and resemble those from marmoset monkeys. This can be explained by the relatively small number of neurons simulated in the computational model [Parasuram2016]. Therefore, for the forthcoming analyses, only will be considered.

Figure 7: Extracellular activity simulated by the computational models resulting from the parameters optimised (i.e., genotypes computed with the evolutionary algorithm), modelled as local field potentials (LFP) at the centre of the regions involved in the BG-T-C circuit. The clusters and were computed by applying the K-means technique directly to the genotypes, hence were not influenced by the neurophysiological activity simulated. (a) Example of simulated LFPs for the highest silhouette evolved individual from cluster , PD condition. For a clearer visualisation, signals were bandpass filtered to the [8,50] Hz interval, only for this panel. Compare with Figure 3a. (b) Mean PSD for healthy (blue) and PD (red) conditions from the 50 models with the highest silhouette of each cluster, normalised by maximum PSD value for each time-window, followed by the ratio between PD and healthy PSD for each frequency (see Equation 5). To improve visualisation, is set to the percentile

of the healthy spectrum. (c) Box-plot regarding the beta band (13-30 Hz) of the LFP from the 50 models with the highest silhouette of clusters 1 (left) and 2 (right). Outliers were represented by black diamonds. Unpaired t-tests were applied to evaluate statistical significance against the null hypothesis that H and PD values are drawn from the same underlying distribution (p-value notation:

; ; ; ; ). a.u.: arbitrary units.

3.4 Spike activity from healthy models are significantly different from those of PD models

Regarding spike activity, the marmosets’ dataset was not provided with a representative set of spike trains from all regions of the circuit, hence they were not a suitable ground-truth for validating the activity from the computational model. For this reason, the spikes synthesised by the computational model were analysed based on evidence from the literature [Prescott2006AProcessing].

First, we assessed the differences in mean firing rates (MFR) between the healthy and PD conditions for the marmoset-based computational models in cluster . Figure 8a shows the simulated MFR in each brain region for ms, considering the 50 models in with the highest silhouette with respect to the cluster. Results indicate a counter-intuitive relationship between the MFR and the LFP power spectra observed in Figure 7c. Consider, for instance, the GPe and GPi. Both regions show a higher beta-band LFP magnitude in PD condition, but while GPi MFR in PD condition is higher than that from healthy condition, GPe MFR is the opposite.

From Figure 8b and Figure 8c, we observe that neuronal trajectories are intertwined, with no clear difference in the reduced-dimension state space. This is justified by the relatively mild, though statistically significant, differences in MFR (Figure 8a). As described in Section 2.5, neuronal trajectories were compared with DTW in three scenarios: healthy vs healthy models (HxH), PD vs PD models (PDxPD), and healthy vs PD models (HxPD). As , the number of pairs from which the DTW was computed was for each scenario. The results from this analysis are shown in Figure 8d, in which the scalar outputs of the DTW algorithm are considered for all possible pairs within groups, for the HxH and PDxPD comparisons, or between groups, for the HxPD comparisons. Since two trajectories generated by the same individual were not compared on any of the analyses, we have computed statistical significance using unpaired tests, differently from the remaining analyses in the paper.

Figure 8: Firing rates and dynamics regarding the spike activity simulated with the computational models derived by the parameter sets from cluster . Simulations were ran for . (a) Mean firing rates for each region in cluster (means and standard deviations). (b) Projection of three principal components of the most representative individual (i.e., highest silhouette) of cluster , where , and are the principal components with the highest variance. (c) Representation of those components using contour lines. (d) Box-plot of the DTW between the dynamics of one simulation of all genotypes belonging to cluster

. All simulations were performed with the same seed for the generation of random numbers. Higher DTW values mean that the pairs of trajectories being compared are less similar to each other. Unpaired t-tests were applied to evaluate statistical significance in (a), against the null hypothesis that H and PD MFR values at each region are drawn from the same underlying distribution, and in (d), against the null hypothesis that a given pair of DTW vectors is drawn from the same distribution as each of the others (p-value notation:

; ; ; ). a.u.: arbitrary units.

Trajectories from the HxH scenario were statistically more similar than trajectories from the other conditions. Thus, the intertwined trajectories observed in PCA (Figure 8b and Figure 8c) in fact relate to significant differences between healthy and PD neuronal dynamics. Interestingly, PDxPD trajectories differ more than those from HxH, which can be interpreted as a less homogeneous, regarding neuronal dynamics, genotype to phenotype mapping.

3.5 Healthy and PD spike coherence from the computational model resembles that from marmoset monkeys

To conclude our model validation, we selected the top five genotypes with highest silhouette from cluster and calculated the magnitude-squared coherence (MSC) within and between each simulated brain region (Section 3.5) for healthy and PD conditions. Results revealed that computational models ran in the healthy condition provided a lower peak MSC in the 13-30 Hz band when compared to that from the PD condition (Figure 9a), with two important observations: genotype I has higher peak MSC when compared to the other 4 genotypes in the healthy condition, and genotypes II and III have a lower widespread peak MSC when in the PD condition compared to that from other genotypes in the same condition. Statistical analysis confirmed the significant differences in all five genotypes when comparing the global coupling metric (see Section 2.6 and Equation 6) between healthy and PD conditions (Figure 9b), that is, PD models present a higher widespread coherence in the 13-30 Hz band than that observed in healthy models.

Figure 9: Coherence analyses computed for the spike activity of the five parameter sets, optimised through the evolutionary algorithm, with the highest silhouettes with respect to cluster . These parameter sets could were employed to derive the healthy and PD computational models employed to these analyses. (a) Peak magnitude-squared coherence (MSC) in the 13-30 Hz band within and between each simulated brain region for the top five genotypes with highest silhouette from cluster . Only connections whose peak MSC values are above significance level are shown. (b) Global coupling metric (median value of the MSC matrix) between brain regions for healthy and PD conditions (see Section 2.6 and Equation 6). (p–value notation: )

4 Discussion

Marmoset monkeys are prominent in neuroscience research [Cyranoski2009, Kishi2014, Mitchell2015, Miller2016]. Although there are anatomical and physiological differences between BG-T-C circuit in rodents and primates, neurophysiological data from rodents are far more available than from primates. Considering that the structure of the BG-T-C circuit presents similar characteristics among all vertebrates [Koprich2017AnimalDevelopment], we assumed that the rat model presented by Kumaravelu et al. [Kumaravelu2016ADisease] was a suitable starting point to build a computational model of those structures in primates. The core hypothesis was that, by keeping the same brain regions and connectivity patterns of the rat model and modifying a set of parameters, the computational model could reproduce neural dynamics of healthy and PD marmoset conditions.

Our dataset comprised simultaneous LFP recordings from regions of the BG-T-C network and power density spectra (PSD) analysis revealed significantly higher 13 to 30 Hz LFP PSD magnitudes for PD marmosets on all regions. This result might be interpreted cautiously, given that one healthy marmoset is being compared to three PD marmosets. Also, results refer to a broad range of frequencies, hence different interval choices may influence the analysis. Nonetheless, one would expect a widespread significant increase in LFP power centred in (but not limited to) the beta band in PD affected brains [Santana2014SpinalDisease, Tinkhauser2017BetaMedication].

Regarding the MFR results from the computational model (Figure 8a), there are significant differences between the healthy and PD conditions. Single-neuron firing rates vary considerably depending on animal species, whether the animal is fully awaken, engaged in behavioural tasks, or anaesthetised [Hardman2002, VanAlbada2009, McGregor2019CircuitDisease]). Data from human subjects, even though scarce, are in line with animal results [Du2018]. Moreover, there is a great neuronal diversity within the BG-T-C network, both in terms of neuronal physiology and connectivity, which have been shown to have a non-trivial relationship with field potentials [Sharott2009, Bolam2000, Holgado2010, Benhamou2012]. Our model partially takes into account this diversity, nevertheless the reported MFR are in agreement with the literature: comparing PD with H conditions, a higher MFR in GPi, STN, and Str, and a lower MFR in GPe, TH, and CTX.

The data-driven modelling strategy adopted in this paper is consolidated in computational neuroscience literature [Nowke2018], but often leads to multiple models fitting a particular data set [Bahuguna2017]. Therefore, model optimisation should be followed by a model selection phase. We clustered high fitness solutions with respect to evolved parameters and obtained two clusters, and found two clear sets of parameters that reproduce the increased beta-band oscillations observed in PD marmosets [Santana2014SpinalDisease]. However, when perturbing the model to shift from PD to healthy dynamics, only one of the clusters fitted the marmoset data. Notably, we evolved solutions based on LFP data but computational model firing rates resemble those reported in previous works [VanAlbada2009, Li2015, Du2018]. Nevertheless, as data becomes available, future works should explore different fitness functions based on single-neuron activities or other features of LFP. Lastly, in this context, our simulated neurons are formed by a single cylindrical compartment, thus future works should consider using neurons with more and more complex compartments and connections, possibly including multiple dendritic branches and active ionic channels. This would lead to more realistic simulated LFP signals [Parasuram2016], but at the expense of heavier computing resources.

One of the great challenges in neuroscience is to link the activity of large neural populations to motor and cognitive behaviours. One strategy is to study the intrinsic high-dimensional dynamics of neural populations from its low-dimensional dynamics given by time-varying trajectories [Cunningham2014, Sussillo2014], thus emphasising circuit over single-neuron function. For example, Humphries et al. [Humphries2018] showed that neural low-dimensional dynamics given by PCA of neuronal activity can explain Aplysia rhythmic movement control and propose that only the low-dimensional dynamics are consistent within and between nervous systems. Also, the shape and amplitude of neural trajectories can explain different behavioural outcomes [Gamez2019]. Combining PCA and DTW, we found that neural trajectories from high-fitness models are more similar in healthy conditions than in PD conditions. This is in line with results from Russo et al. [Russo2019], who demonstrated, using computer simulations, later confirmed by data from the supplementary motor area in monkeys, that low trajectory divergence is essential in neural circuits involved in action control. PCA is a simple, established method for dimensionality reduction, but other computational tools tailored to neuronal data, such as Gaussian-Process Factor Analysis (GPFA) [Yu2009] and jPCA [Churchland2012], should be considered in further analyses. Another possible approach is to use more advanced machine learning methods to identify PD-related features from neural data, likewise Ranieri et al. [Ranieri2020]

, who employed a deep learning framework to unveil PD features from marmoset data.

Finally, as part of our model validation, we assessed functional coupling within and between simulated brain regions by means of coherence between spike trains. In contrast to structural coupling, characterised by physical neuronal connections, functional connectivity is an emergent phenomenon commonly linked to synchronisation in neural rhythms in diverse spatiotemporal scales and is the basis of neural communication and cognitive processing [Singer1999, Buzsaki2010, Fries2015, Lakatos2019]. Several neural disorders, including PD, present a disruption in functional connectivity [Uhlhaas2010, Mathalon2015, Halje2019]. In particular, Santana et al. [Santana2014SpinalDisease] showed that 6-OHDA marmoset models of PD have a widespread coherence peak in the beta band when compared to healthy individuals. Our computational model is in line with this result, which is relevant not only as further evidence of its biological plausibility, but also because one of the established therapies to alleviate PD motor symptoms is the use of deep brain stimulation (DBS) [Poewe2017]. Thus, we believe that the work presented here can be used to test hypotheses that employ DBS. For instance, Romano et al. [Romano2020] performed a comprehensive analysis of frequency-dependent effects of DBS on the same model that we used here, tuned for rodent data [Kumaravelu2016ADisease], and found that neural oscillatory modulations were similar to those observed in electrical brain and spinal cord stimulation of primates [Wang2018, Santana2014SpinalDisease].

Certain simplifications inherent to our approach may be worth a mention, as they may serve as inspiration for improvements in future research. In our work, LFP generation followed the method described in Parasuram et al. [Parasuram2016], and implemented in NetPyNE, which does not consider the influence of sinks. Despite being a simplification, the method has been able to reproduce features of real LFP data, and is computationally feasible. In this approach, LFP peaks and valleys are directly related to transmembrane ionic currents from each neuronal source, which in turn relate to neuronal firing rates, and electrode position. As we have assigned coordinates to the simulated electrodes corresponding to the centre of each simulated region, we can assume that simulated LFP dynamics is due to altered spiking activity in multiple neuronal sources from different brain regions.

In our work, likewise Kumaravelu et al. [Kumaravelu2016ADisease] and previous seminal BG-T-C modelling works such as Humphries et al. [Humphries2006], and van Albada and Robinson [VanAlbada2009], we did not model any structural synaptic plasticity mechanisms. Our synapses were modelled as bi-exponential and alpha synapses, including transmission delays. Nevertheless, as model dynamics unfold, functional plasticity mechanisms may take place in the sense that the closed-loop, recursive network architecture could lead to single neurons and brain regions whose electrical activity are sensitive to past network states. In fact, the depletion of dopamine, one of the hallmarks of PD, affects structural and functional plasticity. Our model considers the loss of dopaminergic neurons (see Section 2.1, for a complete description), thus we believe that the model is suited for the investigation of functional plasticity phenomena. This analysis is beyond the scope of our work, but the reader can relate the change in oscillatory neural dynamics we described to different functional states. For instance, Humphries et al. [Humphries2006] show that action selection in the BG is closely linked to oscillatory activity. In future works, we plan to use this model in a neurorobotics context, in which sensory inputs and motor responses can be used to highlight functional plasticity mechanisms differences between healthy and PD states.

5 Conclusion

Computational models are invaluable tools for advancing our knowledge on the neural dynamics of our brain, either under healthy conditions or with neurological disorders. In this work, we created the first computational model of PD based on data from Marmoset monkeys both in healthy and parkinsonian conditions. Our data-driven approach used simultaneous, multisite electrophysiological recordings from healthy and 6-OHDA+AMPT marmoset models of PD. Even though the physiopathology underlying PD share similarities across vertebrate species, there are important, species-specific differences in the anatomy and neural dynamics of the BG-T-C circuit. Hence, the design of a primate computational model of PD is of paramount importance.

Electrophysiological datasets from animal models often do not include comprehensive biophysical data such as single-neuron membrane conductances and neuronal cell densities. These parameters are central for building a biophysical computational model. Thus, to address this gap, we implemented a DE to search the multidimensional model parameter space for solutions that could reproduce features of the animal LFP recordings. Our model was based on a well known rat model of PD [Kumaravelu2016ADisease]. The main novelty aspects of our model are: 1) we use a marmoset monkey BG-T-C electrophysiological database; 2) we added LFP simulations to the model, in addition to spike dynamics; and 3) we developed a DE-based optimisation to search for unknown parameters. With this framework, we were able to reproduce several of the previously reported PD electrophysiological biomarkers observed and recorded from the Marmoset monkeys.

Our computational model present beta-band LFP power spectra differences between the healthy and the PD conditions, which Wang et al. [wang2016subthalamic] also found in human patients with dystonia. This is in line with a body of literature that shows that beta-band LFP modulations are not a PD-specific biomarker (see Poewe et al. [Poewe2017] and references therein). Although our model is focused on PD, the electrophysiological features we use are known to be related to other neural disorders and thus should not be considered as exclusive to PD. Also, based on the study of Wang et al. [wang2016subthalamic], we suggest as future work to conduct the phase amplitude coupling in the STN experiment in our computational model.

Most PD computational models do not consider brain-body-environment interactions. Embodied cognitive science studies have provided solid evidence that neural activity is shaped by such interactions [Engel2001, Pfeifer2006, Badcock2019, Musall2019]. In PD and other neural disorders, body-environment interactions influences motor control [Snijders2010, Santos2017], but its impact on neural dynamics remains unclear. Moreover, the BG-T-C neuronal network is clearly related to action selection and decision making [Humphries2006, Mink2018, Suryanarayana2019]. Therefore, we believe that our marmoset-based computational model associated with robotics may offer an alternative approach to elucidate the mechanisms underlying brain-body-environment interactions in PD [Gurney2004, Prescott2006, Krichmar2018, Pronin2021NeuroroboticReview]

. A possible approach would be to employ this computational model in a sensorimotor loop based on visual inputs from video cameras and motor outputs to actuators such as the robot’s motors. In this scenario, computer vision algorithms would transform the images into stimuli for the computational model, so that the resulting currents and action potentials would be used to generate perturbations that would govern the behaviours of the actuators. The resulting framework could become a new tool for studying the underlying mechanisms of PD and the effects of different interventions regarding the simulated circuit.


This work is part of the Neuro4PD project, granted by Royal Society and Newton Fund (NAF\R2\180773), and São Paulo Research Foundation (FAPESP), grants 2017/02377-5 and 2018/25902-0. Moioli and Araujo acknowledge the support from the National Institute of Science and Technology, program Brain Machine Interface (INCT INCEMAQ) of the National Council for Scientific and Technological Development(CNPq/MCTI), the Rio Grande do Norte Research Foundation (FAPERN), the Coordination for the Improvement of Higher Education Personnel (CAPES), the Brazilian Innovation Agency (FINEP), and the Ministry of Education (MEC). Romano was the recipient of a master’s scholarship from FAPESP, grant 2018/11075-5. Elias is funded by a CNPq Research Productivity Grant (312442/2017-3). This research was carried out using the computational resources from the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP, grant 2013/07375-0. Additional resources were provided by the Robotics Lab within the Edinburgh Centre for Robotics, and by the Nvidia Grants program.