Parkinson’s disease patient-specific neuronal networks carrying the LRRK2 G2019S mutation unveil early functional alterations that predate neurodegeneration

ARTICLE OPEN Parkinson’s disease patient-specific neuronal networks carrying the LRRK2 G2019S mutation unveil early functional alterations that predate neurodegeneration G. Carola 1,2,14 , D. Malagarriga 1,2,14 , C. Calatayud 1,2,3,14 , M. Pons-Espinal 1,2,14 , L. Blasco-Agell 1,2 , Y. Richaud-Patin 3 , I. Fernandez-Carasa 1,2 , V. Baruffi 1,2 , S. Beltramone 1,2,4 , E. Molina 5 , P. Dell’Era 4 , J. J. Toledo-Aral 6 , E. Tolosa 7 , A. R. Muotri 8 , J. Garcia Ojalvo 9 , J. Soriano 10,11 ✉,A.Raya 3,12,13 ✉and A. Consiglio 1,2,4 ✉ A deeper understanding of early disease mechanisms occurring in Parkinson’s disease (PD) is needed to reveal restorative targets. Here we report that human induced pluripotent stem cell (iPSC)-derived dopaminergic neurons (DAn) obtained from healthy individuals or patients harboring LRRK2 PD-causing mutation can create highly complex networks with evident signs of functional maturation over time. Compared to control neuronal networks, LRRK2 PD patients’networks displayed an elevated bursting behavior, in the absence of neurodegeneration. By combining functional calcium imaging, biophysical modeling, and DAn-lineage tracing, we found a decrease in DAn neurite density that triggered overall functional alterations in PD neuronal networks. Our data implicate early dysfunction as a prime focus that may contribute to the initiation of downstream degenerative pathways preceding DAn loss in PD, highlighting a potential window of opportunity for pre-symptomatic assessment of chronic degenerative diseases. npj Parkinson’s Disease (2021) 7:55 ; https://doi.org/10.1038/s41531-021-00198-3 INTRODUCTION Parkinson’s disease (PD) is the most common neurodegenerative movement disorder, with an estimated prevalence in industrialized countries of 0.3% in the general population, which increases to 1.0% in people older than 60 years and to 3.0% in people older than 80 years 1 . Clinically, PD is characterized by classical motor syndrome linked to a progressive loss of dopamine-containing neurons (DAn) in the substantia nigra pars compacta, and disabling non-motor symptoms related to extranigral lesions. Current therapies for PD are symptomatic and do not limit the progression of disability with time. It has been proposed that early intervention might slow down or even stop disease progression, by preserving neurons from the undergoing irreversible neurodegeneration 1,2 . However, early treatment relies on early diagnosis, which unfortunately is especially complicated in the case of PD. Current diagnostic modalities in PD are based on the presence of motor symptoms, a stage at which up to 70% of DAn have been lost 3 . Even though pre-motor symptoms are known to precede clinical diagnosis of PD by as much as a decade, they are rather unspecific and unsuitable as stand-alone biomarkers of the disease 4 . Therefore, the identification of early diagnostic or progression markers of PD represents an urgent medical need. Although the majority of PD cases are of unknown cause, socalled idiopathic PD, around 5% have been shown to have a genetic basis, with mutations in the LRRK2 gene accounting for the largest number of patients of familial PD 5 . Interestingly, LRRK2 polymorphisms are also considered a relevant genetic determinant for sporadic PD 6 , and LRRK2 function appears dysregulated in sporadic cases of PD 7 , even in the absence of LRRK2 mutations/ polymorphisms. These findings, together with the fact that PD associated with mutations in LRRK2 (L2-PD) is clinically indistinguishable from sporadic PD, position LRRK2 as an essential player for understanding both genetic and idiopathic PD 8 . LRRK2 is a highly complex protein with multiple enzymatic domains, involved in a variety of intracellular signaling pathways and cellular processes such as cytoskeleton dynamics, vesicle trafficking and endocytosis, autophagy, reactive oxygen species, mitochondrial metabolism, and function of immune cells. However, the exact physiological role of LRRK2 and its implication for PD pathogenesis remains unknown 8 . Of especial relevance for the investigation of early disease markers, transgenic mouse models of L2-PD display, before any events of neurodegeneration, an abnormally elevated excitatory activity and altered spine morphology in dorsal striatal spiny projection neurons 9 . Moreover, experimental models for other neurodegenerative conditions such as Alzheimer’s disease 10,11 and amyotrophic lateral sclerosis 12 , have been shown to exhibit neuron hyperexcitability before the disease onset. It has also been demonstrated that the combination of PD with dementia often correlates with a disruption of both functional and effective connectivity in the cortex 13 . In contrast, the association of PD with depression correlates with disrupted functional connectivity 1 Department of Pathology and Experimental Therapeutics, Bellvitge University Hospital-IDIBELL, Barcelona, Spain. 2 Institute of Biomedicine (IBUB) of the University of Barcelona (UB), Barcelona, Spain. 3 Regenerative Medicine Program, Bellvitge Biomedical Research Institute (IDIBELL), and Program for Clinical Translation of Regenerative Medicinein Catalonia (P-CMRC), Hospital Duran i Reynals, Barcelona, Spain. 4 Cellular Fate Reprogramming Unit, Department of Molecular and Translational Medicine, University of Brescia, Brescia BS, Italy. 5 Sanford Consortium for Regenerative Medicine, UC San Diego, La Jolla, CA, USA. 6 Instituto de Biomedicina de Sevilla (IBiS), Hospital Universitario Virgen del Rocío/CSIC/Universidad de Sevilla, Sevilla, Spain. 7 Department of Neurology, Hospital Clinic de Barcelona, Institut d’Investigacions Biomediques August Pi i Sunyer (IDIBAPS), University of Barcelona (UB), Barcelona, Spain. 8 Departments of Pediatrics and Cellular & Molecular Medicine, University of California San Diego, La Jolla, CA, USA. 9 Department of Experimental and Health Sciences, Universitat Pompeu Fabra, Parc de Recerca Biomèdica de Barcelona, Barcelona, Spain. 10 Departament de Física de la Matèria Condensada, Universitat de Barcelona, Barcelona, Spain. 11 Universitat de Barcelona Institute of Complex Systems (UBICS), Barcelona, Spain. 12 Centre for Networked Biomedical Research on Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), Madrid, Spain. 13 InstitucioCatalana de Recerca i Estudis Avancats (ICREA), Barcelona, Spain. 14 These authors contributed equally: G. Carola, D. Malagarriga, C. Calatayud, M. Pons-Espinal. ✉email: jordi.so[email protected]; [email protected]; [email protected] www.nature.com/npjparkd Published in partnership with the Parkinson’s Foundation 1234567890():,; between the median cingulate cortex and the prefrontal cortex and cerebellum 14,15 . The development of induced pluripotent stem cell (iPSC) technologies enables the generation of patient-specific, diseaserelevant, cell-based experimental models of human diseases. Importantly, iPSC-based models can recapitulate some of the earliest signs of disease, even at pre-symptomatic stages 12,16 .In this study, we used an experimental platform based on DAnenriched neuronal cultures derived from L2-PD patients, their gene-edited isogenic counterparts, or from healthy individuals. Such cultures formed active neuronal networks, the functionality of which was analyzed by calcium imaging. After multiple iterations of experimental characterization and biophysical modeling of neuronal network behavior, we could identify early alterations in PD neuronal function that were not present in control networks, and that predated the onset of neuron degeneration. RESULTS Generation and characterization of iPSC-derived DA neurons A total of seven iPSC lines representing L2-PD patients and healthy aged-matched controls, along with gene-edited counterparts and fluorescent TH reporters, were used for the current studies (see Table 1and “Methods”for further details). Some of these iPSC lines have been previously generated and fully characterized in our laboratories 17–19 , whereas two additional TH reporter lines were generated for this study (Supplementary Fig. 1). iPSC differentiation toward DAn fate was performed using a modified version of the previously established midbrain floorplate protocol 20 , which enabled the maintenance of differentiated cells for up to 10 weeks 21 .Briefly, we first cultured iPSCs on Matrigel with mTeSR medium until they reached 80% confluence, then we induced specification toward the ventral midbrain (VM) fate using a combination of knockout serum medium and neural induction medium (Fig. 1a). At day 12 postplating (D12) the cells exhibited a homogeneous morphology and marker profile of VM floor-plate progenitors, expressing FOXA2+/LMX1A+and the VM NPC markers such as OTX2 and EN1 paired with the neuroectodermal stem cell marker Nestin (Fig. 1b). These progenitors were then cultured in neuronal differentiation medium supplemented with growth factors including BDNF, GNDF, TGF-β,andDAPT,withtheaimtofoster neuronal differentiation and survival (Fig. 1a). At D50, the majority of the cells were positive for the dendritic marker MAP2, and quantitative immunolabeling for tyrosine hydroxylase (TH) and FOXA2 revealed that 30–40% of them were also committed to DA neuronal fate (Fig. 1c, d). By the same day of differentiation (D50), about 35% of TH+ neurons also expressed the A9 domain-specific marker G proteinactivated inward rectifier potassium channel 2 (GIRK2) (Fig. 1c, d), and displayed Dopamine Transporter (DAT), the essential marker of mature DAn (Fig. 1e). TH+neurons also showed expression of pre-synaptic marker Synapsin and post-synaptic marker PSD95, indicating their capability to form synapses (Fig. 1f). Under these conditions, controls (CTR) and PD-iPSCs gave rise to VM DAn that were morphologically homogeneous and showed the expected features of mature VM DAn, including complex dendritic arborization (Fig. 1g). Although the protocol used was comparably effective in all iPSC lines analyzed, we found some variability in the number of VM DAn across lines at D35, ranging from 10 to 20% VM DAn ratio concerning all differentiated cells. This ratio did not depend on the presence or type of disease (Fig. 1h). Instead, it seemed related to the specific iPSC clone used and to the evolution at the early stages of differentiation. After establishing equivalent DAn-enriched cultures from both control and patient iPSC lines, we evaluated whether there were any differences in cell viability between DAn derived from controls and patient iPSC lines. We first noted that control and L2-PD iPSCderived DAn were indistinguishable, based on immunofluorescence (IF) images (Fig. 1g) up to 80 days in culture, the latest timepoint investigated. By counting the number of TH-expressing cells (Fig. 1h) we found no decline when cultured over time up to D80, strongly suggesting that DAn are not degenerating under these conditions (Fig. 1g, h). We also found no differences in the percentage of cells with pyknotic nuclei in patient lines compared to control lines (data not shown; D50: 12–15% in all the lines; D80: 15–20% in all the lines). We next analyzed whether changes in transcriptome profiles suggestive of neurodegeneration appeared under these conditions. For this purpose, we measured the expression of a panel of 770 genes relevant for this process (NanoString Human Neuropathology Panel) in neuronal cultures differentiated from control and L2-PD iPSC at D50. Independent differentiation experiments from the same iPSC lines displayed highly similar gene expression profiles (Fig. 1i, Supplementary Fig. 2a), highlighting the robustness of the differentiation protocol. A comparison of gene expression levels between control and L2-PD cultures identified 157 differentially-expressed transcripts with ≥2-fold change (50 upregulated and 107 downregulated), but no differences reached statistical significance with p-Adj ≤0.1 (Supplementary Dataset), nor did they correlate with PD-related hallmarks (Supplementary Fig. 2b–d). Thus, even though previous studies have found overt signs of neurodegeneration upon long-term culture of DAn differentiated from L2-PD iPSC 17,18,22,23 , under the culture conditions used in the present work (that include neurotrophic factors), they appeared to form overall healthy neuronal cultures similar to the ones generated by control iPSC. Characterization of neuronal activity Calcium activity recordings were then performed across all the seven independent lines to examine the functional maturity of the iPSC-derived neuronal networks after 35, 50, and 80 days of differentiation (D35-D50-D80) (Fig. 2a). We noted that our calcium fluorescence assay enabled tracking the behavior of ~500 neurons in the field of view with a high spatial and temporal resolution, allowing us to resolve single cells and their dynamic interactions. Table 1. Summary of the healthy controls and patients used in this study. Cell line code Status Code used in this study Sex Age at donation Age at onset LRRK2 mutation Isogenic control SP-11 Control CTR F 48 No SP-11TH Control CTR TH F 48 No SP-12wt/wt Isogenic Control isoPD1 F 63 G2019S LRRK2 G2019S corrected SP-12 L2-PD PD1 F 63 49 G2019S SP-12TH L2-PD PD1 TH F 63 49 G2019S SP-12THwt/wt Isogenic Control Iso PD1TH F 63 G2019S LRRK2 G2019S corrected SP-13 L2-PD PD2 F 68 57 G2019S G. Carola et al. 2 npj Parkinson’s Disease (2021) 55 Published in partnership with the Parkinson’s Foundation 1234567890():,; Sharp increases in the fluorescence traces (Fig. 2a) revealed spontaneous neuronal activations, which were analyzed to extract the onset times of elicited action potentials. With these data we first computed the average neuronal activity along culture time (Fig. 2b) and observed that all lines evolved similarly, indicating that any anomalies in PD networks would depend on the structure of the activity patterns, rather than on the strength of activity itself. Thus, we next computed the Global Network Activity (GNA), Fig. 1 Generation of ventral midbrain (VM) dopaminergic neurons from human iPSCs. a Scheme showing VM dopaminergic neurons differentiation protocol. bRepresentative immunofluorescence (IF) images of CTR (SP11) and PD1 (SP12) NPCs at day 12 of the differentiation process. iPSC-derived neural cultures express floor plate progenitor markers, such as Lmx1A, FoxA2, Otx2, Nestin, and EN1. cRepresentative IF images of CTR (SP11) and PD1 (SP12) differentiated neuronal cultures expressing neuronal markers (MAP2, TH) and midbrain-type DA markers (FoxA2 and Girk2) at day 50. Scale bar is 50 μm. dPercentage of TH/FoxA2 and TH/Girk2 vmDAn at day 50 of all the lines [CTR (SP11), CTR TH (SP11 TH), gene-edited isoPD1 (SP12 wt/wt), PD1 (SP12), PD1 TH (SP12 TH), PD2 (SP13)]. A number of independent experiments n=3. e–fRepresentative images of neuronal cultures at D50 expressing specific markers of maturation: eDAT (dopamine transporter) (scale bar, 50 μm); and fPSD95 (post-synaptic marker, red) and synapsin (synaptic marker, green). Orthogonal views show colocalization (scale bar, 10 μm). gRepresentative images of CTR (SP11) and PD1 (SP12) neuronal culture at D35-D50 and D80 showing the expression of TH and MAP2 markers at the different timepoints. hPercentage of TH+cells over DAPI at the three different timepoints (D35-D50 and D80). iHeatmap of gene expression profiles of neuronal cultures of CTR (SP11), PD1 (SP12), PD2 (SP13), and gene-edited isogenic PD1 line (isoPD1) at D50, with dendrograms showing the strong similarities between independent experiments and the absence of clustering between control (CTR and isoPD lines, light brown) and PD (PD1 and PD2 lines, blue) conditions. Shown are transcripts with ≥2-fold change in expression, grouped as upregulated (50 transcripts, green bar) or downregulated (107 transcripts, purple bar). No statistically significant differences were found at p-Adj ≤0.1 when comparing control and PD conditions. Each type of culture was analyzed a minimum of 2 independent times, except for the gene-edited isogenic PD1 line, for which only one RNA preparation passed the quality control and could be analyzed. G. Carola et al. 3 Published in partnership with the Parkinson’s Foundation npj Parkinson’s Disease (2021) 55 defined as the fraction of neurons in the network that coactivated in a time window of 1 s. As shown in Fig. 2c, the GNA captures the level of neuronal synchronization present in the raster plot. For control cultures and also for the rescued isogenic PD line (left and central panels), collective activity encompassed synchronous activations of moderate size, between 10 and 40% of the network. In contrast, PD neuronal cultures (Fig. 2c, right panels) displayed a two-state dynamics, with strong whole-network synchronous events combined with quiescent intervals. Importantly, the distinct GNA patterns of control and PD neuronal cultures hint at the existence of intrinsically different network mechanisms in the two systems, which orchestrate a markedly different collective behavior. The GNA also informed that average neuronal activity by itself was not sufficiently informative to reveal alterations due to disease. To quantify the differences between cell lines exposed by the GNA analysis, we next analyzed the amplitude of the GNA events and extracted those that exceeded a given threshold (Fig. 2c, red dotted line), termed ‘extreme event’. We observed that the ‘ratio of extreme events’, i.e., the number of large bursting episodes relative to all detected episodes, was much higher in PD lines than in CTR or in genetically-corrected isogenic control (isoPD) lines, particularly at late stages of maturation (D80) (Fig. 2d), which could render them more susceptible to stress by creating an overly rigid, synchrony-locked network despite their continued viability in culture. Functional connectivity of control and PD neuronal networks The different structure of global activity of CTR, isoPD, and PD lines hints at the existence of distinct functional connectivity traits between them. To shed light on network functionality, we used transfer entropy (see “Methods”for details) to compute the functional connectivity among all pairs of active neurons in a given network. As shown in Fig. 3a (top panels), the functional networks displayed abundant connections, with a combination of short-range and long-range links that extended across the network. Even though the functional networks seemed similar in spatial organization, we observed significant differences between control and PD networks. The first difference was a lower density of connections in the PD line, suggesting an overall degradation of functionality. A second difference concerns the structure of functional communities, i.e., groups of neurons that tend to connect among themselves more strongly than with the rest of the network. As shown in the bottom panels of Fig. 3a, functional communities were abundant in CTR and isoPD cultures (highlighted blue boxes), indicating cross-talking of neurons in small groups, in agreement with the moderate size of collective activations. For PD cultures, however, the communities were Fig. 2 Network dynamics of vmDA neurons. a Representative image of a bright field recording at D80 (isoPD1) using the calcium imaging assay. Scale bar is 100 μm. Regions of interest were manually selected for each neuron (diameter 10 μm; color circles) to obtain the normalized calcium fluorescence time series of spontaneous activity, DFF (%) 100 ðFF0Þ=F0, with F 0 the fluorescence signal of the neuron at rest. The green box illustrates the fast rise of fluorescence upon activation that procures the spiking onset time (arrowhead). The dashed black boxes illustrate coordinated neuronal activity that shape network bursts when several neurons are involved. bAverage neuronal activity along maturation for the different cell lines. Data points are expressed as mean±SD. Trend lines are linear regressions. The number of cultures used in each condition and timepoint were: D35 (CTR: n=3; isoPD1: 3; PD1: 4; PD2: 4); D50 (3; 9; 9; 5); D80 (3; 9; 8; 5). cRepresentative raster plots (top) and global network activity (GNA, bottom) for CTR (SP11), isoPD1, and PD1 (SP12) neuronal cultures at D80. Each plot shows 5 min of recording. Peaks in the GNA reveal network bursts (blue dots). Extreme bursting events (red dots) are those that are above a threshold (red dotted line) set as 95% confidence interval of CTR bursts’distribution. CTR and isoPD1 networks show a relative low percentage of extreme events, which contrasts with the rich abundance of them in PD1 networks. dRatio of extreme events for all studied cell lines at D50 and D80. The colored boxes are a guide to visualize the distributions, which show the mean ± SD and the individual realizations (dots). For panels (b) and (d), the number of independent experiments for each condition and timepoint are D35 (CTR: n=3, isoPD1: 3, PD1: 4, PD2: 4); D50 (3, 9, 9, 5); D80 (3, 9, 8, 5). ***p< 0.001 (ANOVA with multiple comparison analysis). G. Carola et al. 4 npj Parkinson’s Disease (2021) 55 Published in partnership with the Parkinson’s Foundation Fig. 3 Functional connectivity of CTR, isoPD, and PD lines. a Top, Representative functional connectivity maps at D80 for CTR, isoPD1, and PD1 cultures. Dots are neurons and lines functional connections. The diameter of a dot and its opacity is proportional to the connectivity of the neuron it represents. Bottom, corresponding functional connectivity matrices. Black points are functional connections and colored boxes are functional communities. bComparison of the number of communities and community statistic Q for the three lines at D50 and D80. PD1 networks are excessively integrated, with a relatively small number of communities strongly interconnected (low Q) as compared to CTR and isoPD1 networks. The colored boxes are a guide to visualize the distributions, which show the mean±SD and the individual realizations (dots). Number of independent experiments: D50 (CTR: n=3, isoPD1: 4, PD1: 6); D80 (3, 5, 6). ****p< 0.001; **p< 0.01; *p< 0.05 (ANOVA with multiple comparison analysis). cSketch of functional differences between normal and PD networks according to the data, with PD networks displaying on average a lower connectivity combined with fewer and excessively connected communities. dCumulative distribution functions (CDF) of functional connections (in % of network) at D50 for CTR, PD1, PD2, and gene-edited isogenic PD1 line (isoPD1). The CDFs of the two PD lines exhibit a trend toward a lower connectivity as compared to CTR and isoPD1 ones, which are similar. eCorresponding distributions at D80. The two PD lines maintain their low connectivity profiles and strengthen their differences relative to CTR and isoPD1, which remain similar. The number of independent experiments used in panels (d) and (e) are D50 (CTR: n=3, isoPD1: 4, PD1: 6, PD2: 5); D80 (3, 5, 6, 5). G. Carola et al. 5 Published in partnership with the Parkinson’s Foundation npj Parkinson’s Disease (2021) 55 much larger (orange boxes), indicating not only a failed formation of functional microcircuits, but a tendency toward excessively strong network synchronicity, as also seen in the raster plots. These results were reproducible across realizations. As shown in Fig. 3b, PD networks exhibited a clear trend toward a small number of communities, which was accompanied by a tendency toward a stronger bond between these communities, i.e., an abnormal excessive integration. The latter was captured by the ‘community statistic Q’. The lower the Q value, the stronger the integration in the network. These differences appeared at D50 and substantially strengthened at D80. While CTR and isoPD cultures showed a moderate integration, PD cultures exhibited an abnormally strong integration. The combination of these results is summarized in Fig. 3c, which depicts two toy networks with the same number of neurons but with different functional organization. The normal, healthy network is characterized by a high average connectivity and well-defined small communities, while the PD network is characterized by a lower average connectivity and large, strongly linked communities that effectually almost shape a unique structure. To complete the functional analysis, we compared in more detail the statistics of functional connections. Figure 3d, e shows, for D50 and D80 stages of development, the cumulative distribution function of connections, CDF(k), for control, isoPD, and PD lines. This distribution portrays the probability that a neuron in the network has a number of connections less or equal than k. For D50, all distributions were similar and relatively close to one another, illustrating that they originated from similar dynamics, namely a combination of individual activity and network bursting. The distributions also showed a subtle growth, indicating that weakly connected neurons were rare. However, the PD distributions at D50 revealed a tendency toward a lower connectivity. Indeed, although the distributions were similar in shape, an analysis of the distance between distributions (Kullback–Leibler divergence and Kolmogorov–Smirnov test, Supp. Fig. 3) showed statistical differences among them, indicating that D50 could be the characteristic timepoint at which functional alterations in PD lines start to be detectable. The differences among distributions accentuated at D80, with PD cultures exhibiting a more pronounced trend toward a lower connectivity (a sharp increase of CDF(k) at low k values) that markedly departed from CTR and isoPD cultures. Interestingly, there were no statistically significant differences between CTR and isoPD cultures at this timepoint, a result that suggests the successful rescue of affected cell lines through correction of the LRRK2 mutation by gene edition. In light of these results, we considered the control network as the reference for a healthy development and function, and the departure from it as a signature of the pathology. Thus, we hypothesize that the LRRK2 mutation undermines the development of neuronal circuitry to such an extent that it alters collective network activity and functional organization. Contrasting dynamics of TH and non-TH neurons in control and PD lines To investigate the origin of the functional impairment found in PD DAn, we took advantage of the genetic TH-reporter tool created in our lab 19 . This reporter allows identifying TH and non-TH neuronal populations in the networks and analyzing their functional characteristics separately. Figure 4a, b exemplifies such a construction for representative isoPD and PD networks at D80, in which two neuronal layers, one corresponding to TH neurons (red) and another one corresponding to non-TH neurons (blue), interact functionally. The spontaneous activity of each subpopulation of neurons at D80 is shown in the accompanying raster plots. A comparison of the activity patterns between isoPD and PD cultures reveals the strong contrast in their dynamics. While the non-TH population in isoPD cultures shows a sustained activity with weak collective events, PD cultures show strong synchrony episodes that extend to both subpopulations. In addition, as shown in Fig. 4c, a comparison of the ratio of extreme events for TH and non-TH subpopulations indicates that, on average, PD cultures exhibited a higher number of extreme events in both subpopulations when compared to controls. However, along development, TH and non-TH ratios for controls were similar at D50 and D80, whereas the ratios for the PD line switched from mostly TH at D50 to mostly non-TH at D80. Altogether, these results unfold an abnormal subpopulation dynamics in PD networks, with an overall excessive bursting and a reversing leadership of TH and non-TH subpopulations along development. We also noted that PD cultures at D50 showed important variability among realizations, with extreme events in the non-TH population being absent in some experiments and abundant in others. Given this variability, we hypothesize that the D50 timepoint might mark the onset of structural alterations that later translate into dynamic and functional deficits. Since about 70% of the neurons in our networks are non-TH (Fig. 1g, h), the above results suggest that, in healthy control cultures, non-TH neurons drive spontaneous activity in the network. In contrast, TH neurons play a regulatory yet essential role by facilitating the coexistence of small neuronal coactivations and whole-network bursting events. This regulatory role appears to be lost in PD cultures, in which TH neurons frame spontaneous activity patterns with excessive synchrony that translate into an abundance of extreme bursting episodes. We, therefore, hypothesize that the LRRK2 mutation alters the balance between neuronal subpopulations by degrading the physical coupling among neuronal types. An in silico model captures the dynamic alterations in PD networks Next, we carried out numerical simulations to better understand the impact of TH+cells structural failure on PD network dynamics. Specifically, we evaluated whether a reduction in TH cells connectivity was sufficient to switch global network dynamics from balanced to excessively synchronous. We first reproduced the behavior of control networks. As shown in Fig. 5a, we used the same neuronal spatial arrangement of the cultured networks and considered a mixed population of 55% excitatory non-TH neurons, 25% excitatory TH+neurons, and 20% inhibitory neurons. These values were selected to concord with both the typical 80% excitation of cortical circuits in vitro 24 and the TH/DAPI ratio found in our cultures (Fig. 1). Neuronal dendritic trees and axons were then incorporated according to realistic biological rules, so that a connection was established whenever an axon crossed the dendritic tree of a neuron (Fig. 5b, top). Once the network structure was set, dynamics was incorporated through an extended Hodgkin–Huxley model with parameters adjusted to replicate the activity in controls. As shown in the top panels of Fig. 5c, normal network dynamics was qualitatively similar to the experiments shown in Fig. 2c (left and center) and characterized by an intense background activity in combination with coordinated activity episodes. The inspection of the different populations (Fig. 5c, top, right panels) shows that non-TH neurons were also the drivers of network dynamics. It should be noted that inhibition was necessary in the simulations to ensure a sufficiently high spontaneous activity. To explore PD dysfunction, we used the same control networks and modeled the shortening of neurites in TH neurons, effectually reducing the connectivity probability in the TH population (Fig. 5b, bottom). Such a construction signified that a single TH cell lost about 80% of its connections. We then explored the minimum fraction of TH cells that had to be affected to observe an impact on the network dynamic. As shown in the bottom panels of Fig. G. Carola et al. 6 npj Parkinson’s Disease (2021) 55 Published in partnership with the Parkinson’s Foundation 5c, PD simulations were characterized by an exceedingly synchronous behavior that captured well the experimental observations of PD cultures (Fig. 2c, right), with a much lower background activity and a similar dynamic in both TH and non-TH populations. Simulations also demonstrated that the affectation of ~10% TH cells sufficed to drive the networks toward a chronic bursting behavior with an abundance of extreme, whole-network synchronous events (Fig. 5c). As shown in Fig. 5d, we compared three different ratios of affected TH neurons. An increase in this ratio mimics disease progression, i.e., developmental time in the experiments. The results show that the ratio and occurrence of extreme events are similar for different damage rates, possibly because the remaining excitatory and inhibitory neurons are sufficiently large to maintain activity. Interestingly, these results also indicate that a small damage suffices to drive the network toward an excessively synchronous scenario. Finally, an analysis of functional connectivity showed that the model also reproduces well the experimental observations, with a tendency for simulated PD cultures to shift toward lower connectivity states (Fig. 5e). TH+neurons carrying LRRK2 mutation feature a lower number of dendrites To verify the prediction of the in silico model above and expose physical alterations in the circuitry of PD lines, we investigated the morphology of TH neurons in CTR, isoPD, and PD neuronal cultures at D50 and D80 (Fig. 6a and Supplementary Fig. 4a). We found that DAn differentiated from PD iPSC lines showed a lower number of TH neurites compared to those derived from CTR or isoPD lines (1.2 ± 0.1 neurites for PD vs. 4.5 ± 0.2 for CTR and 4.1 ± 0.3 for isoPD, Fig. 6b). In addition, the number of neurites in CTR and isoPD cultures increased along with development while it decreased in PD cultures, indicating a progressive deterioration of network structure for the latter (Fig. 6b). In contrast, the number of neurites in MAP2+neurons from control, isoPD, and PD lines did not show any significant differences (Fig. 6c), confirming that the dynamic and functional deficits of PD lines are localized to the TH subpopulation, and that these cells experience a gradual structural failure in the form of neurite loss. Next, we evaluated DA release in the media collected from CTR, isoPD, and PD neuronal cultures. Supernatants of PD cultures revealed decreased dopamine levels at D50 and D80 compared with those of CTR and isoPD cultures (Fig. 6d), indicating that the reduction in TH+neurites has functional consequences on DA release. Finally, to investigate whether the early sign of overall alteration in the network structure and dynamics could contribute to neurodegeneration, we maintained iPSC-derived DAn over culture times longer than D80. Interestingly, in aged (110 days, latest timepoint analyzed) cultures, PD DAn showed morphological alterations, including reduced number and length of neurites, and significantly decreased cell survival compared with isoPD DAn (Supplementary Fig. 5a–d). In contrast, neuronal degeneration was Fig. 4 Functional connectivity and dynamics of TH+and non-TH+neurons. a Left, representative functional network of an isoPD1 culture at D80, signaling the location of TH+neurons (red) and non-TH+neurons (blue). The diameter of a dot is proportional to the connectivity of the neuron it represents. Only 10% of the connections are shown for clarity. Right, corresponding raster plots of the two populations. bCorresponding analysis for a PD1 culture, in which the non-TH+subpopulation exhibits a much stronger synchronous behavior as compared to controls. cRatio of extreme events for each subpopulation of neurons at two maturation stages, D50 and D80. PD cultures show in general a higher ratio of extreme events as compared to controls. The non-TH+population in the PD network at D50 shows a strong variability in the ratio of extreme events across realizations, indicating the onset of malfunction. The same population at D80 is dominated by extreme events that reflect the strong synchronous behavior. The colored boxes are a guide to visualize the distributions, which show the mean ± SD and the individual realizations (dots). The number of independent experiments in panel (c) are D50 (CTR: n=3, isoPD1: 5, PD1: 5); D80 (3, 5, 5). ***p< 0.001; **p< 0.01; *p< 0.05 (ANOVA with multiple comparison analysis). G. Carola et al. 7 Published in partnership with the Parkinson’s Foundation npj Parkinson’s Disease (2021) 55 not evident in non-TH+cells, as judged by the percentage of MAP2+/TH−neurons in the neuronal cultures (Supplementary Fig. 5e). DISCUSSION Since early diagnosis of PD is expected to dramatically improve the outcome of therapies under current development, in the present study we interrogated a human neuronal cell-based model of PD for the earliest detectable functional alterations. We found that DA neurons derived from iPSC representing healthy individuals or PD patients harboring LRRK2 mutation developed appropriate physiological characteristics forming complex and mature networks during the differentiation process. However, PD neuronal networks developed abnormal hypersynchrony at the latest timepoint analyzed (D80), in contrast with healthy or geneedited isogenic PD networks. These new data from PD-affected human DAn indicate that early dysfunction may contribute to the initiation of downstream degenerative pathways that ultimately lead to DAn loss in PD. A general limitation of human iPSC-based disease modeling strategies that should be taken into account when interpreting the results of our studies is the notorious Fig. 5 Numerical simulations of normal and PD networks. a Representative CTR culture at day 80, used as a reference for the spatial allocation of neurons in the numerical model. The enlarged area depicts a detail of the neuronal spatial arrangement. 300 neurons are used in the simulation, and are randomly assigned as dopaminergic (DA, red, 25% of network), inhibitory (green, 20%), and excitatory (blue, 55%). bRandom pruning algorithm, in which PD damage is simulated by a shortening of either axons (red curve) or dendritic arbor (dotted red circle) in DA neurons (red hexagons). In the sketch, normal connectivity is established when the dendritic arbor of any neuron (dotted circles) intersects an axon. Neurite pruning disconnects neurons either because of axonal shortening (neuron #1) or dendritic loss (neuron #2) in DAn. cRaster plots of simulated normal and PD networks, with the latter corresponding to an axonal pruning on 10% of DA neurons, each one losing about 80% of connections. In the raster plots, the left panels show the dynamics of the entire network, whereas the right panels show the dynamics in each subpopulation. All raster plots range from 0 to 300, but the number of active neurons in each plot vary according to the population monitored. PD simulations show a markedly synchronous behavior of the excitatory population that translates onto the entire network. dRatio of extreme events for different percentages of pruned DAn subpopulation, showing that the presence of extreme bursting events increases as network affectation progresses. The colored boxes are a guide to visualize the distributions, which show the mean ± SD and the individual realizations (dots). The number of individual simulations is n=4 for each condition. ***p< 0.001; *p< 0.05 (ANOVA with multiple comparison analysis). eCumulative distribution of functional connections (CDF) for two PD networks (at 10 and 30% pruning) and a normal, non-pruned network. The distributions show the same trend as in the experiments, with PD departing from normal toward a low degree functional connectivity scenario. G. Carola et al. 8 npj Parkinson’s Disease (2021) 55 Published in partnership with the Parkinson’s Foundation variability described among iPSC lines and clones 25 . To rule out the impact of interline/interclone variability in our findings, we used several iPSC lines and iPSC clones to represent each experimental condition. Specifically, controls included two independent iPSC lines generated from two independent individuals, and a total of four iPSC clones, whereas the PD condition was similarly represented by two independent iPSC lines generated from two independent individuals, and a total of three iPSC clones (Table 1). Moreover, the confirmation of our findings in neuronal cultures differentiated from gene-edited isogenic controls of PD iPSC provides further reassurance that the alterations identified in PD samples in our studies are indeed related to the disease condition, rather than reflecting the specific iPSC line/clone utilized. Our previous work showed that PD patient-specific iPSC-derived DAn are particularly susceptible to undergo neurodegeneration Fig. 6 Biological confirmation of the in silico hypothesis. a Gene-edited isoPD1 (SP12 wt/wt) and PD1 (SP12) TH+neuronal processes are traced (orange) to determine the number and the structure of neurites at D50 (left panels) and D80 (right panels). The top panels for each cell line show the original image of TH+neurons. A blue triangle highlights a representative traced neuron shown in the middle panels. The bottom panels for each cell line show MAP2+neuronal processes. Scale bar is 50 μm in all images. bThe quantification of the number of TH neurites shows significant differences between CTR and PD1 both at D50 and D80. cNo differences are detected on the number of MAP2 neurites. dLevels of released DA measured in culture media from CTR (SP11), PD1 (SP12), PD2 (SP13), and isoPD1 (SP12 wt/wt) at D50 and D80 of differentiation, relative to CTR levels. *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001 (ANOVA with multiple comparison analysis). Number of independent experiments n=3. G. Carola et al. 9 Published in partnership with the Parkinson’s Foundation npj Parkinson’s Disease (2021) 55